利用永久散射體雷達干涉技術進行太原市地面沉降監測

王 琪

(1. 山西省煤炭地質物探測繪院，山西 晉中 030600; 2. 太原市基礎地理數據中心，山西 太原 030009)

一、引 言

地面沉降是指由于自然因素或人為活動引發地殼表層松散土層壓縮并導致地面標高降低的地質現象。近年來，差分干涉測量(DInSAR)技術已成為地面沉降監測的重要手段，具有全天候、大范圍監測的優勢，理論上可以達到厘米級的精度。特別是近幾年發展的時間序列InSAR技術(time series InSAR)，如永久散射體雷達干涉測量PS-InSAR方法和短基線集雷達干涉測量SBAS InSAR方法，是在傳統DInSAR基礎上發展起來的更具潛力的新的監測方法，其克服了常規DInSAR中的時空失相干，并減弱了大氣效應等因素的影響，大大提高了監測的精度，在城市等緩慢形變區域具有極大的應用潛力和優勢。

本文將覆蓋太原市2003-08—2010-09年間39景Envisat ASAR影像作為試驗數據，采用Andy Hooper開發的StaMPS軟件的PS方法對太原市的地面沉降進行試驗研究，提取該時期的沉降速率圖，研究該區域地表形變特征，并對地面沉降成因進行探討分析。

二、StaMPS永久散射體雷達干涉技術(PS-InSAR)

A Ferretti，C Prati和F Rocca等于1999年共同提出并驗證了永久散射體(PS-InSAR)方法，即在一組長時間序列SAR影像中，選取保持高相干性且較為穩定的點為PS點目標，其受時間空間去相干的影響小[1]。通過分析這些點目標的相位組成，利用各相位分量的時頻特性，可以分離出大氣相位，從而獲得較為精確的地表形變信息。這種算法比較適合有大量人工建筑的城區，而在非城鎮區域或地表形變速率不規則的地區選取的可靠PS點較少。Andy Hooper于2004年提出了一種新的PS算法，該算法采用幅度離散指數和相位空間相關性相結合的方法來識別PS點，不需要先驗知識，可以在任何地形區域選取PS點，極大擴展了PS技術的應用范圍。本文正是基于這個方法對太原市的地面沉降進行了試驗。

假設N+1幅SAR影像，選擇一幅各個方面(時間基線、空間基線、多普勒質心頻率差等)最優的作為主影像，生成N幅干涉條紋圖，引入外部參考DEM去除地形相位的影響,生成差分干涉圖。首先基于Ferretti提出的幅度離散指數DA(amplitude dispersion index)提取PS候選點

式中，σA，μA為N幅干涉條紋圖振幅的標準差和均值。設置DA閾值(0.4～0.42)初步選取PS候選點，對選出的PS候選點集利用相位的空間相關性分析其相位穩定性，據此篩選最佳PS點，基本模型如下：第i幅干涉圖中的第x個PS候選點的差分干涉相位可以表示為

ψx,i=W{φD,x,i+φA,x,i+ΔφS,x,i+Δφθ,x,i+φN,x,i}

(2)

(3)

為此，可以設定用以篩選PS點的相關系數γx為

(5)

式中，N為干涉圖的數量；γx即為判定PS點目標的指標。但僅僅如此還不夠，StaMPS結合信噪比、γx和PS概率進行循環迭代來選取最終的PS點。

(6)

三、太原市地面沉降監測試驗

1． 試驗區概況和試驗數據

太原位于山西省太原斷陷盆地的北緣，屬于汾渭盆地，是我國中西部地區地面沉降最為嚴重的城市之一。太原市從20世紀50年代末發現地面沉降，80年代至90年代為地面沉降快速發展階段，至2004年形成了西張和城區兩個沉降區，西張、萬柏林、下元、吳家堡4個沉降漏斗中心。沉降區北起上蘭鎮，南至劉家堡鄉郝村，西抵西鎮，東達榆次西河堡村；南北長約39 km，東西寬約15 km，沉降面積達548 km2。90年代以來，沉降范圍逐年向盆地邊緣擴展，沉降漏斗面積逐年擴大，南部有向晉中盆地延伸趨勢，已影響到該地區社會經濟的可持續發展[3]。

試驗選取時間跨度從2003年8月17日至2010年9月19日的39景歐空局30 m分辨率的Envisat ASAR單視復數影像(SLC)。ASAR的中心工作頻率為C波段(頻率為5.331 GHz)，成像中心入射角為23°，影像的覆蓋范圍如圖1所示。

圖1 ASAR數據覆蓋范圍

2. 差分干涉處理

首先利用Doris軟件對39景影像經過配準、干涉、去平等步驟進行差分干涉處理，在去除地形相位時，引入了3弧秒分辨率的SRTM DEM數據，平面精度為90 m，高程精度約為10 m。影像間的空間基線越長，對應的地形相位誤差越大。以2009年7月26日為主影像，39景影像的空間基線分布最大856 m，最小40 m，平均為197 m，可以看出基線分布比較均勻，且太原城區地形平坦，從而由地形引起的誤差相位將非常小。經過差分干涉獲得了38對干涉圖，如圖2所示。

圖2 差分干涉圖

3. PS干涉處理

StaMPS的PS處理分為以下幾個步驟：①PS候選點初步選取；②估計PS候選點相位噪聲；③根據相位噪聲特性選取PS點；④進一步篩選出最終PS點；⑤對PS點的纏繞相位進行DEM誤差校正；⑥三維相位解纏；⑦估計空間相關相位誤差；⑧估計其他空間相關誤差項[4]。

步驟①中，設置DA閾值為0.42，取較大的值可以選出更多的PS候選點，但同時存在較多的偽PS點，通過后續的步驟將不斷進行篩選，最終將偽PS點去除。步驟②中，γx是通過估計空間相關相位項和空間不相關地形誤差來定義的(式(5))，將像素采樣到50 m格網空間從而使其滿足空間相關條件，然后通過迭代方法估計空間相關相位項，設置迭代收斂條件為γx的變化值小于0.005。步驟③選取PS點的γx閾值是離差指數DA及偽PS點密度的函數，通過概率統計的方法設定，在處理中設置可接受的偽PS點密度為18/km2。步驟④首先剔除那些由于地理編碼誤差導致的偽PS點，以及受鄰近PS點的影響而被認為是PS點的偽PS點，然后估計各PS點相位噪聲，刪除噪聲比較低的點，經過這一步篩選之后剩余36 895個點，即為最終的PS點。第⑤步對上述的PS點的纏繞相位進行空間非相關地形誤差校正。第⑥步利用對上述校正后的相位進行三維解纏，解纏的參考為以地理坐標(112.530 2°，37.959 5°)為中心，200 m范圍內的29個PS點。可以看出，第1、3、4、9景影像空間垂直基線比較大，致使其基線失相關嚴重，從解纏的相位圖中也可以看出這4景影像存在較大解纏誤差，因此在解纏處理時通過“drop_ifg_index”參數將這4景影像排除，重新解纏并估計空間相關相位項。這樣獲得的解纏相位圖如圖3所示。利用解纏的相位觀測值，反演了地表雷達視線向平均形變速率(MLV)，如圖4所示；同時得到了形變速率的標準差，如圖5所示。

圖3 解纏相位圖

圖4 地表形變速率(mm/a) 圖5 形變速率標準差(mm/a)

4. 地面沉降空間分布特征

根據太原市沉降監測資料及沉降歷史特征(如圖6所示)，PS-InSAR方法獲得的地表形變(如圖5所示)在空間分布上與以往的研究成果基本一致[5-6]。將PS-InSAR監測結果投影到Google影像上并標注了太原市沉降漏斗中心的分布(如圖6所示)。由圖6可以看出，在2003-08—2010-09時間段內太原市地面沉降主要分布在中部和南部地區。北部平原地區沉降速率較小，多數不超過-10 mm/a，部分地區存在反彈回升現象。東西兩側太行山、呂梁山山緣洪積扇沖積平原的沉降速率也較小，大多不超過-25 mm/a。總體來看，太原市地面沉降空間分布從下元沿著汾河往南呈喇叭狀展開，沉降速率由北向南逐漸增大，研究區范圍內最大沉降速率位于劉家堡鄉附近，達到-45 mm/a；其次為孫家寨，達到-40 mm/a。地面沉降漏斗分布由北到南依次為：西張，大部分地區已經出現反彈現象，平均反彈速率為2 mm/a；萬柏林，平均沉降速率為-21.5 mm/a；下元，平均沉降速率-32.4 mm/a；吳家堡，平均沉降速率為-25.9 mm/a；小店鎮，平均沉降速率為-30.9 mm/a;孫家寨，平均沉降速率為-34.8 mm/a；劉家堡，平均沉降速率為-36.9 mm/a。從圖5形變速率的標準差可以看出，最大的形變速率標準差發生在小店沉降中心，最大為2 mm/a，這可能是該地區存在大量耕地，由耕地的季節性翻耕導致地面散射特性變化引起的；其他地區則更小，這說明估計的形變結果是可靠的。

圖6 地面沉降速率在Google地圖上的分布(mm/a)

與2000年之前相比，該時間段內太原的地面沉降無論是在空間分布還是在沉降快慢程度上都發生了變化。1956—2000年以前太原的地面沉降如圖7所示，含西張、萬柏林、下元、吳家堡4個沉降中心，沉降中心沿汾河呈串狀分布，最大沉降中心位于吳家堡，南部的小店鎮和孫家寨沉降漏斗則沒有形成。在研究時間范圍2003-08—2010-09內，太原地面沉降分布向南移動，出現了小店和孫家寨等新的沉降中心，發展成了一條小店至孫家寨的帶狀沉降區域。過去的西張沉降中心趨于穩定，部分地區出現了反彈回升現象；萬柏林沉降速率大為減少，沉降面積也縮小很多；下元和吳家堡仍然為較明顯的沉降中心，且兩者有連成一片的趨勢，但沉降速率與1989—2000年相比已經明顯減緩。表1對比了太原市各沉降中心在1956—2000年與2003—2010年各時期內沉降速率的變化情況。總之，該時期地面累積沉降量由北往南逐漸減小，沉降速率逐漸增大。這說明隨著2003年太原市采取關井壓采及引黃入晉等一系列措施以來，原有地面沉降中心得到了一定程度的控制。

圖7 太原不同歷史時段累積地面沉降演化(太原市地面沉降監測報告，2004)

表1 不同時期太原地面沉降中心平均沉降速率對比 (mm/a)

注：1956—2000年沉降數據來自文獻[7]。

5. 地面沉降時間演變特征

同時，在沉降中心(下元、萬柏林、吳家堡、小店、孫家寨)提取PS點的地面沉降時間序列，通過簡單的線性回歸模型擬合出形變值隨時間的變化，如圖8—圖11所示。小圓圈代表PS點，實線表示擬合的線性形變時間序列，兩條虛線表示兩個連續的PS點之間的相位差不能超過π，這是由于相位存在2π模糊度的原因。在這里所用數據為Envisat ASAR數據，波長5.6 cm，因此虛線代表±14 mm。從時間序列圖可以看到，2007-04-11、2007-12-09、2008-01-13、2008-02-17這4景干涉圖上PS點之間都超出了虛線，發現這4景干涉圖都是在冬天季節，很有可能由于下雪覆蓋的原因使得相干性下降，從而探測到的信號變得不可靠。

圖8 下元沉降時間序列

圖9 萬柏林沉降時間序列

圖10 吳家堡沉降時間序列

圖11 小店鎮沉降時間序列

四、結論與討論

本文利用2003-08—2010-09的39景Envisat ASAR 數據基于StaMPS的PS-InSAR技術計算了該時期太原市的地面沉降，獲得了沉降速率分布圖。結果表明，StaMPS不僅在城區能夠探測到較多可靠的PS點，并且在有植被的地區也能提取到PS點目標，但是在這些地區提取的形變結果并不是很可靠。通過與原有的沉降資料對比驗證，證實該方法在城市地面沉降監測中的可靠性。但是，該方法僅能夠探測到雷達視線向(LOS)的線性形變，未能對非線性形變進行估計，這對于StaMPS的推廣應用將是不利的。由于缺少必要的數據，沒有獲得如(U,E,N)等方向的形變值。未來可聯合GPS連續運行參考站(CORS)的數據進行解算，從而全面掌握更為精細的地表三維形變場。

參考文獻：

[1] 廖明生，林暉.雷達干涉測量—原理與信號處理基礎[M].北京：測繪出版社，2003.

[2] HOOPER A,SEGALL P,ZEBKER H. Persistent Scatterer Interferometric Synthetic Aperture Radar for Crustal Deformation Analysis, with Application to Volcán Alcedo,Galápagos[J].Journal of Geophysical Research,2007,112(B7).DOI:10.1029/2006JB004763.

[3] 孫自永，馬騰，馬軍，等.太原市地層空間異質性對地面沉降分布的影響[J].巖土力學,2007,28(2):399-405.

[4] FERRETTI A,PRATI C,ROCCA F. Permanent Scatterers in SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing, 2001,39(1):8-21.

[5] 吳宏安，張永紅，陳曉勇，等.基于小基線DInSAR的技術監測太原市2003～2009年地表形變場[J].地球物理學報，2011,54(3):673-680.

[6] 郭清海. 山西太原盆地孔隙地下水系統演化與相關環境問題成因分析[D].武漢：中國地質大學，2005.

[7] 王艷，廖明生，李德仁，等.利用長時間序列相干目標獲取地面沉降場[J].地球物理學報，2007，50(2)：598-604.