一種顧及永久散射體的SBAS InSAR時(shí)序地表沉降提取方法

王舜瑤，盧小平，劉曉幫，付睢寧

(河南理工大學(xué)礦山空間信息技術(shù)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室，河南 焦作 454003)

地表沉降是制約城市規(guī)劃與發(fā)展的重要因素，常規(guī)的地表沉降手段如水準(zhǔn)測(cè)量、GPS測(cè)量等方法，存在監(jiān)測(cè)點(diǎn)過(guò)于稀疏、效率低等不足。而新興的合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(differential interferometric synthetic aperture radar,D-InSAR)能夠大范圍、無(wú)接觸、全天候、全天時(shí)地獲取地表的毫米級(jí)形變信息[1-4]，尤其是永久散射體(persist scatterers InSAR,PS InSAR)和小基線(xiàn)集(small BAseline subsets InSAR,SBAS InSAR)技術(shù)[5-7]，前者能有效提取非線(xiàn)性形變并消除大氣相位影響，后者可以顯著增加相干影像對(duì)的數(shù)目，解決了時(shí)間、空間失相關(guān)和大氣效應(yīng)對(duì)InSAR的影響；同時(shí)，隨著高分辨率SAR衛(wèi)星的成功運(yùn)行，D-InSAR技術(shù)在地表沉降監(jiān)測(cè)中發(fā)揮著越來(lái)越重要的作用[8]。文獻(xiàn)[9]利用PS InSAR對(duì)香港國(guó)際機(jī)場(chǎng)沉降的提取，監(jiān)測(cè)結(jié)果與水準(zhǔn)測(cè)量一致。文獻(xiàn)[10]借助PS InSAR及高分辨率X波段Cosmo-SkyMed雷達(dá)影像，成功獲取了上海市地鐵隧道沿線(xiàn)及部分高速公路的沉降情況。文獻(xiàn)[11]利用SBAS InSAR技術(shù)并結(jié)合GPS和地下水水位數(shù)據(jù)，得到了洛杉磯地表時(shí)序形變及含水層貯水系數(shù)等物理參數(shù)。文獻(xiàn)[12]利用SBAS InSAR技術(shù)監(jiān)測(cè)凍土形變，所得到的結(jié)果與物理變化規(guī)律相吻合。PS/SBAS InSAR各有優(yōu)點(diǎn)，這兩種技術(shù)的交叉融合也是當(dāng)前D-InSAR研究的熱點(diǎn)。文獻(xiàn)[13]結(jié)合SBAS和PS的相干點(diǎn)選取方法，提高了相干點(diǎn)的數(shù)量，并提取了某火山的形變信息。文獻(xiàn)[14]融合PS/SBAS InSAR算法，通過(guò)構(gòu)建局部Delaunay三角網(wǎng)，從某一參考點(diǎn)開(kāi)始對(duì)形變速率增量作積分，提取了徐州大區(qū)域地表形變。由于在SBAS InSAR的軌道精煉步驟中，引入穩(wěn)定的地面控制點(diǎn)可以有效估計(jì)殘留相位并去除平地效應(yīng)，從而提高結(jié)果的可靠性。然而，在沒(méi)有地面控制點(diǎn)信息的情況下，盲目地選擇地面控制點(diǎn)反而可能引入誤差。

針對(duì)該問(wèn)題，本文提出一種顧及永久散射體的SBAS InSAR處理方法，通過(guò)相干性、振幅離差指數(shù)、形變速率3閾值提取穩(wěn)定的永久散射體，并將其作為地面控制點(diǎn)引入SBAS InSAR處理流程，最后通過(guò)實(shí)例驗(yàn)證了本文方法的有效性。

1 本文方法

1.1 基于振幅離差和相干性雙閾值的永久散射體探測(cè)方法

對(duì)于復(fù)數(shù)SAR影像，若其實(shí)部和虛部分別存在標(biāo)準(zhǔn)差為σs的高斯噪聲，則振幅值A(chǔ)服從Rice分布，即

(1)

式中，I0為Bessel函數(shù)；g為目標(biāo)反射量，是正實(shí)數(shù)。當(dāng)比值g/σs較小時(shí)，即該目標(biāo)為低信噪比目標(biāo)時(shí)，Rice分布演化為Rayleigh分布。而對(duì)于高信噪比目標(biāo)(通常g/σs>4)，該分布接近于高斯分布。當(dāng)σs?g時(shí)，存在以下關(guān)系

(2)

式中，σp為相位標(biāo)準(zhǔn)差；σsI為虛部標(biāo)準(zhǔn)差；σA為時(shí)序振幅標(biāo)準(zhǔn)差；mA為時(shí)序振幅均值；DA為振幅離差指數(shù)。

相干系數(shù)的數(shù)學(xué)表達(dá)式為

(3)

式中,γ為相干系數(shù)；M、S分別為構(gòu)成干涉對(duì)的SAR影像；*為復(fù)數(shù)共軛；m、n分別為局部窗口的大小；i、j則為像素坐標(biāo)。

設(shè)置振幅離差閾值TD對(duì)永久散射體目標(biāo)進(jìn)行初次探測(cè)，而后設(shè)置相干性閾值Tγ對(duì)永久散射體目標(biāo)進(jìn)行二次探測(cè)，濾除大部分失相干嚴(yán)重的像素點(diǎn)，從而獲得最終探測(cè)結(jié)果。

1.2 基于SBAS InSAR的時(shí)序地表沉陷提取

在(t0，t1，…，tn)時(shí)間內(nèi)，從獲取的同一地區(qū)N景SAR影像中選取1景圖像作為主圖像，將其他SAR圖像與主圖像配準(zhǔn)。假設(shè)每期SAR圖像都滿(mǎn)足臨界基線(xiàn)要求，則可生成M組干涉像對(duì)，并滿(mǎn)足下列條件

N/2≤M≤N(N-1)/2

(4)

假設(shè)前后兩景影像獲取時(shí)間分別為tA和tB，這兩景影像差分生成第j景干涉條紋圖，在不考慮大氣延遲相位、殘余地形相位和噪聲相位的情況下，干涉相位可表示為

(5)

式中，φ為干涉相位；j為差分干涉圖的景號(hào)，j∈(1,…,M)；λ為電磁波波長(zhǎng)；d(tA,x,r)和d(tB,x,r)分別為tA和tB時(shí)刻相對(duì)于d(t0,x,r)=0的雷達(dá)視線(xiàn)方向累積形變量。

相位信息是通過(guò)相位解纏獲得，已知像素點(diǎn)的形變量對(duì)應(yīng)的未知相位值為

φ=[φ(t0)φ(t1)…φ(tN-1)]T

(6)

解纏后的相位值為

δφdef=[δφdef(Δt0)δφdef(Δt1)…δφ(ΔtM-1)]T

(7)

將形變累積相位與干涉相位序列的線(xiàn)性關(guān)系用矩陣表示為

δφdef=Aφ

(8)

式中，A為M×N階系數(shù)矩陣。若M≥N，通過(guò)最小二乘法可得

φ=(ATA)-1ATδφdef

(9)

若M

vj=(φj-φj-1)/(tj-tj-1)

(10)

用式(10)替換式(5)中的相位，可得

(11)

各時(shí)段速度在主、從圖像時(shí)間間隔上的積分用矩陣可表示為

δφdef=Bv

(12)

式中,B為M×N階矩陣。第j行位于主輔影像獲取時(shí)間之間的列，B(j,k)=tk-tk-1，其他情況B(j,k)=0。此時(shí)，將奇異值分解應(yīng)用于矩陣B，就可得到速度矢量v的最小范數(shù)解，獲得相位的線(xiàn)性形變信息。對(duì)殘余相位的時(shí)間和空間進(jìn)行濾波，可分離出大氣相位、非線(xiàn)性形變相位，將這兩項(xiàng)形變疊加即可求出精確的地表形變量。

2 試驗(yàn)及處理流程

試驗(yàn)數(shù)據(jù)為21景Sentinel 1A衛(wèi)星C波段單視復(fù)數(shù)影像，成像時(shí)間為2016年11月至2017年11月，覆蓋鄭州全市域。本文選擇鄭州市四環(huán)路以?xún)?nèi)作為研究區(qū)域，并對(duì)數(shù)據(jù)進(jìn)行裁剪，試驗(yàn)區(qū)域如圖1所示。區(qū)域左上角和右下角的地理坐標(biāo)分別為(34.847 5°N，113.565 1°E)、(34.659 8°N，113.790 8°E)。區(qū)域內(nèi)有3個(gè)三等水準(zhǔn)點(diǎn)SZ01、SZ02和SZ03，可作為精度驗(yàn)證的地面控制點(diǎn)。用于去平地處理的DEM數(shù)據(jù)為美國(guó)航空航天局(NASA)發(fā)布的30 m分辨率SRTM(shuttle radar topography mission，SRTM)V4版本[15]。精密衛(wèi)星軌道數(shù)據(jù)由Sentinel-1 Quality Control網(wǎng)站提供。試驗(yàn)流程如圖2所示。

試驗(yàn)處理軟件采用瑞士sarmap公司研發(fā)的SARScape，該軟件架構(gòu)于專(zhuān)業(yè)的ENVI遙感圖像處理軟件之上，提供用戶(hù)友好的圖形化操作界面，具有專(zhuān)業(yè)雷達(dá)圖像處理和分析功能[16]，能夠滿(mǎn)足本次試驗(yàn)的需求。

2.1 地面控制點(diǎn)的提取

首先對(duì)原始SLC數(shù)據(jù)進(jìn)行多視處理，得到強(qiáng)度信息圖，并根據(jù)試驗(yàn)區(qū)域?qū)?qiáng)度信息圖進(jìn)行裁剪，去掉不需要的區(qū)域，以提高后續(xù)處理的效率。對(duì)裁剪后的數(shù)據(jù)進(jìn)行PS點(diǎn)提取，其中探測(cè)PS點(diǎn)的振幅離差閾值TD設(shè)置為3.2。由于研究區(qū)為城市區(qū)域，存在大量永久散射體(如大型建筑物、路燈、鐵路、公路等)，相干性閾值設(shè)置過(guò)低會(huì)導(dǎo)致PS點(diǎn)數(shù)目劇增，大大增加數(shù)據(jù)計(jì)算的時(shí)間，從而增加后續(xù)篩選工作的難度，通過(guò)多次試驗(yàn)最終確定相干性閾值Tγ為0.97，共探測(cè)到3580個(gè)PS點(diǎn)，如圖3所示。通過(guò)空間域低通、時(shí)間域高通濾波的方法減小干涉相位中的大氣延遲相位誤差，獲取更加可靠的PS點(diǎn)形變信息，將形變速率閾值Tv設(shè)為±0.1 mm/a，再次對(duì)這些PS點(diǎn)進(jìn)行篩選，獲得33個(gè)3閾值探測(cè)PS點(diǎn)，并將這些點(diǎn)作為地面控制點(diǎn)引入后續(xù)的SBAS InSAR處理流程。由于SBAS InSAR需要的地面控制點(diǎn)為雷達(dá)斜距坐標(biāo)系，因此需要對(duì)3次PS探測(cè)結(jié)果進(jìn)行投影變換，將地理坐標(biāo)系下的PS點(diǎn)轉(zhuǎn)換為雷達(dá)斜距坐標(biāo)備用，如圖4所示。

2.2 改進(jìn)的SBAS InSAR處理方法

對(duì)裁剪后的數(shù)據(jù)進(jìn)行SBAS InSAR處理，空間基線(xiàn)、時(shí)間基線(xiàn)分別如圖5、圖6所示，其中空間基線(xiàn)閾值為5%，時(shí)間基線(xiàn)閾值為100 d。

對(duì)干涉處理后得到的干涉圖快視圖進(jìn)行目視解譯，手動(dòng)去除相干性較差的干涉相對(duì)，避免對(duì)形變提取結(jié)果產(chǎn)生影響。在軌道精煉和重去平步驟中，導(dǎo)入之前通過(guò)3次閾值篩選生成的地面控制點(diǎn)，代替人工目視選取地面控制點(diǎn)，再通過(guò)兩次反演去除殘余地形和大氣相位，經(jīng)過(guò)地理編碼后獲得最終形變速率和沉降時(shí)間序列圖。

3 試驗(yàn)結(jié)果及說(shuō)明

試驗(yàn)獲得了測(cè)區(qū)的形變速率圖(如圖7所示)及相對(duì)于第一期SAR圖像(2017年10月26號(hào))的沉降量時(shí)間序列圖(如圖8所示)。試驗(yàn)結(jié)果表明，該地區(qū)出現(xiàn)區(qū)域性形變，形變速率從-61～27 mm/a不等，市中心區(qū)域出現(xiàn)毫米級(jí)小幅抬升，南部抬升程度較大，最大值為49 mm，而西南部則呈下降趨勢(shì)，最大下沉量為-82 mm，西北部和東部也出現(xiàn)不規(guī)則下沉。

為驗(yàn)證本文方法的優(yōu)越性，使用相同原始數(shù)據(jù)，通過(guò)SBAS InSAR及PS InSAR方法分別獲取測(cè)區(qū)時(shí)序地表沉降，并將結(jié)果與本文方法所獲結(jié)果進(jìn)行對(duì)比，同時(shí)引入3個(gè)三等水準(zhǔn)測(cè)量點(diǎn)位的兩次測(cè)量結(jié)果進(jìn)行驗(yàn)證，具體數(shù)據(jù)情況見(jiàn)表1。假設(shè)水準(zhǔn)實(shí)測(cè)數(shù)據(jù)為真值，并將兩次水準(zhǔn)測(cè)量結(jié)果與時(shí)間最接近的SAR影像時(shí)間段內(nèi)(2017年3月31日—2017年9月15日)的沉降情況進(jìn)行對(duì)比，結(jié)果見(jiàn)表2，其中正值為抬升，負(fù)值為下降。

表1 三等水準(zhǔn)點(diǎn)高程測(cè)量結(jié)果及形變量統(tǒng)計(jì)

表2 本文方法與SBAS InSAR、PS InSAR和水準(zhǔn)測(cè)量的形變量對(duì)比 mm

由表1和表2可知：①PS InSAR和SBAS InSAR提取的時(shí)序地表形變精度均達(dá)到亞厘米級(jí)，在該試驗(yàn)區(qū)內(nèi)，PS InSAR結(jié)果比SBAS InSAR結(jié)果更接近真值，這是由于城市中的大型建筑物、樓房、高架橋、路燈等可作為永久散射體的情況較多，能充分發(fā)揮PS InSAR的優(yōu)勢(shì)；②本文方法提取結(jié)果比SBAS InSAR和PS InSAR提取結(jié)果更接近真值。

但本文方法與水準(zhǔn)測(cè)量產(chǎn)生了一定的誤差，可能的原因如下：①SAR數(shù)據(jù)分辨率的限制，數(shù)據(jù)的單個(gè)像元實(shí)地并不是一個(gè)點(diǎn)，覆蓋了一定的范圍，無(wú)法精確與水準(zhǔn)點(diǎn)重合，其沉降情況受到該區(qū)域整體沉降情況的綜合影響；②兩次水準(zhǔn)測(cè)量的時(shí)間與SAR數(shù)據(jù)的獲取時(shí)間存在一定差異，無(wú)法精準(zhǔn)重合；③水準(zhǔn)測(cè)量本身可能存在誤差，而且時(shí)序數(shù)據(jù)量不足，無(wú)法充分驗(yàn)證地表時(shí)序沉降量，若有時(shí)序GPS高程點(diǎn)測(cè)量結(jié)果支撐驗(yàn)證，則更能體現(xiàn)InSAR方法的時(shí)效性。

4 結(jié) 語(yǔ)

為彌補(bǔ)PS InSAR、SBSA InSAR單一處理方法的不足，本文提出了一種顧及永久散射體的SBAS InSAR時(shí)序地表形變提取方法，并通過(guò)鄭州市區(qū)的實(shí)例研究得出以下結(jié)論：①PS InSAR和SBAS InSAR在城市區(qū)域的時(shí)序地表形變提取精度均可達(dá)到亞厘米級(jí)，而且由于城市中存在較多永久散射體，因此PS InSAR比SBAS InSAR更適用；②本文方法提取結(jié)果優(yōu)于SBAS InSAR和PS InSAR，可以更有效地提取地表時(shí)序形變。