升降軌PSInSAR觀測反演沉降與水平向位移試驗

王艷， 張玲， 葛大慶， 張學東， 李曼

(1.中國國土資源航空物探遙感中心,北京 100083; 2.北京建筑大學測繪與城市空間信息學院，北京 100044)

0 引言

雷達干涉測量(InSAR)所獲取的形變監測結果為視線向形變量，可依據形變方向與視線向的投影關系，進行不同方向形變量的求解。上下、東西和南北3個方向上形變量求解的前提是具備多維(3個方向)觀測值，但實際上，受制于雷達的成像方式，難以同時滿足3個方向的獨立觀測。因而，直接利用InSAR監測結果反演三維形變量存在一定的難度。實際應用中，如地震等大尺度形變場，一般存在不同方向上的形變特征，可以利用升降軌、左右視觀測[1-3]，以及圖像偏移量跟蹤方法等實現不同方向形變量的求解。

依據地表形變三維矢量的構成模型，升降軌InSAR在二維觀測量的條件下難以直接反演三維變化，其根本原因在于InSAR升降軌觀測量個數不夠。同時，受制于InSAR觀測對于形變的敏感程度和誤差影響，南北方向求解的穩定性最差，其次為東西方向，垂直方向最為準確[4-6]。在3個不同方向的InSAR觀測模式下具備了3組觀測數據，但解析模型仍然不存在多余觀測作為約束條件。盡管滿足解方程的條件，其結果的可靠性仍需要其余約束。因而，當前對于三維形變量估算仍以觀測量的構成模型，結合地表形變活動模型進行反演。根據雷達視線向實際觀測值對各個方向的敏感程度，在PSInSAR觀測條件下，對三維模型進行簡化可求解水平和垂直方向上的移動量。

本研究中以升降軌觀測下的形變量模型構成和相位觀測在各個方向上對形變量的敏感程度，將其簡化處理以求解上下和東西2個方向的移動量。試驗中利用ENVISAT衛星獲取的升降軌數據，對PSInSAR處理得到的形變速率進行沉降和水平向移動量反演試驗，分別求解了水平移動和垂直變形速率，結果表明試驗區內具有垂直移動占主導，水平移動極其微小的特點。

1 升降軌反演地表形變量的簡化模型

融合升降軌InSAR觀測的形變量反演方法可以區分垂直方向和水平方向(東西方向)形變移動量。考慮到升降軌2組觀測是從2個幾乎相反的方向獲得數據(圖1)，對于地表同一目標點而言，東西方向的水平移動在2組觀測數據中的表現是基本相反的形變，而垂直方向上則表現為基本相同的形變。視線向形變量(dLOS)的構成[4]為

dLOS=(UNsinφ-UEcosφ)sinθ+UVcosθ。

(1)

圖1 升降軌觀測模式下的垂向與東西向形變分解示意圖

式中:UN，UE和UV分別表示北南、東西和垂直方向的地表移動分量；φ為方位角(正北與傳感器飛行方向夾角，自北順時針為正)，按照習慣降軌方向設為正，升軌方向為負；θ為入射角。對于ENVISAT雷達衛星掃描模式下的IS2數據而言，升降軌條件下其方位角分別為-11.9°和191.9°，由此可近似得到sinφ≈0(為-0.19)，cosφ≈±1(升軌時為0.98，降軌時為-0.98)； 代入式(1)中，視線向形變量表達式可進一步簡化為

dLOS≈UEsinθ+UVcosθ

，

(2)

式中右式第1項在降軌時取“+”號，升軌時取“-”號。由于入射角在感興趣區的變化極其微小，因此可用其均值來估計地表形變的水平和垂直分量。對應于ENVISAT IS2數據，整景覆蓋的入射角變化低于6°，一般的監測范圍內均小于該值。由此，假定視線向產生1 cm的形變量，而由采用入射角均值所引起的最大誤差低于1 mm。

顯然，根據式(2)不能區分地表形變的南北方向的水平移動。實際上，南北方向的水平移動對垂向和東西向分量估計的影響非常微小，垂向分量估計幾乎不受南北向移動的影響。假設目標在南北方向產生1 cm的水平移動，將其錯誤估計為垂向移動的誤差為0.8 mm，而在東西向的分量錯誤估計為2 mm。

因而，由式(2)可推出升降軌觀測模式下東西和垂直方向的形變量求解公式為

,

(3)

。

(4)

式中:dLOS(desc)表示降軌數據的視線向形變量；dLOS(asc)表示升軌數據的視線向形變量。

利用式(3)(4)分別解算東西和垂直方向的地表移動分量的前提是要求升降軌模式下同一觀測區域內具有相同的觀測對象，即分別在升軌和降軌解算數據中具有同名像元。這在相干性良好的差分干涉圖解算地震形變場移動分量中較為常見。而對于非線性形變場中相干目標的解算，該方法的缺點是，在地形起伏較大的地區，相對于單個數據集可用的像元數目，用來估計東西向和垂向的地表移動分量的像元數目較少。這是由每個像元在升軌和降軌數據集中固有的成像方式及其對相干性的影響所引起的。由于單個目標的解算依賴于對同名相干像元的識別，而實際上升降軌下完全同名點的數目有限，對于估計整體形變場不利。為此,針對升降軌PSInSAR觀測值，可將升降軌中的一景作為主影像，另一景為輔影像，以此來插值提取相同位置的變形量，從整體上反演地表形變的水平和垂直位移。

2 PSInSAR數據處理關鍵步驟

永久散射體干涉測量(PSInSAR)技術的核心是對相干像元，即永久散射體(permanent scatter)的差分干涉相位序列進行分析，根據差分干涉相位各分量的時空特征，估算大氣波動影響、DEM誤差以及噪聲等，將其從差分干涉相位中逐個分離，最終獲取每個PS的線性和非線性形變速率以及DEM誤差等參數。

干涉像對組合是PSInSAR處理的重要過程，其目的是將干涉像對進行優化組合。由于地形相位誤差和相干性均受空間基線的影響，而短基線條件有利于相干性的保持，因而選擇空間基線小于給定閾值，如ENVISAT可選300 m的像對生成干涉圖。同時，顧及到地表形變速率的大小，對于一般地面沉降監測而言，選擇時間間隔小于2 a，大于2～3個衛星重復周期的像對生成干涉紋圖，在多數情況下是適用的。

對識別出的相干目標候選點的差分相位序列的互差進行處理是PSInSAR處理的關鍵所在，其過程是以Delaunay三角網構建不規則格網連接將其進行連接，對鄰近點pm和pn的干涉相位時間序列進行二次求差，應用二維周期圖估計點間形變速率和高程誤差。式(5)為利用點間相位差估計形變速率和高程誤差的函數模型，在完成頻率估計后，利用最小二乘法實現整個三角網的積分，生成形變速率圖，以反映形變場分布特征。

(5)

其中

(6)

式中：γmodel為模型相關系數；pm和pn分別為2景影像同名點對應的第m和第n個點；N為干涉像對個數；φdiff表示差分相位；Ti為第i景干涉紋圖的時間基線；φmodel表示模型相位；λ為雷達波波長；γi為第i景干涉紋圖的相關系數；bi為第i景干涉紋圖的垂直基線；θ為雷達入射角；υ(pm)和υ(pn)分別表示第m點和第n點的形變速度；ε(pm)和ε(pn)分別為第m點和第n點的DEM誤差。

3 升降軌PSInSAR形變量反演與分析

3.1 升降軌下PSInSAR觀測量

研究需選擇同時獲取了升軌和降軌模式下ENVISAT-ASAR數據的地區(如圖2中紅色框和粉色框的重疊區)。本研究選定的圖像(圖2中的藍色框區域)覆蓋范圍為30 km×30 km，位于蘇州市區西北周邊，分別獲取了2006—2010年間升軌(Track-39)和降軌(Track-275)雷達數據24景和27景。

圖2 研究區及升降軌數據分布

PSInSAR處理得到的形變速率圖為離散點模式，而難以直接提取同名點，因而需對2幅圖進行插值處理，生成連續分布速率圖。圖3所示為升降軌模式下Track-39和Track-275的PSInSAR處理結果。在插值處理之前，先進行坐標系的統一，以消除

圖3 升(左)降(右)軌模式下PSInSAR視線向形變速率(空間插值后)

由于解算過程中參考位置不同而產生的基準偏差[7-8]。坐標系的統一采用地形校正方法，然后在相同的地面坐標系下同時完成2幅SAR影像的地理編碼。地理編碼的誤差均小于0.15個像元，2圖像之間不再進行二次匹配。參考基準偏差的補償利用同名點統計整體偏差值，進而補償形變場整體差異。最終得到的2景圖像中各形變中心分布位置相同，監測結果量值相近(圖3)。

3.2 地表變形量求解

根據式(2)(3)對入射角進行歸一化(以正弦為例)處理，按照同一模式下的入射角進行分析。圖4和圖5分別為研究區入射角歸一化結果和方位角歸一化結果。

圖4 升(左)降(右)軌模式下歸一化入射角

圖5 升(左)降(右)軌模式下歸一化方位角

在升降軌解算時，對于入射角的考慮為2組數據集里每個像元的入射角，而實際上這種情況下難以進行參數估計。因此，結合本研究區范圍較小的實際，對每個像元均按照統一入射角來計算。圖4中入射角的最大差別為0.05 rad，圖5中方位角的最大差別為0.002 rad，表明二者的變化極其微小，對入射角變化的簡化和對方位角變化的忽略處理是正確的。

利用升降軌視線向觀測量解算垂向與東西向形變速率(圖6)。解算結果表明，研究區垂向變形占主導，即以地面沉降為主，而水平方向變化極其微小。主要沉降變形區(圖中A區)水平變形較為顯著，這主要是在地勢平坦地區由于不均勻沉降引發地表曲率變化，進而產生向沉降中心方向的輕微滑動。對照圖6中左圖所示，水平變形相對顯著的地區與沉降中心位置相同。這種形變類型在礦山沉陷

等大變形中表現較為普遍; 在地形起伏較顯著的地區(圖6B區)，其水平變形量相對較明顯，這主要與地形有關，地形影響或主導了該地區發生形變的活動方向，這類變形在滑坡等單一方向為主導的移動中較為普遍，而在地勢平坦地區不顯著。

4 結論

試驗研究表明，利用簡化的二維模型可實現對升降軌模式下InSAR觀測量的分解，從而可利用升降軌數據求解地表垂直方向和東西方向的變形。由于模型過于簡化，適用性會受到很大的局限。對于文中研究區地勢起伏變化小的平原地區，升降軌模式對于沉降監測影響均不明顯，其監測結果都能較好地反映地表形變的實際情況，受雷達觀測方向影響不明顯。這為大范圍地面沉降監測提供了依據，即利用統一模式觀測可滿足沉降監測需要。

實際上，升降軌觀測僅提供了二維觀測量，仍然不能滿足多參數反演的需要。受制于地形和地表形變類型的影響，結合高分辨率雷達觀測數據，針對形變模型的參數反演將更為有效。這也是今后需要進一步研究的方向。

參考文獻(References)：

[1] Fialko Y,Sandwell D,Simons M,et al.Three-dimensional deformation caused by the Bam,Iran,earthquake and the origin of shallow slip deficit[J].Nature,2005,435(7040):295-299.

[2] Funning G J,Parsons B,Wright T J,et al.Surface displacements and source parameters of the 2003 Bam(Tran)earthquake from Envisat advanced synthetic aperture radar imagery[J].Journal of Geophysical Research,2005,110(B9):B09406.

[3] Bonforte A F,Guglielmino M,Coltelli A,et al.Structural assessment of mount etna volcano from permanent scatterers analysis[J].Geochem Geophys Geosyst,12(2),doi:10.1029/2010GC003213.

[4] Hanssen R F.Radar Interferometry-data Interpretation and Error Analysis[M].New York:Kluwer Academic Publishers,2002:12-65.

[5] 夏耶.巴姆地震地表形變的差分雷達干涉測量[J].地震學報,2005,27(4):423-430.

Xia Y.Bam earthquake:Surface deformation measurement using Radar interferometry[J].Acta Seismologica Sinica,2005,27(4):423-430.

[6] 孫建寶,梁芳,徐錫偉,等.升降軌道ASAR雷達干涉揭示的巴姆地震(Mw6.5)3D同震形變場[J].遙感學報,2006,10(4):489-496.

Sun J B,Liang F,Xu X W,et al.3D co-seismic deformation field of the Bam earthquake(Mw6.5) from ascending and descending pass ASAR Radar interferometry[J].Journal of Remote Sensing,2006,10(4):489-496.

[7] 葛大慶,王艷,張玲,等.德州地面沉降-回彈及地下水位波動的InSAR長時序監測——以德州市為例[J].國土資源遙感,2014,26(1):103-109.

Ge D Q,Wang Y,Zhang L,et al.Seasonal subsidence-rebound and ground water level changes monitoring by using coherent target InSAR technique: A case study in Dezhou,Shandong[J].Remote Sensing for Land and Resources.2014,26(1):103-109.

[8] 王艷,葛大慶,張玲,等.升降軌PSInSAR地面沉降監測互檢驗與時序融合[J].國土資源遙感,2014,26(4):125-130.

Wang Y,Ge D Q,Zhang L,et al.Inter-comparison and time series fusion of ascending and descending PSInSAR data for land subsidence monitoring[J].Remote Sensing for Land and Resources,2014,26(4):125-130.