CRInSAR與PSInSAR融合探測(cè)地表線性變形速率

邢學(xué)敏 賀躍光 聞德保 周訪濱

1 長(zhǎng)沙理工大學(xué)特殊環(huán)境道路工程湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙市萬(wàn)家麗南路960號(hào)，410014

2 長(zhǎng)沙理工大學(xué)交通運(yùn)輸工程學(xué)院，長(zhǎng)沙市萬(wàn)家麗南路960號(hào)，410014

3 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙市麓山南路932號(hào)，410083

永久散射體差分干涉測(cè)量技術(shù)（PSInSAR）廣泛應(yīng)用于區(qū)域地表形變的探測(cè)中［1－3］。然而，受區(qū)域相干性的限制，可能導(dǎo)致一些低相干區(qū)域內(nèi)探測(cè)不到足夠數(shù)目的PS點(diǎn)。人工角反射器技術(shù)（CR）是通過(guò)在研究區(qū)域內(nèi)安裝人工角反射器，通過(guò)對(duì)這些點(diǎn)上的相位進(jìn)行建模分析，以求解其地表變形速率。人工角反射器安裝靈活，當(dāng)研究區(qū)域內(nèi)探測(cè)不到PS點(diǎn)或當(dāng)PS點(diǎn)密度不足時(shí)可發(fā)揮較強(qiáng)的優(yōu)勢(shì)［4］。此外，在CR 點(diǎn)上的GPS接收機(jī)還可獲取CR 點(diǎn)的高程及坐標(biāo)信息，一方面可大大避免在差分干涉處理過(guò)程中由于地理編碼而引入的誤差，另一方面可作為PS 網(wǎng)絡(luò)解纏的起算數(shù)據(jù)，避免由于人為選取假設(shè)穩(wěn)定點(diǎn)而引入的不確定性。本文將CRInSAR 技術(shù)與PSInSAR融合起來(lái)用于解算地表線性變形速率。設(shè)計(jì)了詳細(xì)的算法流程，分別利用模擬實(shí)驗(yàn)及真實(shí)數(shù)據(jù)實(shí)驗(yàn)對(duì)算法進(jìn)行驗(yàn)證。

1 PSInSAR 算法

假設(shè)有N＋1幅時(shí)間序列影像，利用二軌差分干涉處理方法可對(duì)應(yīng)生成N幅差分干涉圖。選取一定數(shù)量的PS點(diǎn)后，將任意相鄰的兩個(gè)PS點(diǎn)進(jìn)行連接，可組建成PS 基線網(wǎng)絡(luò)。針對(duì)網(wǎng)絡(luò)中任意一條PS基線（假設(shè)包含第i和第j個(gè)PS點(diǎn)），其相位差可以表示為［5－6］：

其中m為干涉對(duì)序號(hào)；分別為兩PS點(diǎn)的線性變形速率增量和高程改正增量；分別為相鄰兩PS點(diǎn)在第m幅干涉圖中對(duì)應(yīng)的相位增量和整周模糊度增量；為高程相位改正系數(shù)，其中分別表示時(shí)空基線，θ為影像的入射角，Ri表示雷達(dá)傳感器與PS點(diǎn)之間的距離；表示殘余相位，主要由大氣延遲相位、非線性形變相位和噪聲組成［7］。

針對(duì)上式未知參數(shù)Δvi，j、ΔδHi，j及的求解可利用LAMBDA算法進(jìn)行［8］，ΔδHi，j及Δvi，j的估值結(jié)果可作為空間維相位解纏的觀測(cè)數(shù)據(jù)，利用間接平差的方法求解出各PS 點(diǎn)上的絕對(duì)線性變形速率v及高程改正值δH。

2 CRInSAR 算法

利用CRInSAR 技術(shù)進(jìn)行形變監(jiān)測(cè)時(shí)，需要在影像上先提取出CR 點(diǎn)的行列號(hào)信息，然后建立如下函數(shù)關(guān)系模型［9］：

式中p表示CR 點(diǎn)序號(hào)，m表示干涉對(duì)序號(hào)，標(biāo)有p和m的各分量均表示第p個(gè)CR 點(diǎn)與參考CR點(diǎn)的對(duì)應(yīng)參數(shù)差；為干涉相位差；為整周模糊度差；為地形引起的相位差；vp為相對(duì)于參考點(diǎn)的形變速率差；Tm為時(shí)間基線；是殘余相位，與（1）式中ε含義相同；地形相位，其中λ是雷達(dá)波長(zhǎng)，Bm為垂直基線長(zhǎng)度，Rp為CR 點(diǎn)與衛(wèi)星位置間距離，ΔHp是CR 點(diǎn)與參考點(diǎn)的相對(duì)高程。

實(shí)際計(jì)算中，假設(shè)有N＋1 幅SAR 影像，有M＋1個(gè)CR 點(diǎn)，選取一個(gè)CR 點(diǎn)為參考點(diǎn)，一幅SAR 影像為主影像，可對(duì)應(yīng)生成N幅干涉圖。當(dāng)CR 點(diǎn)間距在1 000 m 范圍內(nèi)時(shí)，大氣相位影響可忽略，因此，可不考慮殘余相位影響［10］。而均為已知量，則在式（2）中未知參數(shù)僅為整周模糊度增量Δk、形變速率v。利用LAMBDA 算法，可逐步分離出未知參數(shù)，進(jìn)而得出最終的CR 點(diǎn)上的形變速率值vp。

3 PSInSAR 與CRInSAR 融合解算算法

利用LAMBDA 算法獲取了相鄰PS點(diǎn)的線性變形速率增量及高程改正增量后，即可作為輸入數(shù)據(jù)，用于實(shí)現(xiàn)任意PS 點(diǎn)上絕對(duì)線性變形速率及高程改正值的求解，這一過(guò)程又被稱(chēng)為空間維相位解纏。傳統(tǒng)空間維解纏會(huì)受主觀選取穩(wěn)定點(diǎn)的制約，帶有很大的不確定性［11－12］。我們預(yù)先在研究區(qū)域內(nèi)安裝CR 點(diǎn)，獲取其線性變形速率和高程改正值作為PS基線網(wǎng)絡(luò)的空間維解纏的約束數(shù)據(jù)，利用間接平差方法求解各PS 點(diǎn)的絕對(duì)線性變形速率值和高程改正值，進(jìn)而實(shí)現(xiàn)PS基線網(wǎng)絡(luò)的空間維解纏。圖1為融合算法中CR點(diǎn)與PS點(diǎn)相對(duì)位置示意圖，以線性變形速率為例，針對(duì)任一條基線邊，我們都可以將Δvi，j作為間接平差函數(shù)模型的觀測(cè)值，線性變形速率參數(shù)vi、vj則為待求未知參數(shù)。間接平差函數(shù)模型可表示為：

圖1 PS與CR 點(diǎn)相對(duì)位置示意圖Fig.1 Diagrammatic sketch of reference position between CRs and PSs

上式中，Hcr為CR 點(diǎn)上的線性變形速率，可通過(guò)§2介紹的算法預(yù)先計(jì)算得出；δHcr為CR 點(diǎn)上的高程改正結(jié)果，可以通過(guò)式δHcr＝Hcr－Hdem進(jìn)行計(jì)算。其中Hcr表示CR 點(diǎn)上的GPS接收機(jī)獲取的高程值，Hdem表示CR 點(diǎn)對(duì)應(yīng)的外部DEM 高程結(jié)果。由此，可利用最小二乘原理將研究區(qū)域內(nèi)所有PS點(diǎn)上的線性變形速率及高程改正參數(shù)求解出來(lái)。根據(jù)上述討論，將PSIn－SAR與CRInSAR 融合解算算法流程圖表示在圖2中。

圖2 融合算法流程Fig.2 Flowchart of combined calculation algorithm

4 模擬實(shí)驗(yàn)

4.1 實(shí)驗(yàn)流程

為驗(yàn)證前述算法流程，設(shè)計(jì)一套模擬實(shí)驗(yàn)。首先利用peaks函數(shù)來(lái)模擬線性變形速度場(chǎng)，如圖3 所示。模擬速度場(chǎng)的總面積為100 像素×100像素，區(qū)域內(nèi)以無(wú)形變?yōu)橹鳎冃螀^(qū)域有下沉、抬升，呈對(duì)稱(chēng)分布。從中隨機(jī)選取200個(gè)像素作為參考PS 點(diǎn)，12 個(gè)像素作為CR 點(diǎn)。為增大PS點(diǎn)密度，再隨機(jī)選取600個(gè)點(diǎn)位作為擴(kuò)展PS點(diǎn)。對(duì)參考PS點(diǎn)和擴(kuò)展PS點(diǎn)共同建立起分級(jí)PS基線網(wǎng)，類(lèi)似于大地測(cè)量中分級(jí)控制網(wǎng)的方式（圖4）。參考PS點(diǎn)類(lèi)似于大地測(cè)量控制網(wǎng)中一級(jí)控制點(diǎn)，而擴(kuò)展PS點(diǎn)為二級(jí)控制點(diǎn)，需要將參考PS點(diǎn)與其鄰近的擴(kuò)展PS點(diǎn)相連接構(gòu)成基線邊。高程改正值采用高斯模型進(jìn)行模擬，其均值為0m，標(biāo)準(zhǔn)偏差為±3m。相位整周模糊度利用隨機(jī)整數(shù)模擬器模擬。這樣，PS點(diǎn)與CR 點(diǎn)的模擬真實(shí)速率值和模擬真實(shí)高程改正值均已給出。

圖3 模擬線性變形速率場(chǎng)Fig.3 The simulated linear velocity field

圖4 參考PS點(diǎn)與擴(kuò)展PS點(diǎn)位置示意圖Fig.4 Distribution of reference PSs and extended PSs

由式（1），差分干涉相位可以根據(jù)上面介紹的各模擬分量組建而成，其中涉及到的影像參數(shù)可直接在SAR 影像的頭文件中讀取（如時(shí)間基線、空間基線、高程相位轉(zhuǎn)換系數(shù)）。由于現(xiàn)有數(shù)據(jù)主要為ALOS 衛(wèi)星影像，因此本實(shí)驗(yàn)中選取的是PALSAR 影像各類(lèi)參數(shù)。設(shè)置隨機(jī)噪聲均方差為0.5rad。

完成各項(xiàng)模擬工作后，我們利用融合算法求解各PS 點(diǎn)上絕對(duì)下沉速率及高程改正參數(shù)值。首先，對(duì)所有PS及CR 點(diǎn)進(jìn)行基線網(wǎng)絡(luò)布設(shè)。采用Delaunay算法將這200 個(gè)PS 點(diǎn)與12 個(gè)CR點(diǎn)共同組成的212個(gè)一級(jí)控制點(diǎn)建立成參考基線網(wǎng)絡(luò)，共生成608條基線邊。再將600個(gè)擴(kuò)展PS點(diǎn)與212個(gè)一級(jí)控制點(diǎn)構(gòu)建成擴(kuò)展PS 網(wǎng)，在布網(wǎng)過(guò)程中限制二級(jí)網(wǎng)絡(luò)基線邊長(zhǎng)不超過(guò)5像素，則共生成2 295條基線邊。如圖5所示，圖中三角形表示CR 點(diǎn)位置。然后，依據(jù)式（1）構(gòu)建基線邊上任意兩個(gè)PS點(diǎn)的差分相位，并利用LAMBDA算法對(duì)式（1）進(jìn)行時(shí)間維相位解纏，求解出所有基線邊上的線性變形速率增量及高程改正值增量。最后，將CR 點(diǎn)上預(yù)先給出的線性變形速率及高程改正值模擬值作為基線網(wǎng)絡(luò)空間維解纏的起算數(shù)據(jù)，根據(jù)間接平差方法，求解出所有PS點(diǎn)上線性變形速率和高程改正結(jié)果。

圖5 模擬的參考基線網(wǎng)和擴(kuò)展基線網(wǎng)Fig.5 Simulated network of referenced PSs and extended PSs including CRs

4.2 實(shí)驗(yàn)結(jié)果

模擬實(shí)驗(yàn)預(yù)先給出了所有PS點(diǎn)的線性變形速率和高程改正模擬值，將其作為真值，與利用融合解算方法求解出的線性變形速率和高程改正計(jì)算值進(jìn)行比較，驗(yàn)證算法的精度。對(duì)參考PS 點(diǎn)和擴(kuò)展PS點(diǎn)上的線性變形速率計(jì)算值與模擬真值的差異作為其對(duì)應(yīng)精度，將精度分布直方圖表示在圖6中。從圖6明顯看出，無(wú)論參考PS點(diǎn)還是擴(kuò)展PS 點(diǎn)，線性變形速率計(jì)算值與模擬真值差異總體分布在－4～＋4mm 之間，說(shuō)明其與模擬真實(shí)結(jié)果十分接近。根據(jù)統(tǒng)計(jì)，無(wú)論參考網(wǎng)還是擴(kuò)展網(wǎng)，都有接近60%的PS點(diǎn)線性變形速率差異分布在－1～＋1mm，僅有個(gè)別PS點(diǎn)的差異超過(guò)10mm。計(jì)算得出，參考網(wǎng)點(diǎn)的線性速率差異均方根誤差為±0.37mm／a。

圖6 參考網(wǎng)PS點(diǎn)及擴(kuò)展網(wǎng)PS點(diǎn)線性變形速率計(jì)算值與模擬真值差異分布直方圖Fig.6 Statistics bars of Linear velocity rates errors at the reference PSs and extended PSs

圖7為真實(shí)PS點(diǎn)線性變形速率場(chǎng)與計(jì)算獲取的PS 點(diǎn)線性變形速率場(chǎng)的對(duì)比效果圖，從中可以看出二者顏色分布非常接近，但仍有部分區(qū)域存在細(xì)微的顏色差異。如圖中紅色方框標(biāo)出的A 區(qū)內(nèi)，計(jì)算結(jié)果中有約一半的點(diǎn)出現(xiàn)隆起，顏色表現(xiàn)為桔黃色，對(duì)應(yīng)的變形速率為0.02～0.04 cm／a，而模擬真值（圖7（b））中這一區(qū)域內(nèi)只有少量點(diǎn)出現(xiàn)隆起，大部分點(diǎn)表現(xiàn)為無(wú)形變。在B 區(qū)內(nèi)，計(jì)算結(jié)果（圖7（a））中出現(xiàn)了零星分布的下沉點(diǎn)，而模擬真值結(jié)果（圖7（b））中整個(gè)區(qū)域內(nèi)所有點(diǎn)都表現(xiàn)為隆起。出現(xiàn)這種情況主要是由于在模擬過(guò)程中加入了噪聲，另外，在算法求解過(guò)程中沒(méi)有考慮非線性形變及大氣延遲等因素的影響。在實(shí)驗(yàn)過(guò)程中獲取的高程改正值結(jié)果均方誤差為±0.5m，說(shuō)明利用這一算法進(jìn)行求解，可以通過(guò)增加高程改正的方式修正低精度的外部DEM，且修正結(jié)果精度為亞m級(jí)，可以減少在外部DEM 訂購(gòu)方面的成本，甚至可以不使用外部DEM 去除地形信息，直接獲取待求PS 點(diǎn)上的高程信息及線性變形速率。

圖7 計(jì)算獲取的PS點(diǎn)線性變形速率與模擬真值對(duì)比Fig.7 Comparison of linear velocities at all PSs with the simulated real value

5 結(jié) 語(yǔ)

本文將CRInSAR 與PSInSAR 兩種時(shí)間序列InSAR 地表變形監(jiān)測(cè)技術(shù)融合，用于解算地表線性變形速率，設(shè)計(jì)了詳細(xì)的算法流程。算法總體思想是利用CR 點(diǎn)上獲取的線性變形速率結(jié)果和高程改正值作為PS基線網(wǎng)絡(luò)的空間維解纏間接平差函數(shù)模型的起算數(shù)據(jù)，利用最小二乘原理求解出所有PS點(diǎn)上的線性變形速率與高程改正參數(shù)估值。這一算法可在一定程度上避免人為選取假設(shè)穩(wěn)定點(diǎn)的不確定性，融合了CR 與PS兩種技術(shù)的優(yōu)勢(shì)，可應(yīng)用于PS 點(diǎn)較為稀少且缺少先驗(yàn)變形信息的研究區(qū)域。為驗(yàn)證算法的可行性，本文設(shè)計(jì)并實(shí)現(xiàn)了模擬實(shí)驗(yàn)流程。從模擬實(shí)驗(yàn)的結(jié)果可以看出，這一融合算法解算線性變形速率可達(dá)到mm級(jí)精度，具備一定的可行性及可靠性。下一步將重點(diǎn)研究不同CR 約束網(wǎng)型對(duì)融合算法精度的影響，指導(dǎo)CR 點(diǎn)的安裝與布設(shè)，以起到對(duì)PS基線網(wǎng)絡(luò)的最佳約束效果。

［1］Ferretti A，Savio G，Barzaghi R，et al.Submillimeter Accuracy of InSAR Time Series：Experimental Validation［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（5）：1 142－1 153

［2］Crosetto M，Monserrat O，Iglesias R，et al.Persistent Scatterer Interferometry：Potential，Limits and Initial Cand X－Band Comparison［J］.Photogrammetric Engineering and Remote Sensing，2010，76（9）：1 061－1 069

［3］Osmano?lu B，Dixon T H，Wdowinski S，et al.Mexico City Subsidence Observed with Persistent Scatterer InSAR［J］.International Journal of Applied Earth Observation and Geoinformation，2011，13（1）：1－12

［4］Xia Y，Kaufmann H，Guo X.Landslide Monitoring in the Three Gorges Area Using D－InSAR and Corner Reflectors［J］.Photogrammetric Engineering and Remote Sensing，2004，70（4）：1 167－1 172

［5］Ferretti A，Prati C，Rocca F.Permanent Scatterers in SAR Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39（1）：8－20

［6］Ferretti A，Prati C，Rocca F.Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，2000，38（5）：2 202－2 212

［7］Kampes B M，Hanssen R F.Ambiguity Resolution for Permanent Scatterer Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing， 2004， 42（11）：2 446－2 453

［8］Xia Y，Kaufmann H，Guo X.Differential SAR Interferometry Using Corner Reflectors［C］.Geoscience and Remote Sensing Symposium，2002

［9］Liu G，Luo X，Chen Q，et al.Detecting Land Subsidence in Shanghai by PS－networking SAR Interferometry［J］.Sensors，2008，8（8）：4 725－4 741

［10］陳強(qiáng)，丁曉利，劉國(guó)祥，等.雷達(dá)干涉網(wǎng)絡(luò)的基線識(shí)別與解算方法［J］.地球物理學(xué)報(bào)，2009，52（9）：2 230－2 236（Chen Qiang，Ding Xiaoli，Liu Guoxiang，et al.Baseline Recognition and Parameter Estimation of Persistent－scatterer Network in Radar Interferometry［J］.Chinese Journal of Geophysics，2009，52（9）：2 230－2 236）

［11］陳強(qiáng)，丁小利，劉國(guó)祥.永久散射體雷達(dá)差分干涉應(yīng)用于區(qū)域地表沉降探測(cè)［J］.地球物理學(xué)報(bào)，2007，50（3）：737－743（Chen Qiang，Ding Xiaoli，Liu Guoxiang.Radar Differential Interferometry Based on Permanent Scatterers and Its Application to Detecting Regional Ground Subsidence［J］.Chinese Journal of Geophysics，2007，50（3）：737－743）

［12］Xia Y.CR－Based SAR－Interferometry for Landslide Monitoring［C］.Geoscience and Remote Sensing Symposium，2008