郭炳躍，何敏，劉建東

（1.江蘇省地質(zhì)礦產(chǎn)調(diào)查研究所，南京210018；2.河海大學(xué)測量工程系，南京210098）

利用PS－InSAR技術(shù)監(jiān)測南通市區(qū)地面沉降

郭炳躍1，何敏2，劉建東1

（1.江蘇省地質(zhì)礦產(chǎn)調(diào)查研究所，南京210018；2.河海大學(xué)測量工程系，南京210098）

以南通市區(qū)2006年至2007年的SAR影像為數(shù)據(jù)源，利用PS－InSAR技術(shù)對南通市區(qū)進(jìn)行地面沉降研究。結(jié)果顯示，南通市區(qū)存在多個沉降漏斗，但沉降量不大，大部分區(qū)域的線性沉降速率不超過11mm/a。

PS－InSAR技術(shù)；遙感監(jiān)測；地面沉降；南通市

1 引言

過渡開采含水層和城市建設(shè)引發(fā)的地面沉降災(zāi)害正影響著社會的發(fā)展。非常明顯的地質(zhì)現(xiàn)象是，大面積地區(qū)的地面在多年內(nèi)持續(xù)出現(xiàn)幾毫米甚至幾米的垂向位移。地面變形對城市的發(fā)展而言是非常重要的一個問題，地面沉降常常能夠破壞基礎(chǔ)設(shè)施，引發(fā)嚴(yán)重的經(jīng)濟損失。據(jù)估計，全世界有150多個城市出現(xiàn)了因地下水過量開采而引發(fā)的嚴(yán)重的地面沉降問題［1～4］。

本文將利用PS－InSAR技術(shù)技術(shù)方法與ERS和ENVISAT裝載的SAR傳感器獲得的圖像，重點分析發(fā)生南通市區(qū)的地面沉降。利用不同時期的干涉測量時間序列數(shù)據(jù)對地面沉降進(jìn)行分析，以顯示其在監(jiān)測地面沉降發(fā)展方面的潛力。

2 地質(zhì)背景

南通市位于新長江三角洲平原，地形平坦，地面高程4m左右，沿江一帶地勢較高，地面高程5m左右。地層為全新世亞粘土或亞砂土，成陸時間距今約2～4ka。覆蓋層為第四系的河流沖積與濱海相松散沉積物，松散沉積物厚度很大。從巖土工程角度來看，這些新近沉積物的壓縮性很高，處于欠固結(jié)狀態(tài)，問題最大。

南通市地下水類型主要是松散孔隙水，含水層比較多，各含水層的變化也比較復(fù)雜，由上到下分別為潛水含水層組、第Ⅰ承壓水、第Ⅱ承壓水、第Ⅲ承壓水、第Ⅳ承壓水和第Ⅴ承壓水。其中潛水含水層組由第四系全新統(tǒng)地層組成，含水層為亞砂土、粉砂，水位埋深0.3～3.5m，水位變化受降水控制，地層富水性較差。第Ⅰ承壓水由上更新統(tǒng)海、陸交替相沉積物組成，含水層為亞粘土與粉砂互層，富水性較差，單井涌水量50～60m3/d。第Ⅱ承壓水由中更新統(tǒng)地層組成，含水層以粉、細(xì)砂為主，厚15～85 m，含水層頂板埋深80～120m，變化較復(fù)雜。受沉積環(huán)境影響，富水性差異較大，一般單井涌水量500～1 000m3/d。第Ⅲ承壓水由下更新統(tǒng)河湖相地層組成，含水層為粉砂、細(xì)砂及含礫砂，單井涌水量一般為2 000m3/d左右，高的可達(dá)4 000m3/d，含水層富水性較好，是南通市主要開采含水層。第Ⅳ承壓水由上新統(tǒng)（N2）河湖相地層組成，含水層為粘性土夾粉砂、細(xì)砂，底部含少量礫石，厚度25～40m，富水性較好，單井涌水量一般為2 000m3/d［5］。

3 IPTA方法原理及處理流程

IPTA方法是基于PS技術(shù)的，可以說是PS技術(shù)思路的一種實現(xiàn)。假設(shè)有覆蓋同一地區(qū)的K＋1幅SAR影像，選擇其中一幅影像作為主影像，其余的所有影像都配準(zhǔn)并采樣到主影像像素空間，這樣K＋1幅SAR影像總共可形成K個干涉對。IPTA方法僅跟蹤成像區(qū)域內(nèi)雷達(dá)散射特性較為穩(wěn)定的目標(biāo)，而放棄那些失相關(guān)嚴(yán)重的分辨率單元，這些目標(biāo)幾乎不受失相關(guān)噪聲影響，即使在多年時間間隔的干涉對中，仍能保持較高的干涉相關(guān)性，因此被稱為永久散射體（PS）?；诠庾V離差小，后向散射強度大，不同時間SAR影像上的幅度差距小等原則，選取出候選PS點目標(biāo)，對K個干涉對上的PS點目標(biāo)逐一進(jìn)行干涉處理，再與外部DEM進(jìn)行差分，得到每個候選PS點的差分干涉相位。每個點目標(biāo)的差分干涉相位φdint包括：DEM高程誤差導(dǎo)致的相位項φdem，沿LOS方向的地表形變相位φdef，大氣影響相位φatm，及噪聲相位φnoise［6～9］。

其中：

式（2）中，λ為雷達(dá)波長；Rm為主影像成像時雷達(dá)至地面目標(biāo)的斜距；θ為雷達(dá)入射角；B⊥為干涉對垂直基線；δH為DEM高程改正。式（3）將地表形變分成線性和非線性兩部分；v為地表線性形變速率；t是干涉對時間基線。將式（2）、式（3）代入式（1）中，則

式（4）就是IPTA方法的相位模型，式中φres包含3種成分：大氣影響、地表非線性形變及噪聲，稱之為殘余相位。K1與DEM高程改正δH成正比，K2與LOS方向地表線性形變速率成正比，因此求出式（4）的系數(shù)K1、K2便可得到DEM高程改正及地表線性形變速率。

IPTA方法對干涉對垂直基線B⊥及時間基線t進(jìn)行回歸分析，求得公式（4）的系數(shù)，具體處理流程如圖1所示。首先將所有影像配準(zhǔn)、采樣到同一像元空間，基于配準(zhǔn)后的SLC影像，挑選出候選PS點目標(biāo)，當(dāng)研究區(qū)內(nèi)的SAR影像足夠多（＞25幅）時，以后向散射的穩(wěn)定性為準(zhǔn)則確定候選PS點目標(biāo)；當(dāng)SAR影像較少（＜25幅）時，以后向散射強度及光譜離差為準(zhǔn)則。然后計算各干涉對候選PS點目標(biāo)的差分干涉相位，對各PS點目標(biāo)的差分干涉相位時間序列進(jìn)行回歸分析，求得相位模型（公式（4））中的系數(shù)K1、K2。回歸模型計算得到的相位值標(biāo)準(zhǔn)差表征了回歸模型的精度，并根據(jù)此精度更新候選點目標(biāo)，剔除回歸模型精度較低的點目標(biāo)?；貧w分析的結(jié)果包括：高程改正，線性形變速率，解纏相位，回歸模型精度及殘余相位，利用這些結(jié)果對DEM高程、線性形變速率、基線長、大氣影響等不斷更新，使用更新后的差分干涉相位重新進(jìn)行回歸分析。如此迭代，不斷精化相位模型，直至得到滿意的結(jié)果，最終結(jié)果為PS點目標(biāo)改正后的高程、線性形變速率、形變歷史及精度信息［10，9］。

圖1 IPTA方法實現(xiàn)流程圖Fig.1 Flow chart of IPTA

4 南通地面沉降監(jiān)測數(shù)據(jù)

采用歐空局ENVISAT衛(wèi)星C波段雷達(dá)傳感器ASAR從2006年到2007年獲得的覆蓋南通地區(qū)的15幅SAR影像為源數(shù)據(jù)，這些影像均為單視復(fù)數(shù)影像（Single Look Complex，SLC）。根據(jù)使所有干涉對的時間基線及空間基線都盡量小的原則，選取2007年8月5日獲得的影像作為主影像，其余的14幅影像為從影像，形成14個干涉對。表1列出了所有影像的成像日期及相對于主影像的垂直基線B⊥和時間基線T。

實驗所需的ENVISAT衛(wèi)星精密軌道數(shù)據(jù)由荷蘭Delft大學(xué)空間研究中心（DEOS）提供，軌道徑向精度達(dá)到5～6cm，法向精度優(yōu)于20cm［12］，DEM數(shù)據(jù)使用SRTM－3數(shù)據(jù)。從SAR影像中截取以南通市區(qū)為中心，覆蓋范圍約為68km×57km的區(qū)域進(jìn)行處理，包括整個南通市區(qū)及通州，還覆蓋了海門、如皋、如東部分區(qū)域，圖2中方框所圍區(qū)域即為SAR影像截取范圍。

表1 南通地區(qū)ASAR影像數(shù)據(jù)和基線參數(shù)Table 1 ASAR image data and baseline parameters for Nantong

圖2 進(jìn)行IPTA處理的SAR影像覆蓋范圍Fig.2 SAR imaging coverage under IPTA processing

5 處理結(jié)果

5.1 IPTA處理結(jié)果

利用IPTA技術(shù)對表1中列出的ASAR數(shù)據(jù)進(jìn)行處理，得到研究區(qū)域（圖2）的地表形變速率，如圖3所示，一個圓點表示一個PS點目標(biāo)，圓點的顏色表示線性形變速率的大小，各種顏色表示的具體值域見圖示，形變量的單位是mm/a，負(fù)值表示沉降，正值表示抬升，顏色偏紅的地區(qū)即為沉降速率相對較大的地區(qū)。圖4為PS點目標(biāo)線性形變速率的分布直方圖，橫坐標(biāo)為線性形變速率區(qū)間，縱坐標(biāo)為落在各區(qū)間內(nèi)的PS點個數(shù)。

由圖4可知，在所選研究區(qū)域內(nèi)，共探測出25 047個PS點目標(biāo)，區(qū)域內(nèi)最大的形變速率為－16.821mm/a，平均形變速率為－4.155 6mm/a。從圖3中可看出PS點目標(biāo)多分布在市區(qū)（南通市區(qū)、通州市區(qū)、海門市區(qū)），而在其他鄉(xiāng)村范圍內(nèi)PS點則分布非常稀疏，且南通市區(qū)及海門地區(qū)的沉降速率相對較大。

5.2 PS點目標(biāo)線性形變速率空間分析結(jié)果

由于Kriging插值需要滿足正態(tài)分布及二階平穩(wěn)的前提假設(shè)，為保證插值的準(zhǔn)確性，在ArcGIS中選取PS點分布較密集、年沉降速率較大的南通市區(qū)進(jìn)行kriging插值。通過直方圖檢驗，南通市區(qū)內(nèi)PS點目標(biāo)的線性形變速率近似服從正態(tài)分布，且不存在全局離群值。南通市區(qū)內(nèi)PS點目標(biāo)的聚類Voronoi圖如圖5所示，灰色多邊形包圍的PS點目標(biāo)與其周圍所有PS點目標(biāo)的線性形變速率均不在同一級區(qū)間內(nèi)，則認(rèn)為此PS點為局部離群值，利用相鄰PS點目標(biāo)線性形變速率的平均值代替其原始值。為檢驗減弱離群值影響方法的有效性，將南通市區(qū)的PS點分割成兩部分，75%的點用來空間結(jié)構(gòu)建模及生成表面，25%的點用來驗證預(yù)測的質(zhì)量，稱為測試數(shù)據(jù)。分別對局部離群值處理前后PS點目標(biāo)的線性形變速率進(jìn)行內(nèi)插，測試數(shù)據(jù)的精度結(jié)果如表2所示，可見以相鄰點的均值代替局部離群值的方法有效減弱了局部離群值對空間插值的影響，提高了內(nèi)插精度。

表2 消除離群值影響前后插值結(jié)果比較Table 2 Comparison between interpolations before and after removing outliers

數(shù)據(jù)分析完成后，利用ArcGIS地統(tǒng)計分析向?qū)Чぞ哌M(jìn)行Kriging插值，插值結(jié)果如圖6所示，不同顏色代表了不同的線性形變速率區(qū)間，可見南通市區(qū)有多個沉降漏斗，沉降速率基本上都小于11 mm/a。圖中以星號為中心的沉降漏斗，覆蓋范圍及沉降量相對較大。區(qū)域1主要覆蓋了港閘區(qū)的南通船舶配套集中工業(yè)區(qū)，區(qū)域2為唐閘鎮(zhèn)街辦區(qū)域，區(qū)域4在港閘區(qū)政府區(qū)域，區(qū)域5覆蓋了以南通城市博物館為中心的區(qū)域，區(qū)域6是以紫瑯醫(yī)院為中心的區(qū)域。南通地區(qū)地面沉降主要是由于地下水開采引起地下水位下降，導(dǎo)致地層內(nèi)部壓力失衡，含水層本身及頂、低板粘性土層失水壓密而引起的。圖6中星號標(biāo)定的區(qū)域是工業(yè)園區(qū)，這些地方都是地下水開采最為集中的區(qū)域，所以沉降也就最為嚴(yán)重。

圖3 地表線性形變速率圖Fig.3 Surface linear deformation rate

圖4 PS點目標(biāo)線性形變速率分布直方圖Fig.4 Linear deformation rate distribution for PS point targets

圖5 利用ArcGIS生成的南通市區(qū)PS點目標(biāo)線性形變速率聚類Voronoi圖Fig.5 Linear deformation rate clustering Voronoi of PS point targets produced with ArcGIS

圖6 南通市區(qū)地表形變速率圖Fig.6 Surface deformation rate for Nantong City

6 結(jié)論

本文利用永久散射體雷達(dá)差分干涉測量技術(shù)提取地表線性形變速率的方法，以南通地區(qū)2006～2007年的SAR影像研究了南通市區(qū)的地面沉降。在此獲取南通市區(qū)沉降監(jiān)測數(shù)據(jù)的基礎(chǔ)上，利用ArcGIS地統(tǒng)計分析模塊對PS點的線性形變速率進(jìn)行空間分析，識別出其中的離群值并對其進(jìn)行處理，有效地減弱了離群值對空間插值的影響；選用Kriging方法對線性形變速率進(jìn)行內(nèi)插，獲得了南通市區(qū)高分辨的形變速率圖。研究結(jié)果顯示，南通市區(qū)存在多個沉降漏斗，但沒有出現(xiàn)沉降量非常大的沉降漏斗，在2006～2007年間大部分地區(qū)的線性沉降速率不超過11mm/a。

［1］Farr T G，Kobrick M.Shuttle radar topography mission produces a wealth of data［J］.EOS，TRANSACTIONS，American Geophysical Union，2000，81（48）：583－585.

［2］張詩玉，李陶，夏耶.基于InSAR技術(shù)的城市地面沉降災(zāi)害監(jiān)測研究［J］.武漢大學(xué)學(xué)報（信息科學(xué)版），2008，33（8）：850－853.

［3］羅小軍，黃丁發(fā)，劉國祥.基于永久散射體雷達(dá)差分干涉測量的城市地面沉降研究——以上海地面沉降監(jiān)測為例［J］.測繪通報，2009，（4）：4－8.

［4］羅海濱，何秀鳳.基于DInSAR方法檢測南京地表沉降的結(jié)果與分析［J］.高技術(shù)通訊，2008，18（4）：418－421.

［5］田立，錢宇紅.南通市地下水開采現(xiàn)狀及開發(fā)利用研究［J］.地下水，2008，30（3）：66－68.

［6］陳強，劉國祥，丁曉立，李永樹.永久散射體雷達(dá)差分干涉應(yīng)用于區(qū)域地表沉降探測［J］.地球物理學(xué)報，2007，50（3）：737－743.

［7］Ferretti A，Prati C，Rocca F.Permanent Scatterers in SAR Interferometry［J］.IEEETransactionsonGeoscienceandRemoteSensing，2001，39（1）：8－19.

［8］Ferretti A，Prati C，Rocca F.Non－linear subsidence rate estimation using permanent scatters in differential SAR interferometry［J］.IEEETransactionsonGeoscienceandRemote Sensing，2002，38（5）：2202－2212.

［9］Charles Werner，et al..Interferometric Point Target Analysis for Deformation Mapping［C］.IGARSS’03，Toulouse，F(xiàn)rance，2007，July：4362－4364.

［10］Qing Zhao，Hui Lin，et al..A Study of Ground Deformation in the Guangzhou Urban Area with Persistent Scatterer Interferometry［J］.Sensors，2009，9（1）：503－518.

［11］湯國安，楊昕.ArcGIS地理信息系統(tǒng)空間分析實驗教程［M］.北京：科學(xué)出版社，2006.

［12］陳強，劉國祥，李永樹.粗/精軌道數(shù)據(jù)對衛(wèi)星InSAR DEM精度影響的對比分析［J］.遙感學(xué)報，2006，10（4）：475－481.

MONITORING THE GROUND SUBSIDENCE IN NANTONG CITY WITH THE PS－INSAR TECHNOLOGY

Guo Bing－yue1，He Min2，Liu Jian－dong1
（1.Jiangsu Institute of Geology and Mineral Survey，Nanjing 210018，China；2.Department of Measuring Engineering，Hohai University，Nanjing 210098，China）

The ground subsidence in Nantong city is studied by using the PS－InSAR technology and the SAR image data from 2006to 2007.The result reveals several subsidence funnels with small settlements and the linear settling velocity of no more than 11mm/a in most regions.

PS－InSAR；remote sensing；ground subsidence；Nantong

P642.26

A

1006－4362（2011）04－0103－05

2011－03－09 改回日期：2011－09－26

郭炳躍（1980－ ），男，碩士，工程師，從事地質(zhì)環(huán)境與地質(zhì)災(zāi)害防治研究工作。