利用永久散射體雷達(dá)干涉技術(shù)進(jìn)行太原市地面沉降監(jiān)測

王 琪

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

一、引 言

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

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

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

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

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

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

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

(2)

(3)

為此，可以設(shè)定用以篩選PS點(diǎn)的相關(guān)系數(shù)γx為

(5)

式中，N為干涉圖的數(shù)量；γx即為判定PS點(diǎn)目標(biāo)的指標(biāo)。但僅僅如此還不夠，StaMPS結(jié)合信噪比、γx和PS概率進(jìn)行循環(huán)迭代來選取最終的PS點(diǎn)。

(6)

三、太原市地面沉降監(jiān)測試驗(yàn)

1． 試驗(yàn)區(qū)概況和試驗(yàn)數(shù)據(jù)

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

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

圖1 ASAR數(shù)據(jù)覆蓋范圍

2. 差分干涉處理

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

圖2 差分干涉圖

3. PS干涉處理

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

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

圖3 解纏相位圖

圖4 地表形變速率(mm/a) 圖5 形變速率標(biāo)準(zhǔn)差(mm/a)

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

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

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

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

圖7 太原不同歷史時(shí)段累積地面沉降演化(太原市地面沉降監(jiān)測報(bào)告，2004)

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

注：1956—2000年沉降數(shù)據(jù)來自文獻(xiàn)[7]。

5. 地面沉降時(shí)間演變特征

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

圖8 下元沉降時(shí)間序列

圖9 萬柏林沉降時(shí)間序列

圖10 吳家堡沉降時(shí)間序列

圖11 小店鎮(zhèn)沉降時(shí)間序列

四、結(jié)論與討論

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

參考文獻(xiàn)：

[1] 廖明生，林暉.雷達(dá)干涉測量—原理與信號處理基礎(chǔ)[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] 孫自永，馬騰，馬軍，等.太原市地層空間異質(zhì)性對地面沉降分布的影響[J].巖土力學(xué),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的技術(shù)監(jiān)測太原市2003～2009年地表形變場[J].地球物理學(xué)報(bào)，2011,54(3):673-680.

[6] 郭清海. 山西太原盆地孔隙地下水系統(tǒng)演化與相關(guān)環(huán)境問題成因分析[D].武漢：中國地質(zhì)大學(xué)，2005.

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