邢學敏 賀躍光 聞德保 周訪濱

1 長沙理工大學特殊環(huán)境道路工程湖南省重點實驗室，長沙市萬家麗南路960號，410014

2 長沙理工大學交通運輸工程學院，長沙市萬家麗南路960號，410014

3 中南大學地球科學與信息物理學院，長沙市麓山南路932號，410083

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

1 PSInSAR 算法

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

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

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

2 CRInSAR 算法

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

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

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

3 PSInSAR 與CRInSAR 融合解算算法

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

圖1 PS與CR 點相對位置示意圖Fig.1 Diagrammatic sketch of reference position between CRs and PSs

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

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

4 模擬實驗

4.1 實驗流程

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

圖3 模擬線性變形速率場Fig.3 The simulated linear velocity field

圖4 參考PS點與擴展PS點位置示意圖Fig.4 Distribution of reference PSs and extended PSs

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

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

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

4.2 實驗結果

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

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

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

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

5 結 語

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

［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］陳強，丁曉利，劉國祥，等.雷達干涉網(wǎng)絡的基線識別與解算方法［J］.地球物理學報，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］陳強，丁小利，劉國祥.永久散射體雷達差分干涉應用于區(qū)域地表沉降探測［J］.地球物理學報，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