王艷， 張玲， 葛大慶， 張學(xué)東， 李曼

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

0 引言

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

依據(jù)地表形變?nèi)S矢量的構(gòu)成模型，升降軌InSAR在二維觀測(cè)量的條件下難以直接反演三維變化，其根本原因在于InSAR升降軌觀測(cè)量個(gè)數(shù)不夠。同時(shí)，受制于InSAR觀測(cè)對(duì)于形變的敏感程度和誤差影響，南北方向求解的穩(wěn)定性最差，其次為東西方向，垂直方向最為準(zhǔn)確[4-6]。在3個(gè)不同方向的InSAR觀測(cè)模式下具備了3組觀測(cè)數(shù)據(jù)，但解析模型仍然不存在多余觀測(cè)作為約束條件。盡管滿足解方程的條件，其結(jié)果的可靠性仍需要其余約束。因而，當(dāng)前對(duì)于三維形變量估算仍以觀測(cè)量的構(gòu)成模型，結(jié)合地表形變活動(dòng)模型進(jìn)行反演。根據(jù)雷達(dá)視線向?qū)嶋H觀測(cè)值對(duì)各個(gè)方向的敏感程度，在PSInSAR觀測(cè)條件下，對(duì)三維模型進(jìn)行簡(jiǎn)化可求解水平和垂直方向上的移動(dòng)量。

本研究中以升降軌觀測(cè)下的形變量模型構(gòu)成和相位觀測(cè)在各個(gè)方向上對(duì)形變量的敏感程度，將其簡(jiǎn)化處理以求解上下和東西2個(gè)方向的移動(dòng)量。試驗(yàn)中利用ENVISAT衛(wèi)星獲取的升降軌數(shù)據(jù)，對(duì)PSInSAR處理得到的形變速率進(jìn)行沉降和水平向移動(dòng)量反演試驗(yàn)，分別求解了水平移動(dòng)和垂直變形速率，結(jié)果表明試驗(yàn)區(qū)內(nèi)具有垂直移動(dòng)占主導(dǎo)，水平移動(dòng)極其微小的特點(diǎn)。

1 升降軌反演地表形變量的簡(jiǎn)化模型

融合升降軌InSAR觀測(cè)的形變量反演方法可以區(qū)分垂直方向和水平方向(東西方向)形變移動(dòng)量?？紤]到升降軌2組觀測(cè)是從2個(gè)幾乎相反的方向獲得數(shù)據(jù)(圖1)，對(duì)于地表同一目標(biāo)點(diǎn)而言，東西方向的水平移動(dòng)在2組觀測(cè)數(shù)據(jù)中的表現(xiàn)是基本相反的形變，而垂直方向上則表現(xiàn)為基本相同的形變。視線向形變量(dLOS)的構(gòu)成[4]為

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

(1)

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

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

dLOS≈UEsinθ+UVcosθ

，

(2)

式中右式第1項(xiàng)在降軌時(shí)取“+”號(hào)，升軌時(shí)取“-”號(hào)。由于入射角在感興趣區(qū)的變化極其微小，因此可用其均值來(lái)估計(jì)地表形變的水平和垂直分量。對(duì)應(yīng)于ENVISAT IS2數(shù)據(jù)，整景覆蓋的入射角變化低于6°，一般的監(jiān)測(cè)范圍內(nèi)均小于該值。由此，假定視線向產(chǎn)生1 cm的形變量，而由采用入射角均值所引起的最大誤差低于1 mm。

顯然，根據(jù)式(2)不能區(qū)分地表形變的南北方向的水平移動(dòng)。實(shí)際上，南北方向的水平移動(dòng)對(duì)垂向和東西向分量估計(jì)的影響非常微小，垂向分量估計(jì)幾乎不受南北向移動(dòng)的影響。假設(shè)目標(biāo)在南北方向產(chǎn)生1 cm的水平移動(dòng)，將其錯(cuò)誤估計(jì)為垂向移動(dòng)的誤差為0.8 mm，而在東西向的分量錯(cuò)誤估計(jì)為2 mm。

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

,

(3)

。

(4)

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

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

2 PSInSAR數(shù)據(jù)處理關(guān)鍵步驟

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

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

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

(5)

其中

(6)

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

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

3.1 升降軌下PSInSAR觀測(cè)量

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

圖2 研究區(qū)及升降軌數(shù)據(jù)分布

PSInSAR處理得到的形變速率圖為離散點(diǎn)模式，而難以直接提取同名點(diǎn)，因而需對(duì)2幅圖進(jìn)行插值處理，生成連續(xù)分布速率圖。圖3所示為升降軌模式下Track-39和Track-275的PSInSAR處理結(jié)果。在插值處理之前，先進(jìn)行坐標(biāo)系的統(tǒng)一，以消除

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

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

3.2 地表變形量求解

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

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

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

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

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

等大變形中表現(xiàn)較為普遍; 在地形起伏較顯著的地區(qū)(圖6B區(qū))，其水平變形量相對(duì)較明顯，這主要與地形有關(guān)，地形影響或主導(dǎo)了該地區(qū)發(fā)生形變的活動(dòng)方向，這類變形在滑坡等單一方向?yàn)橹鲗?dǎo)的移動(dòng)中較為普遍，而在地勢(shì)平坦地區(qū)不顯著。

4 結(jié)論

試驗(yàn)研究表明，利用簡(jiǎn)化的二維模型可實(shí)現(xiàn)對(duì)升降軌模式下InSAR觀測(cè)量的分解，從而可利用升降軌數(shù)據(jù)求解地表垂直方向和東西方向的變形。由于模型過(guò)于簡(jiǎn)化，適用性會(huì)受到很大的局限。對(duì)于文中研究區(qū)地勢(shì)起伏變化小的平原地區(qū)，升降軌模式對(duì)于沉降監(jiān)測(cè)影響均不明顯，其監(jiān)測(cè)結(jié)果都能較好地反映地表形變的實(shí)際情況，受雷達(dá)觀測(cè)方向影響不明顯。這為大范圍地面沉降監(jiān)測(cè)提供了依據(jù)，即利用統(tǒng)一模式觀測(cè)可滿足沉降監(jiān)測(cè)需要。

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

參考文獻(xiàn)(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] 夏耶.巴姆地震地表形變的差分雷達(dá)干涉測(cè)量[J].地震學(xué)報(bào),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雷達(dá)干涉揭示的巴姆地震(Mw6.5)3D同震形變場(chǎng)[J].遙感學(xué)報(bào),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] 葛大慶,王艷,張玲,等.德州地面沉降-回彈及地下水位波動(dòng)的InSAR長(zhǎng)時(shí)序監(jiān)測(cè)——以德州市為例[J].國(guó)土資源遙感,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地面沉降監(jiān)測(cè)互檢驗(yàn)與時(shí)序融合[J].國(guó)土資源遙感,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.