王艷， 葛大慶， 張玲， 李曼， 郭小方， 王毅

(中國國土資源航空物探遙感中心,北京 100083)

0 引言

在觀測時間同步的條件下，同一地區(qū)不同觀測模式下獲取的PSInSAR監(jiān)測結(jié)果是對同一目標(biāo)的多角度觀測，這為升降軌2組觀測值的精度互檢驗提供了基礎(chǔ)。相對于傳統(tǒng)的地面測量檢驗而言，這種檢驗為InSAR觀測的內(nèi)部檢驗。若地表形變方向已知，則可以根據(jù)附加先驗知識的形變模型進(jìn)行分解，求解二維及三維變形量。該思路對地形變化較大地區(qū)的變形觀測更具實際意義[1-3]。在觀測對象移動矢量基本方向已知的條件下，利用單一軌道下的雷達(dá)視線向(light of sight,LOS)變形量也可實現(xiàn)主要移動方向移動量的估計。

由于視線向形變量對垂向變化更敏感，因而用InSAR測量以沉降為主的垂向變形的準(zhǔn)確性最高。對于升降軌獲取的2組視線向形變量，可將其轉(zhuǎn)換為沉降量后進(jìn)行沉降監(jiān)測的精度互檢驗。比較PSInSAR觀測形變量的基本條件是提取相干目標(biāo)的觀測值，而實際上，受雷達(dá)波入射方向和入射角的影響，升降軌圖像中相干目標(biāo)的分布密度和位置不盡相同，絕對準(zhǔn)確地提取每個相干目標(biāo)的觀測值是難以實現(xiàn)的，需采用統(tǒng)計方法分別提取同一組相干目標(biāo)對應(yīng)的2組PSInSAR觀測值(分別為升降軌)，進(jìn)行整體比較和點位比較。整體比較以PSInSAR沉降速率為主，用于估計2組觀測值的監(jiān)測精度； 點位比較以單個相干目標(biāo)的累積沉降量為主，比較2組觀測結(jié)果時間特征的差異。

相對于PSInSAR監(jiān)測獲取的形變速率觀測值，形變序列為另一觀測結(jié)果，其理論精度低于形變速率[4-5]。如果對升降軌模式形變序列數(shù)據(jù)進(jìn)行融合，則加密了單一觀測條件下的形變序列樣本集，能有效地揭示監(jiān)測對象時域上的非線性動態(tài)變化特征。

基于上述思路，本文提出了利用升降軌模式下2組PSInSAR觀測值進(jìn)行地面沉降監(jiān)測精度互檢驗及時序融合加密的方法。該方法以PSInSAR觀測的沉降速率為比較對象，在確定主輔軌道的基礎(chǔ)上對輔軌道觀測數(shù)據(jù)進(jìn)行基準(zhǔn)補償，以2組觀測值互差的均方差為檢驗指標(biāo)進(jìn)行統(tǒng)計檢驗和精度評價。在此基礎(chǔ)上，對輔軌道的觀測時序進(jìn)行基準(zhǔn)偏差補償，進(jìn)而實現(xiàn)2組觀測序列數(shù)據(jù)的時序融合。

1 觀測值坐標(biāo)系的統(tǒng)一

對比升降軌PSInSAR的監(jiān)測結(jié)果，首先需要統(tǒng)一坐標(biāo)系，使不同軌道獲取的形變觀測值(形變速率和形變序列)處于相同的參考坐標(biāo)系中，具有相同的空間基準(zhǔn)[5-6]。升降軌成像模式如圖1所示。

圖1升降軌觀測模式(左)及相干目標(biāo)分布(右)示意圖

Fig.1Ascendinganddescendingobservation(left)anddistributionofcoherenttarget(right)

坐標(biāo)系統(tǒng)的統(tǒng)一可通過2種方式來實現(xiàn)，即雷達(dá)坐標(biāo)系下的統(tǒng)一和地面坐標(biāo)系下的統(tǒng)一。前者是指將具有不同重疊度的SAR圖像直接進(jìn)行精確配準(zhǔn)； 后者則是先分別對圖像進(jìn)行地理編碼，實現(xiàn)正射校正后再進(jìn)行圖像匹配。受雷達(dá)成像方式的影響，在地形起伏較大地區(qū)的直接配準(zhǔn)誤差也較大。對于升降軌模式，由于雷達(dá)波入射方向和入射角的差異而產(chǎn)生的SAR數(shù)據(jù)不均勻變形會使這一方式的配準(zhǔn)精度受到限制。而對于地面坐標(biāo)系下的統(tǒng)一，要先進(jìn)行正射校正，即地理編碼，使得不同軌道下的雷達(dá)圖像位于相同的地面坐標(biāo)系下，以消除因地形起伏引起的畸變影響，進(jìn)而實現(xiàn)升降軌圖像的配準(zhǔn)。這種處理方式需要同時獲取研究區(qū)的DEM數(shù)據(jù)，而DEM的精度決定了地理編碼的精度，也決定了圖像的配準(zhǔn)精度。通常高程精度優(yōu)于10 m的DEM數(shù)據(jù)能保證平坦地區(qū)的配準(zhǔn)精度優(yōu)于1個像元。

本文采用地面坐標(biāo)系統(tǒng)一的方式進(jìn)行升降軌模式下PSInSAR測量結(jié)果的坐標(biāo)轉(zhuǎn)換。將升降軌下的2景SAR影像進(jìn)行地理編碼，并實現(xiàn)升降軌下PSInSAR處理獲取的2景沉降速率圖的坐標(biāo)轉(zhuǎn)換。確定升降軌中的主軌道，在獲取位于地面坐標(biāo)系下的雷達(dá)強度圖像后，為進(jìn)一步提高配準(zhǔn)精度，應(yīng)用多項式糾正完成主輔軌道的精確配準(zhǔn)。轉(zhuǎn)換函數(shù)為

(1)

式中： △x和△y分別為主輔影像x和y方向的相對偏移量；p為多項式階數(shù)，選擇3階多項式；a和b為多項式系數(shù)。

2 觀測值參考基準(zhǔn)的統(tǒng)一

PSInSAR觀測結(jié)果在平面上是“浮動的”，其絕對變化量取決于參考基準(zhǔn)的變化[6]。多軌道下參考基準(zhǔn)統(tǒng)一的本質(zhì)是解決不同參考位置間的整體偏差。受制于雷達(dá)波入射角和入射方向的影響，同一相干目標(biāo)在升降軌SAR觀測下對應(yīng)的位置略有差別。由于PSInSAR是對同一地區(qū)相同時段的觀測，因而測量值是對相同變形信號的描述。依據(jù)單個目標(biāo)進(jìn)行直接比較難以實現(xiàn)，需利用統(tǒng)計方法計算獲取。進(jìn)行參考位置的統(tǒng)一需選擇主軌道，其他軌道的觀測值均參考該條軌道。根據(jù)主軌道上相干目標(biāo)的空間位置，提取輔軌道上對應(yīng)的形變參數(shù)，進(jìn)行整體偏差求解。具體計算式為

(2)

3 觀測值的比較與互檢驗

本研究選擇覆蓋蘇州市西北區(qū)的ENVISAT數(shù)據(jù)。降軌為Track-275，升軌為Track-39。分別對2006年1月—2010年10月間2條軌道下的27景(降軌)和24景(升軌)SAR數(shù)據(jù)(均為IS2觀測角度)進(jìn)行時序分析處理，獲取了地面沉降速率和形變序列。對升降軌下的SAR圖像經(jīng)過坐標(biāo)統(tǒng)一后，得到地面坐標(biāo)系下沉降速率圖(圖2)?？梢钥闯?，研究區(qū)內(nèi)零星分布著多個沉降中心，最大沉降速率為39.1 mm/a，除沉降中心外，其他地區(qū)的沉降速率普遍小于10 mm/a。

圖2 降軌(左，Track-275)和升軌(右，Track-39)下PSInSAR 監(jiān)測的地面沉降速率

為了對升降軌觀測值進(jìn)行比較與互檢驗，以Track-275(降軌)為主軌道，利用該軌道覆蓋范圍內(nèi)相干目標(biāo)提取對應(yīng)的升軌(Track-39)下的PSInSAR觀測值(沉降速率)，共計提取53 361個相干點用于統(tǒng)計比較。按照式(2)求解Track-39 相對于Track-275 的整體均值偏差，并進(jìn)行Track-39 的基準(zhǔn)偏差修正。在數(shù)據(jù)分析過程中，首先對Track-39數(shù)據(jù)進(jìn)行插值處理，生成連續(xù)分布的沉降面?？紤]到相干目標(biāo)沉降值的影響范圍，利用鄰近點插值法，以100 m為半徑進(jìn)行升軌PSInSAR沉降速率結(jié)果的插值。對升降軌共有樣本點提取的2組數(shù)據(jù)進(jìn)行直方圖統(tǒng)計，得到如圖3所示的結(jié)果。

圖3 降軌(左，Track-275)和升軌(右，Track-39)模式下地面沉降速率統(tǒng)計直方圖

圖4為升降軌PSInSAR觀測沉降速率在基準(zhǔn)補償前后的相關(guān)統(tǒng)計散點分布圖。比較二者可知，2組觀測結(jié)果的相關(guān)性在補償前后一致，主要差異為2組觀測值的整體偏差。

圖4降軌(Track-275)和升軌(Track-39)基準(zhǔn)補償前(左)后(右)地面沉降速率相關(guān)統(tǒng)計

Fig.4Correlationofsubsidencevelocitybefore(left)andafter(right)referenceeffsetcompensationofdescending(Track-275)andascending(Track-39)

由圖4可知，升降軌模式下研究區(qū)地面沉降速率的整體分布符合線性函數(shù)特征。因而，可將基準(zhǔn)偏差補償后的沉降速率(Track-39)與主軌道(Track-275)的沉降速率按照

Y=A+BX

(3)

進(jìn)行線性統(tǒng)計。式中：X為升軌(輔軌道，Track-39)基準(zhǔn)補償后的PSInSAR監(jiān)測速率(圖5)；Y為降軌(主軌道，Track-275)下PSInSAR觀測速率；A為整體偏差；B為線性函數(shù)的斜率。

圖5 Track-39升軌基準(zhǔn)補償后PSInSAR監(jiān)測的地面沉降速率圖

經(jīng)過函數(shù)擬合得到A=0.17 mm，B=1.008。B值接近于1，表明二者在修正后沉降速率的參考基準(zhǔn)不存在顯著差異，線性函數(shù)特征明顯。

為檢驗2組數(shù)據(jù)的相關(guān)性，直接利用升降軌PSInSAR觀測值進(jìn)行統(tǒng)計比較，以互差的均方差m’為統(tǒng)計指標(biāo)，計算式為，

(4)

式中：Xi為降軌觀測值；Yi為升軌觀測值(經(jīng)過基準(zhǔn)偏差補償)；P為樣本點個數(shù)。對提取的53 361個相干點進(jìn)行統(tǒng)計，得m’=1.824 mm。結(jié)果表明了研究區(qū)地表以垂向變化為主，水平移動量極為微小，證明了文獻(xiàn)[8]試驗結(jié)果的正確性。同時也表明，PSInSAR方法在監(jiān)測以垂向下沉為主的地面沉降時的可靠性。當(dāng)雷達(dá)數(shù)據(jù)為20～30景時，沉降速率監(jiān)測互檢驗的精度優(yōu)于2 mm。

4 觀測序列的時序融合

升降軌PSInSAR數(shù)據(jù)處理方法如文獻(xiàn)[7]所述。 假設(shè)升軌模式下有M景SAR影像，對應(yīng)的形變觀測序列值為

(5)

降軌模式下有N景SAR影像，對應(yīng)的形變觀測序列值為。

(6)

顯然，2組觀測序列存在參考基準(zhǔn)整體速率偏差和因起始時間差所引起的形變累積量差，需要逐個補償。對輔軌道下(這里以升軌為例)的形變序列進(jìn)行速率整體偏差修正，求得參考基準(zhǔn)補償后相干目標(biāo)的形變(沉降)序列為

(7)

在完成速率偏差修正后，需將2組觀測值統(tǒng)一到相同的時間起點。此時，形變序列的補償量為輔軌道相對于主軌道觀測序列起始時刻的形變量偏差，將之加到序列累積量上即可，即,

(8)

在完成2種參考基準(zhǔn)的補償后，主輔軌道下的形變序列便統(tǒng)一到相同的時間和空間參考基準(zhǔn)上，得到相干目標(biāo)加密后的形變(沉降)序列為

D=d0,d1,…，dM-1,…,dM+N-1。

(9)

經(jīng)上述過程得到每個相干目標(biāo)獨立的觀測形變序列。圖6分別為樣本點P1和P2(見圖5)的沉降

圖6 升降軌下PS點時序融合后的沉降序列

序列。其中，P1點沉降較為緩慢，小于10 mm/a，加密后的結(jié)果顯示2組觀測結(jié)果整體速率基本一致，非線性特征不明顯；P2點在局部時段內(nèi)有顯著的非線性變化特征，其整體沉降速率超過20 mm/a，觀測時段內(nèi)的累積沉降量達(dá)100 mm以上?？傮w而言，升降軌下的2組觀測值的時間變化特征一致，經(jīng)參考基準(zhǔn)整體補償后的時間序列降低了速率偏差的時間積累效應(yīng)，而融合后的加密觀測則凸顯了非線性變化特征，表明多角度InSAR聯(lián)合觀測對于地面沉降監(jiān)測在時序分析上的積極意義。

5 結(jié)論

1)升降軌PSInSAR地面沉降觀測值的比較實現(xiàn)了PSInSAR觀測結(jié)果的內(nèi)部檢驗。在SAR數(shù)據(jù)量為20～30景的條件下，地面沉降速率監(jiān)測結(jié)果的互檢驗精度優(yōu)于2 mm，觀測樣本的差異主要為參考基準(zhǔn)的整體偏差。

2)對時序觀測序列的時間和空間偏差修正實現(xiàn)了升降軌形變序列的時序融合，加密了單一軌道下的觀測序列，可有效揭示非線性變化過程，為時序結(jié)果的檢驗提供了依據(jù)。

3)從廣義InSAR形變監(jiān)測技術(shù)上，本研究完成了形變監(jiān)測結(jié)果精度的互檢驗與2組獨立形變觀測值的時序融合2項內(nèi)容。前者為檢驗InSAR測量的精準(zhǔn)度提供了方法和證據(jù)； 后者有助于InSAR地表形變監(jiān)測結(jié)果精度的提高，并能更好地提取形變的非線性特性。

參考文獻(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,Parson B,Wright T J,et al.Surface displacements and source parameters of the 2003 Bam(Iran)earthquake from Envisat advanced synthetic aperture Radar imagery[J].Journal of Geophysical Research,2005,110(B9):B09406.

[3] Fielding E J,Talebian M,Rosen P A,et al.Surface ruptures and building damage of the 2003 Bam,Iran,earthquake mapped by satellite synthetic aperture Radar interferometric correlation[J].Journal of Geophysical Research,2005,110(B3):B03302.

[4] Mora O,Mallorqui J J,Broquetas A,et al.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Trans Geosci Remote Sens,2003,41(10):2243-2253.

[5] Ketelaar G,Van L F,Marinkovic P,et al.Multi-track PS-InSAR datum connection[C]//Proceeding of IEEE International Geoscience and Remote Sensing Symposium.Barcelona:IEEE,2007:2481-2484.

[6] 葛大慶,王艷,張玲,等.低相干條件下區(qū)域性地面沉降InSAR調(diào)查與監(jiān)測技術(shù)研究[R].北京:中國國土資源航空物探遙感中心,2009.

Ge D Q,Wang Y,Zhang L,et al.Research on key InSAR technology to monitor regional subsidence occurred in low correlated regions[R].Beijing:China Aero Geophysical Surveying and Remote Sensing Center for Land and Resources,2009.

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

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

[8] 王艷,張玲,葛大慶,等.升降軌PSInSAR觀測反演沉降與水平向位移試驗[J].國土資源遙感,2014,26(4):97-102.

Wang Y,Zhang L,Ge D Q,et al.Experimental study on vertical and horizontal displacement retrieval by jointly analysis of ascending and descending PSInSAR data[J].Remote Sensing for Land and Resources,2014,26(4):97-102.