基于高分辨率TerraSAR－X影像的PSInSAR地表形變監(jiān)測(cè)

丁榮榮，徐 佳，林曉彬，許 康

(1.河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098;2.福建省港航管理局勘測(cè)中心，福州 350009;3.江蘇省測(cè)繪地理信息局，南京 210013)

0 引言

近幾十年來(lái)，合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)迅速發(fā)展，特別是合成孔徑雷達(dá)差分干涉測(cè)量(DInSAR)技術(shù)，由于具有全天時(shí)、全天候、覆蓋廣、分辨率高、穿透能力強(qiáng)以及主動(dòng)式遙感等優(yōu)點(diǎn)，被廣泛應(yīng)用于地震、地面沉降、火山和冰川活動(dòng)等的監(jiān)測(cè)［1］。但在利用常規(guī)DInSAR技術(shù)監(jiān)測(cè)長(zhǎng)時(shí)間的地表微小形變時(shí)，存在著時(shí)間、空間去相干和大氣延遲等影響，大大限制了DInSAR技術(shù)的發(fā)展和應(yīng)用。為了克服這些問(wèn)題，以永久散射體干涉測(cè)量(PSIn-SAR)為代表的時(shí)序分析方法不斷發(fā)展，并且以TerraSAR－X［2］為代表的高空間分辨率(3 m)數(shù)據(jù)的使用，在很大程度上促進(jìn)了PSInSAR技術(shù)的發(fā)展。本文利用高空間分辨率條帶模式的TerraSAR－X數(shù)據(jù)，采用干涉測(cè)量點(diǎn)目標(biāo)分析(interferometric point target analysis，IPTA)技術(shù)［3－4］對(duì)常州地區(qū)進(jìn)行地表形變監(jiān)測(cè)，其結(jié)果證明了TerraSAR－X數(shù)據(jù)在地表形變監(jiān)測(cè)中的優(yōu)勢(shì)。

1 研究區(qū)概況及數(shù)據(jù)選擇

1.1 研究區(qū)概況

常州市地處長(zhǎng)江下游南岸，位于江蘇省南部，長(zhǎng)江三角洲中心地帶，與上海、南京、杭州等市皆等距相鄰，與蘇州、無(wú)錫等市聯(lián)袂成片，介于N31°09'～32°04'、E119°08'～120°12'之間。常州市城鎮(zhèn)眾多，城市地區(qū)紡織業(yè)發(fā)達(dá)。自20世紀(jì)70年代始，由于大量集中開(kāi)采地下水，導(dǎo)致市區(qū)出現(xiàn)了嚴(yán)重的地面沉降;同時(shí)，隨著經(jīng)濟(jì)和城市化的發(fā)展，城市周邊地區(qū)的鄉(xiāng)鎮(zhèn)企業(yè)逐漸興起，地下水開(kāi)采進(jìn)一步加劇。據(jù)最新水準(zhǔn)點(diǎn)測(cè)量資料，在常州市五星鄉(xiāng)—青龍鄉(xiāng)一線的西南，湖塘鎮(zhèn)—戚墅埝一線以北的大面積范圍內(nèi)，1980年以來(lái)的地面沉降量超過(guò)了600 mm，其中超過(guò)800 mm的重度沉降區(qū)面積約3 417 km2，最大累計(jì)沉降量超過(guò)1 000 mm。

本文的實(shí)驗(yàn)區(qū)以常州市區(qū)為中心，覆蓋范圍約32 km×53 km，包括整個(gè)常州市中心及武進(jìn)地區(qū)，如圖1所示。

圖1 實(shí)驗(yàn)區(qū)范圍Fig.1 Range of experimental spot

1.2 數(shù)據(jù)選擇

本次選取了2011年10月—2013年7月期間37景TerraSAR－X［5］高分辨率條帶模式SAR數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。其影像均為單視復(fù)數(shù)影像(single look complex，SLC)。根據(jù)使所有干涉對(duì)的時(shí)間基線及空間基線最優(yōu)的原則，以2012年8月5日的影像為主影像，其余的36景影像為從影像，形成36個(gè)干涉對(duì)。表1列出了所有影像的成像日期及相對(duì)于主影像的垂直基線和時(shí)間基線。由表1可見(jiàn)，干涉對(duì)的基線都較小，適合 PSInSAR 分析［6－8］。

表1 37景Terrasar－X影像數(shù)據(jù)參數(shù)Tab.1 Parameters of 37 TerraSAR－X images

2 方法與實(shí)驗(yàn)

2.1 IPTA技術(shù)的原理和處理流程

PSInSAR 技術(shù)由意大利學(xué)者 Ferretti［9－10］等提出，它利用相位穩(wěn)定的點(diǎn)(PS)提取形變信息。由于其具有相干性幾乎不受時(shí)間基線和空間基線影響的特性，從而很好地避免了傳統(tǒng)DInSAR技術(shù)受時(shí)間和空間失相干的限制。點(diǎn)目標(biāo)干涉測(cè)量分析(IPTA)技術(shù)是一種改進(jìn)的PSInSAR技術(shù)，主要優(yōu)點(diǎn)是對(duì)低相干的地面仍能選取一定數(shù)量的PS點(diǎn)。只對(duì)提取的點(diǎn)進(jìn)行時(shí)間維和空間維分析，可以獲取長(zhǎng)時(shí)間序列的大范圍緩慢地表形變場(chǎng)。主要技術(shù)思想如下:

設(shè)研究區(qū)有N幅不同時(shí)相的SAR影像，選取其中一幅作為公共主影像，其余N－1幅為從影像，將從影像分別與主影像進(jìn)行配準(zhǔn)、重采樣和干涉處理，得到N－1幅干涉圖。利用外部DEM數(shù)據(jù)進(jìn)行差分干涉處理，就可得到N－1幅差分干涉相位，從而得到每個(gè)PS點(diǎn)上的N－1個(gè)時(shí)序差分干涉相位。其中，PS點(diǎn)解纏的干涉相位φunw包括DEM高程誤差引起的地形相位φtopo、沿LOS(line of sight)方向的地表形變相位φdef、大氣影響相位φatm及噪聲相位φnoise，IPTA方法的相位模型為

式中:K1=4π/(λRsinθ)，λ 為雷達(dá)波長(zhǎng)，R為雷達(dá)到地面目標(biāo)的斜距，θ為雷達(dá)入射角;B⊥為干涉對(duì)垂直基線;△H為DEM高程改正值;K2=4π/λ;t為干涉對(duì)時(shí)間基線;v為沿雷達(dá)視線方向的線性形變速率;φres是PS點(diǎn)的殘余相位，它包括大氣延遲相位、非線性形變相位和噪聲相位。

根據(jù)各部分相位在時(shí)間域和空間域的不同特性，利用最小二乘擬合、回歸分析等方法分離各不同參數(shù)，同時(shí)精化各參數(shù)，對(duì)殘余相位進(jìn)行分解，使用三維濾波方法，在空間上對(duì)噪聲相位濾波，在時(shí)間上對(duì)大氣相位和非線性形變相位濾波，進(jìn)而精化整個(gè)相位模型，求解出最佳的△H和v，即得到了PS點(diǎn)的DEM高程改正值和LOS方向上的線性形變速率。主要工作流程如圖2所示。

圖2 IPTA方法流程Fig.2 Flow chart of IPTA

2.2 數(shù)據(jù)處理和結(jié)果分析

利用IPTA方法對(duì)37景TerraSAR－X數(shù)據(jù)進(jìn)行處理，通過(guò)獲取PS點(diǎn)，在PS點(diǎn)上進(jìn)行時(shí)序分析，得到了相干點(diǎn)目標(biāo)在雷達(dá)視線(LOS)向的形變(圖3)。

圖3 TerraSAR－X地表形變平均速率(左)及其PS點(diǎn)分布(右)Fig.3 M ean velocity(left)derived from TerraSAR－X data and its PS points distribution(right)

從圖3(左)可以看出，主要沉降區(qū)域分布在武進(jìn)地區(qū)，其中鄒區(qū)鎮(zhèn)—牛塘鎮(zhèn)—廟橋鎮(zhèn)沉降連成一線，最大沉降速率為31.494 mm/a;從沉降圖上可清晰看到鄒區(qū)、牛塘附近均出現(xiàn)明顯沉降，其中牛塘下方地區(qū)沉降較為嚴(yán)重;在廟橋附近有大范圍沉降出現(xiàn);常州市中心城區(qū)比較平穩(wěn)，沒(méi)有明顯沉降?？傮w來(lái)看，常州市中心沉降比較緩慢，較大的沉降集中在武進(jìn)地區(qū)。

在數(shù)據(jù)的處理過(guò)程中，利用 Gamma軟件［11］生成了36幅干涉對(duì)影像，依據(jù)振幅與后向散射強(qiáng)度相結(jié)合的方法獲取PS點(diǎn)，從每幅干涉對(duì)影像中都得到539 927個(gè)PS點(diǎn)。將利用TerraSAR－X數(shù)據(jù)提取的高相干PS點(diǎn)分布部分放大(圖3(右))。通過(guò)PS點(diǎn)的位置可以看出，這些PS點(diǎn)大多密集分布于建筑物、道路、橋梁等具有穩(wěn)定電磁反射特性的物體上，這在某種程度上驗(yàn)證了探測(cè)到的PS點(diǎn)［12］的可靠性。

圖4 Envisat ASAR地表形變平均速率(左)及其PS點(diǎn)分布(右)Fig.4 M ean velocity(left)derived from Envisat ASAR data and its PS points distribution(right)

為了對(duì)比，本文利用2007—2010年期間的27景中等分辨率Envisat ASAR數(shù)據(jù)進(jìn)行處理，得到的常州市地表形變平均速率如圖4(左)所示。在設(shè)定 相同閾值的情況下，對(duì)比圖3(右)和圖4(右)可以看出，應(yīng)用TerraSAR－X數(shù)據(jù)提取的PS點(diǎn)覆蓋了整個(gè)沉降區(qū)域，而利用Envisat ASAR數(shù)據(jù)提取的PS點(diǎn)分布相對(duì)稀疏，不能全面覆蓋整個(gè)區(qū)域。從圖3(右)可以看出，TerraSAR－X PS點(diǎn)的分布主要呈現(xiàn)聚束狀，聚集式分布于建筑物、道路等地物上，可以更好地描述地物細(xì)節(jié)，這是Envisat ASAR數(shù)據(jù)在PSInSAR處理過(guò)程中無(wú)法得到的。圖3(左)上呈現(xiàn)的地表形變平均速率沒(méi)有圖4(左)的大的主要原因是TerraSAR－X數(shù)據(jù)的時(shí)間基線只有660 d，所以沉降速率較Envisat ASAR數(shù)據(jù)的要小。

進(jìn)一步將沉降較大的區(qū)域局部放大，如圖5所示。武進(jìn)地區(qū)有5處明顯沉降區(qū)(A—E)，其中，

圖5 武進(jìn)地區(qū)地表形變速率圖Fig.5 M ean velocity around W ujin area(TerraSAR－X)

A和B區(qū)沉降速率最大，為14.1 mm/a，清楚可見(jiàn)沉降區(qū)PS點(diǎn)分布在建筑物房頂，主要沉降區(qū)多為紡織廠、化工廠分布的地方;在2011—2013年間，D區(qū)為沉降嚴(yán)重地區(qū)，平均沉降速率達(dá)到29.0 mm/a;C區(qū)和E區(qū)處沉降量較大，最大沉降速率分別達(dá)到20.0 mm/a和19.8 mm/a。從2007—2011 年間 Envisat ASAR數(shù)據(jù)提取的沉降速率圖上并沒(méi)有發(fā)現(xiàn)C區(qū)大興化工廠附近的建筑物，而在2011—2013年間TerraSAR－X數(shù)據(jù)提取的沉降速率圖中發(fā)現(xiàn)此處最大沉降速率達(dá)到18.0 mm/a(圖5)，多數(shù)廠房周?chē)嬖诤芏喑两迭c(diǎn)，這也說(shuō)明了過(guò)度開(kāi)采地下水是造成地面沉降的主要原因。

2.3 沈海高速公路沿線地面沉降監(jiān)測(cè)

沈海高速公路橫穿本文實(shí)驗(yàn)區(qū)，長(zhǎng)約50 km。為了更加清楚地說(shuō)明該高速公路附近的地表沉降情況，將其分為Ⅰ，Ⅱ，Ⅲ3段(圖6)。

圖6 沈海高速公路沿線地面沉降Fig.6 Land subsidence along Shenhai highway

TerraSAR－X圖像清晰地顯示出高速公路及其沿線的地物細(xì)節(jié)［13］，PS點(diǎn)聚集分布于高速公路及其沿線周?chē)姆课萆?。在Ⅰ段，高速公路沉降速率較小，為6～14 mm/a;在Ⅱ段，沉降速率逐漸變大，最大達(dá)24.6 mm/a。該區(qū)也正是圖5中的嚴(yán)重沉降區(qū)(D區(qū))，區(qū)段內(nèi)坐落著的許多紡織、化工、印刷等大型工廠都是用水大戶(hù)，正是這些工廠過(guò)度開(kāi)采地下水導(dǎo)致了嚴(yán)重的地面沉降;在Ⅲ段，自西向東沉降速率變小，并逐漸趨于穩(wěn)定。

3 與水準(zhǔn)數(shù)據(jù)的對(duì)比分析

為進(jìn)一步驗(yàn)證實(shí)驗(yàn)結(jié)果精度，在沉降區(qū)附近找到相近的水準(zhǔn)點(diǎn)進(jìn)行比較分析。本文采用歐空局開(kāi)展的PSIC4(Persistent Scatterers Interferometry Codes Cross－Comparison and Certification for Long－Term Differential Interferometry)項(xiàng)目中使用到的2個(gè)指標(biāo)來(lái)衡量測(cè)量精度［14］，水準(zhǔn)數(shù)據(jù)與PSInSAR監(jiān)測(cè)數(shù)據(jù)互差的均值為

水準(zhǔn)數(shù)據(jù)與PSInSAR監(jiān)測(cè)數(shù)據(jù)互差的均方差為

實(shí)驗(yàn)區(qū)域共有70個(gè)同步水準(zhǔn)點(diǎn)，如圖7(a)所示。根據(jù)式(2)和式(3)計(jì)算得到的水準(zhǔn)與PSIn-SAR監(jiān)測(cè)數(shù)據(jù)互差的平均誤差和均方差分別為0.54 mm 和 ±6.40 mm。

圖7 水準(zhǔn)點(diǎn)分布Fig.7 Distribution of bench mark

為了更進(jìn)一步檢驗(yàn)測(cè)量精度，對(duì)沉降明顯的地區(qū)采用沉降區(qū)附近10個(gè)水準(zhǔn)點(diǎn)單獨(dú)比較的方式(圖7(b))，選擇距離該水準(zhǔn)點(diǎn)100m范圍內(nèi)相干點(diǎn)的平均沉降值作為評(píng)定標(biāo)準(zhǔn)(表2)。

表2 PS點(diǎn)沉降誤差Tab.2 Settlement of PS error(mm/a)

由表2可以看出，通過(guò)與區(qū)域內(nèi)分布的水準(zhǔn)點(diǎn)測(cè)量數(shù)據(jù)的比較，P1—P10個(gè)點(diǎn)中有9個(gè)點(diǎn)的誤差都在±3 mm/a以?xún)?nèi)，只有1個(gè)點(diǎn)超過(guò)3 mm/a，但也在±4 mm/a以?xún)?nèi)的。表中結(jié)果與水準(zhǔn)測(cè)量結(jié)果具有較好的一致性，充分說(shuō)明了利用基于TerraSAR－X高分辨率條帶模式SAR數(shù)據(jù)開(kāi)展PSInSAR測(cè)量地表形變的準(zhǔn)確性和可行性。

4 結(jié)論

1)利用TerraSAR－X高分辨率條帶模式SAR數(shù)據(jù)得到了2011年10月—2013年7月期間常州市地表形變平均速率圖，最大沉降速率達(dá)31.494 mm/a，證明常州地區(qū)存在嚴(yán)重的地面沉降。

2)通過(guò)對(duì)比TerraSAR－X與Envisat ASAR的監(jiān)測(cè)結(jié)果發(fā)現(xiàn)，TerraSAR－X提取的PS點(diǎn)分布密度更高，可以更清晰地探測(cè)到紡織廠附近的地物細(xì)節(jié)，分析由于過(guò)度開(kāi)采地下水造成的地面沉降原因，體現(xiàn)了高分辨率TerraSAR－X數(shù)據(jù)的特點(diǎn)，在城市地表形變監(jiān)測(cè)中具有廣闊的應(yīng)用前景。

3)結(jié)合TerraSAR－X數(shù)據(jù)得到的沈海高速公路常州段地表形變信息，分析高速自身沉降及沿線地面沉降，發(fā)現(xiàn)主要沉降地區(qū)正是由于集中開(kāi)設(shè)工廠、過(guò)度開(kāi)采地下水導(dǎo)致的，驗(yàn)證了TerraSAR－X數(shù)據(jù)在監(jiān)測(cè)人工線狀地物方面有很好的應(yīng)用前景。

4)處理過(guò)程中體現(xiàn)了IPTA技術(shù)的顯著優(yōu)勢(shì)。由于它是利用具有穩(wěn)定特性的PS點(diǎn)進(jìn)行形變分析的，有效地避免了常規(guī)DInSAR中存在時(shí)間去相干和大氣效應(yīng)的影響。

5)常州市沉降嚴(yán)重區(qū)域主要集中在武進(jìn)地區(qū)，通過(guò)與水準(zhǔn)數(shù)據(jù)比較，兩者具有很好的一致性，證明了IPTA技術(shù)監(jiān)測(cè)結(jié)果的正確性。

［1］ 廖明生，裴媛媛，王寒梅，等.永久散射體雷達(dá)干涉技術(shù)監(jiān)測(cè)上海地面沉降［J］.上海國(guó)土資源，2012(3):5－10.Liao M S，Pei Y Y，Wang H M，et al.Subsidence monitoring in Shanghai using the PSInSAR technique［J］.Shanghai Land ＆ Resources，2012(3):5－10.

［2］ Even M，Schunert A，Schulz K，et al.Atmospheric phase screen－estimation for PSInSAR applied to TerraSAR－X high resolution spotlight－ data［C］//2010 IEEE International Geoscience and Remote Sensing Symposium(IGARSS).Honolulu，HI:IEEE，2010:2928－2931.

［3］ 周紅滿，胡金星，柳 想，等.基于PS技術(shù)的深圳市地面沉降監(jiān)測(cè)［J］.測(cè)繪工程，2013，22(2):70－74.Zhong H M，Hu JX，Liu X，etal.Ground deformation in the Shenzhen city using permanent scatterer technique［J］.Engineering of Surveying and Mapping，2013，22(2):70－74.

［4］ 郭炳躍，何 敏，劉建東.利用PS－InSAR技術(shù)監(jiān)測(cè)南通市區(qū)地面沉降［J］.地質(zhì)災(zāi)害與環(huán)境保護(hù)，2011，22(4):103－107.Guo B Y，He M，Liu J D.Monitoring the ground subsidence in Nantong city with the PS－ InSAR technology［J］.Journal of Geological Hazards and Environment Preservation，2011，22(4):103－107.

［5］ 李永生，張景發(fā)，羅 毅，等.利用高分辨率聚束模式TerraSARX影像的PSInSAR監(jiān)測(cè)地表變形［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2012，37(12):1452－1455，1514.Li Y S，Zhang JF，Luo Y，et al.Monitoring land deformation using PSInSAR with TerraSAR－X high resolution spotlight SAR images［J］.Geomatics and Information Science of Wuhan University，2012，39(12):1452－1455，1514.

［6］ 何秀鳳，仲海蓓，何 敏.基于PS－InSAR和GIS空間分析的南通市區(qū)地面沉降監(jiān)測(cè)［J］.同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版，2011，39(1):129－134.He X F，Zhong H B，He M.Ground subsidence detection of Nantong city based on PS－ InSAR and GIS spatial analysis［J］.Journal of TongjiUniversity:Natural Science，2011，39(1):129－134.

［7］ 范景輝，李 梅，郭小方，等.基于PSInSAR方法和ASAR數(shù)據(jù)監(jiān)測(cè)天津地面沉降的試驗(yàn)研究［J］.國(guó)土資源遙感，2007，19(4):23－27，124.doi:10.6046/gtzyyg.2007.04.05.Fan JH，LiM，Guo X F，et al.A preliminary study of the subsidence in Tianjin area using ASAR images based on PSInSAR technique［J］.Remote Sensing for Land ＆ Resources，2007，19(4):23－27，124.doi:10.6046/gtzyyg.2007.04.05.

［8］ Strzelczyk J，Porzycka S，Lesniak A.Analysis of ground deformations based on parallel geo－statistical computations of PSInSAR data［C］//Proceedings of the 17th International Conference on Geoinformatics.Fairfax，VA:IEEE，2009:1－6.

［9］ Ferretti A，Prati C，Rocca F L.Permanent scatterers in SAR interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39(1):8－20.

［10］ Ferretti A，F(xiàn)umagalli A，Novali F，et al.A new algorithm for processing interferometric data－ stacks:SqueeSAR［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49(9):3460－3470.

［11］ 姜 江，宋偉東.基于GAMMA的雷達(dá)影像干涉測(cè)量分析［J］.交通科技與經(jīng)濟(jì)，2009，11(1):45－46，53.Jiang J，Song W D.The study on Radar interferometry with GAMMA［J］.Technology ＆ Economy in Areas of Communications，2009，11(1):45－46，53.

［12］ Porzycka－Strzelczyk S，Lesniak A.The preliminary temporal analysis of ground deformations in the area of Dabrowski Coal Basin(south Poland)［C］//IEEE International Geoscience and Remote Sensing Symposium(IGARSS).Honolulu，HI:IEEE，2010:903－905.

［13］ 林 琿，陳富龍，江利明，等.多基線差分雷達(dá)干涉測(cè)量的大型人工線狀地物形變監(jiān)測(cè)［J］.地球信息科學(xué)學(xué)報(bào)，2010，12(5):718－725.Lin H，Chen F L，Jiang LM，et al.Preliminary research on large－scale man－made linear features deformation monitoring using multi－ baseline differential SAR interferometry［J］.Journal of Geo－ Information Science，2010，12(5):718－725.

［14］ 葛大慶.區(qū)域性地面沉降InSAR監(jiān)測(cè)關(guān)鍵技術(shù)研究［D］.北京:中國(guó)地質(zhì)大學(xué)，2013.Ge DQ.Research on the Key Techniques of SAR Interferometry for Regional Land Subsidence Monitoring［D］.Beijing:China University of Geoscience，2013.