彭　鵬， 楊紅磊

(1.安徽省地質(zhì)調(diào)查院，合肥　230001； 2.中國(guó)地質(zhì)大學(xué)(北京)，北京　100083)

?

基于改進(jìn)相干點(diǎn)目標(biāo)技術(shù)的阜陽(yáng)市地面沉降調(diào)查

彭鵬1， 楊紅磊2

(1.安徽省地質(zhì)調(diào)查院，合肥230001； 2.中國(guó)地質(zhì)大學(xué)(北京)，北京100083)

提出了一種改進(jìn)的相干點(diǎn)目標(biāo)時(shí)序InSAR技術(shù)方法。該方法基于子視相干值識(shí)別小數(shù)據(jù)集的相干點(diǎn)目標(biāo)，采用趨勢(shì)項(xiàng)濾波方法去除大尺度殘差相位，并以阜陽(yáng)市2012—2013年間10景Radarsat-2 Wide模式數(shù)據(jù)為例，獲得了高空間密度的時(shí)序高相干點(diǎn)目標(biāo)，提高了結(jié)果的可靠性，獲取的沉降區(qū)數(shù)據(jù)與前人實(shí)測(cè)資料及野外調(diào)查的沉降現(xiàn)象相一致，為開展阜陽(yáng)市地面沉降風(fēng)險(xiǎn)分析和經(jīng)常性預(yù)警監(jiān)測(cè)提供技術(shù)方法和示范應(yīng)用。

相干點(diǎn)目標(biāo)技術(shù)； 地面沉降； 濾波； Radarsat-2

0　引言

InSAR(Interferometric SAR)技術(shù)是利用不同時(shí)期2景SAR影像間的相位差來(lái)獲取地面沉降信息，具有全天時(shí)、全天候和高形變監(jiān)測(cè)精度的特點(diǎn)[1]。利用該技術(shù)已在我國(guó)主要地表形變區(qū)和重大工程區(qū)開展了大量地表形變摸底調(diào)查工作[2]。相干點(diǎn)目標(biāo)技術(shù)(coherent point target，CPT)是通過對(duì)時(shí)間序列中相位保持穩(wěn)定的點(diǎn)進(jìn)行建模分析，獲得高精度的地表形變信息，被廣泛應(yīng)用在地面沉降[3-4]、地震形變[5]、火山運(yùn)動(dòng)[6]以及山體滑坡[7]等領(lǐng)域。目前識(shí)別高質(zhì)量相干點(diǎn)目標(biāo)的常用方法有相干系數(shù)閾值法[8]、振幅離差閾值法[8]、相位離差閾值法[9]和信噪比[10]等。然而，上述方法都基于大數(shù)據(jù)量(>30景)，從統(tǒng)計(jì)的思想選擇相干點(diǎn)目標(biāo)，不適合小數(shù)據(jù)集的情況。雖然相干點(diǎn)目標(biāo)分析(interferometric point target analysis，IPTA)技術(shù)采用了子視相干值來(lái)識(shí)別小數(shù)據(jù)集的相干點(diǎn)目標(biāo)[3]，但在穩(wěn)健性方面仍有不足。針對(duì)這一問題，提出了利用子視相干值離差識(shí)別相干點(diǎn)目標(biāo)的方法。數(shù)據(jù)源采用Radarsat-2數(shù)據(jù)， Radarsat-2 Wide模式具有監(jiān)測(cè)范圍大的特點(diǎn)，適用于地面沉降普查和經(jīng)常性日常監(jiān)測(cè)預(yù)警。但是由于其定軌精度較低，采用常規(guī)的FFT[11]、GCP[12]和二項(xiàng)式擬合方法[13]也難以消除軌道殘余誤差，基于該誤差在空間上表現(xiàn)為大尺度低頻信號(hào)[14]，提出了采用趨勢(shì)項(xiàng)濾波方法糾正其影響。

1　研究方法

改進(jìn)的相干點(diǎn)目標(biāo)技術(shù)共包含識(shí)別相干點(diǎn)目標(biāo)、生成差分干涉圖、相位解纏和估計(jì)形變值4個(gè)部分。

1.1識(shí)別相干點(diǎn)目標(biāo)

首先從若干景SAR影像中，依據(jù)總體相關(guān)性最大原則[14]確定參考影像，然后把其他影像配準(zhǔn)到該影像的成像幾何，配準(zhǔn)精度優(yōu)于0.2個(gè)像元。

與分布式散射體(如農(nóng)田、森林等)相比，點(diǎn)狀目標(biāo)(如建筑物、船舶等)可視為由一個(gè)或者多個(gè)“主散射體”構(gòu)成的反射體，具有強(qiáng)振幅值，不受斑點(diǎn)噪聲的影響，在單視復(fù)數(shù)影像中具有確定的相位值，在頻率域中表現(xiàn)為一種獨(dú)特的光譜特征，可以依據(jù)該屬性來(lái)識(shí)別相干點(diǎn)目標(biāo)。通過在方位向或者距離向降低處理器帶寬，將頻譜分割為若干部分，即為子視分解[15]。由于分布式散射體內(nèi)各獨(dú)立散射體回波的隨機(jī)性，當(dāng)任意子視影像的頻譜不重疊時(shí)，其斑點(diǎn)噪聲像素之間是失相干的，可以采用相關(guān)系數(shù)γ[16]確定相干散射體，即

(1)

式中：X1、X2分別為子視光譜；*代表取共軛；〈〉代表空間求平均。依據(jù)式(1)采用單景SAR影像就可以有效地確定相干點(diǎn)目標(biāo)，但是其中包含了汽車、輪船等移動(dòng)目標(biāo)。在CPT處理中，僅需要在整個(gè)時(shí)間序列內(nèi)相對(duì)穩(wěn)定的相干點(diǎn)目標(biāo)，故提出了采用光譜離差DC衡量時(shí)間序列點(diǎn)的光譜穩(wěn)定性，剔除移動(dòng)相干點(diǎn)目標(biāo)，即

(2)

式中δC和mC分別為時(shí)間序列內(nèi)光譜相關(guān)值的標(biāo)準(zhǔn)差和均值。

1.2生成差分干涉圖

為了弱化時(shí)空失相干的影響，采用短時(shí)空基線的原則確定干涉對(duì)。對(duì)于不同的波長(zhǎng)，設(shè)置不同的基線閾值，一般建議設(shè)為極限基線的1/3。干涉組合確定以后，對(duì)相干點(diǎn)目標(biāo)進(jìn)行干涉，得到每個(gè)點(diǎn)的干涉相位。再依據(jù)軌道和地形數(shù)據(jù)模擬平地效應(yīng)和地形相位，并從干涉圖中去除，得到每個(gè)點(diǎn)的差分干涉相位。此時(shí)的差分干涉相位包含軌道殘余相位、地形殘余相位、形變相位、大氣延遲相位和噪聲相位。除了地形殘余相位和噪聲相位外，其他部分都具有空間相關(guān)性。

1.3相位解纏

由1.2節(jié)獲得的相干點(diǎn)目標(biāo)x的一次差分時(shí)間序列干涉相位為

(3)

從相干點(diǎn)目標(biāo)集中選擇一個(gè)具有穩(wěn)定高相干值且位于地表比較穩(wěn)定區(qū)域的點(diǎn)作為計(jì)算參考點(diǎn)y，x和y再次做差分，得到二次差分干涉相位，表達(dá)式為

(4)

由于相鄰點(diǎn)之間的高程差是垂直基線的函數(shù)，相鄰點(diǎn)之間的形變速率與時(shí)間間隔相關(guān)，式(4)也可以表示為

(5)

式中： λ為雷達(dá)波長(zhǎng)； βk為高程常數(shù)； Δh(x,y)為相鄰點(diǎn)高程差； Tk為相鄰時(shí)刻SAR時(shí)間間隔； Δv(x,y)為時(shí)序形變速率。

采用解空間搜索法求得高程差與形變速率，恢復(fù)相鄰點(diǎn)之間在每個(gè)時(shí)序差分干涉對(duì)上的真實(shí)相位，實(shí)現(xiàn)相干點(diǎn)目標(biāo)時(shí)間域相位解纏。從原始差分干涉相位中減去高程和形變相位后，得到的殘差相位包含軌道殘余誤差、大氣延遲、非線性形變和噪聲相位。由于軌道殘余誤差相位和大氣延遲相位在空間上具有強(qiáng)相關(guān)性，對(duì)解纏后的殘差相位采用趨勢(shì)項(xiàng)濾波方法進(jìn)行分離，該方法由多視平均、空間濾波和空間插值3部分組成。具體流程如下：

(1) 多視平均。假設(shè)M×N的矩陣D對(duì)應(yīng)的矩陣元素為di,j，采用M×N的多視窗口，得到多視平均后的結(jié)果為

(6)

由式(6)可知，對(duì)數(shù)據(jù)進(jìn)行多視平均處理，去除噪聲的同時(shí)，還保持了趨勢(shì)項(xiàng)相位的形態(tài)。

(2) 空間濾波。由于非線性形變相位和軌道殘余相位、大氣延遲相位一樣，在空間上都具有較強(qiáng)相關(guān)性，但是空間尺度較小，當(dāng)采用視數(shù)較大的窗口多視處理后，非線性形變表現(xiàn)為噪聲點(diǎn)，可以采用空間均值、空間中值等濾波方法進(jìn)行去除。

(3) 空間插值。空間濾波后軌道殘余相位和大氣延遲相位的空間分辨率和原始圖像不一致，需要采用空間插值方法(反距離加權(quán)、樣條函數(shù)插值和Kriging插值等)恢復(fù)到原始分辨率。

1.4估計(jì)形變值

1.3節(jié)的處理不僅實(shí)現(xiàn)了差分干涉相位的相位解纏，而且可以分離出軌道殘余相位和大氣延遲相位。從得到的解纏差分干涉圖中剔除軌道殘余相位和大氣延遲相位后，剩下的差分干涉相位只包含了線性形變、非線性形變和高程殘差信息，采用SVD[17]分解便可以獲得時(shí)間序列形變的最小二乘解。

2　阜陽(yáng)市地面沉降監(jiān)測(cè)應(yīng)用

阜陽(yáng)市位于安徽省西北部，黃淮海沖積平原南端，地形平坦開闊。境內(nèi)西北高、東南低，高程為20～39 m，坡降為1/8 000～1/9 000。利用InSAR數(shù)據(jù)尚未開展過系統(tǒng)的地面沉降調(diào)查工作。根據(jù)阜陽(yáng)市一、二等水準(zhǔn)測(cè)量結(jié)果插值繪制的1998—2008年10 a間的沉降等值線圖，沉降中心位于三角游園一帶，沉降速率可達(dá)-40 mm/a[18]，近些年由于限制地下水開采，沉降變化緩慢。

數(shù)據(jù)源為Radarsat-2 Wide模式數(shù)據(jù)，空間分辨率距離向?yàn)?1.8 m，方位向?yàn)?.2 m，極化方式為VV，獲取時(shí)間為2012年2月4日—2013年1月29日，共10景。高程數(shù)據(jù)采用SRTM DEM。

依照總體相干值選擇2012年5月10日為參考影像，把其他影像配準(zhǔn)到該影像的成像幾何。分別對(duì)每景SAR影像進(jìn)行子視分解，求取子視相干值，設(shè)置時(shí)序子視相干均值閾值為0.37，相干離差閾值為1.15， 共識(shí)別 673 456 個(gè)相干點(diǎn)目標(biāo)， 主要

分布在阜陽(yáng)市區(qū)及郊區(qū)的居民地，農(nóng)田區(qū)域也有零星分布，實(shí)地調(diào)查得知，主要是農(nóng)田中的灌溉設(shè)施。設(shè)置空間基線閾值為300 m，時(shí)間基線為200 d，共確定30對(duì)干涉對(duì)，干涉對(duì)組合如圖1所示。

圖1　干涉對(duì)基線組合示意圖Fig.1　Schematic diagram of interferometer baseline combination

由2012年2月4日與2012年4月16日數(shù)據(jù)組合干涉對(duì)生成的差分干涉圖如圖2所示。

圖2　2012-02-04與2012-04-16影像的差分干涉圖Fig.2　Differential interferograms of images on 2012-02-04 and 2012-04-16

從圖2(a)可以看出，由于軌道定軌精度低，生成的基線誤差較大，造成初始差分干涉圖中殘余的平地效應(yīng)和高程誤差較大，甚至淹沒了形變值對(duì)應(yīng)的相位。采用GCP糾正雖然能夠消除一部分基線誤差(圖2(b))，但是仍不能完全消除其影響。圖2(c)為采用趨勢(shì)項(xiàng)濾波方法得到的改進(jìn)后差分干涉圖(多視視數(shù)為10×10，空間濾波方法為均值濾波，濾波窗口為7×7，插值方法為Kriging插值)，可以明顯發(fā)現(xiàn)已基本消除大尺度殘余相位的影響，能夠識(shí)別出形變區(qū)域，說(shuō)明本文方法是可行的。

圖3　阜陽(yáng)市2012—2013年間地面沉降速率Fig.3　Land subsidence rates in Fuyang city from 2012 to 2013

對(duì)所有干涉對(duì)成功相位解纏后，采用SVD分解得到沉降速率如圖3所示。

由圖3可知，阜陽(yáng)市形變主要集中在阜陽(yáng)市區(qū)，存在一個(gè)較明顯的沉降濾斗，市區(qū)南部村鎮(zhèn)也有沉降現(xiàn)象。該漏斗以奎星樓為中心，東至三角洲公園、西至中南現(xiàn)代城、南至區(qū)政府、北至青穎公園，覆蓋范圍約14 km2，最大形變速率為21 mm/a，分布范圍與以往沉降中心一致。由于2008年后該地區(qū)限制開采地下水，沉降有所減緩，本次InSAR調(diào)查速率比參考文獻(xiàn)[18]介紹的有所減小是合理的。

3　結(jié)果驗(yàn)證與機(jī)理分析

通過分析和野外調(diào)查發(fā)現(xiàn)，阜陽(yáng)市和地面沉降直接相關(guān)的現(xiàn)象有抽水井管道彎曲、抬升和橋梁拉伸、錯(cuò)位。

3.1井管彎曲、抬升地表形變現(xiàn)象

開采地下水的深井井管抬升及井臺(tái)地面開裂等地表變形現(xiàn)象如圖4和圖5所示。

圖4 自來(lái)水公司深井井管彎曲Fig.4 BendedpipeofboreintheFuyangtapwatercompany 圖5 井臺(tái)遭受破壞Fig.5 Damagedwellplatform

具體統(tǒng)計(jì)數(shù)據(jù)如表1。

表1　由地面沉降引發(fā)變形井的統(tǒng)計(jì)數(shù)據(jù)

由圖4、圖5和表1可以看到，位于InSAR地面沉降區(qū)的服裝廠深井(現(xiàn)屬自來(lái)水公司)，井臺(tái)較周圍地面有明顯抬升，井臺(tái)地面開裂，裂縫呈放射狀； 該深井井管彎曲明顯，可達(dá)25 cm。該形變主要由于地下水超采造成淺部地層沉降所致，而井管末端所在深部地層較為穩(wěn)定，使最上端的井臺(tái)與井管保持在原始高度，產(chǎn)生類似井管躥升的效果。野外調(diào)查點(diǎn)均位于InSAR形變速率圖中存在明顯沉降的區(qū)域。

3.2橋梁拉伸、錯(cuò)位地表變形現(xiàn)象

阜陽(yáng)市跨度較大的大橋、公路等都出現(xiàn)了橋面開裂、拉伸及錯(cuò)位等現(xiàn)象(如圖6、7所示)，經(jīng)過修繕，已恢復(fù)使用。

圖6 橋梁錯(cuò)位Fig.6 Bridgedislocation圖7 橋梁拉伸Fig.7 Bridgestretching

通過對(duì)InSAR沉降區(qū)進(jìn)行野外實(shí)地調(diào)查，該現(xiàn)象都位于沉降區(qū)范圍，可能為不均勻沉降所致，應(yīng)作為地面沉降重點(diǎn)調(diào)查區(qū)。

3.3地面沉降形成機(jī)理分析

3.3.1地層特征

研究區(qū)區(qū)域地層屬華北地層區(qū)淮河地層分區(qū)阜陽(yáng)地層小區(qū)，發(fā)育有新元古界霍邱群和五河群，古生界寒武系、奧陶系、石炭系和二疊系，中生界白堊系，新生界古近系、新近系和第四系地層。研究區(qū)內(nèi)主要為新生代松散層覆蓋，厚800～900 m，其中第四系厚130～150 m，松散沉積巨厚。

3.3.2水文地質(zhì)特征分析

阜陽(yáng)市地下水類型為單一松散巖類孔隙水，根據(jù)地下水賦存介質(zhì)的性質(zhì)、類別和組合的不同，自上而下劃分為淺層地下水和深層地下水[18-19]。淺層地下水賦存于50 m以淺的全新世、晚更新世地層中，與大氣降水、地表水關(guān)系密切，按上下關(guān)系可稱其為第一含水層組(淺層)； 深層地下水賦存于50 m以下的地層中，與大氣降水、地表水關(guān)系不密切。根據(jù)水文地質(zhì)結(jié)構(gòu)和目前開采現(xiàn)狀，將深層地下水進(jìn)一步劃分為2個(gè)含水層組，即第二含水層組(中深層，埋深50～150 m)和第三含水層組(深層，埋深150～500 m)。隨著城市化的發(fā)展，深井越打越多，越打越深，開采量不斷增加，致使城區(qū)及周邊深層地下水超采，這也成為研究區(qū)發(fā)生地面沉降的主要原因。

1.第一含水層組水位及流向; 2.第二含水層組水位及流向; 3.第三含水層組水位及流向; 4.砂層; 5.黏性土層; 6.壓密釋水及越流層水流向; 7.地表水圖8　阜陽(yáng)市含水層模型剖面示意圖Fig.8　Schematic section of aquifer model in Fuyang

阜陽(yáng)市含水層模型剖面示意圖如圖8所示。

4　結(jié)論

本文分別采用子視相干點(diǎn)目標(biāo)識(shí)別方法和趨勢(shì)項(xiàng)濾波方法，對(duì)傳統(tǒng)相干點(diǎn)目標(biāo)技術(shù)中的相干點(diǎn)目標(biāo)識(shí)別和大尺度殘余誤差剔除2方面進(jìn)行改進(jìn)。以阜陽(yáng)市2012—2013年間10景Radarsat-2 Wide模式數(shù)據(jù)為例進(jìn)行實(shí)驗(yàn)，獲得該市2012—2013年間的形變速率圖，認(rèn)為該市存在一個(gè)較為明顯的形變漏斗。結(jié)合形變區(qū)實(shí)地調(diào)查，初步分析了形變產(chǎn)生的影響及原因，為該市形變監(jiān)測(cè)網(wǎng)的布設(shè)和地面沉降防治提供了決策依據(jù)。

[1] Gabriel A K,Goldstein R M.Crossed orbit interferometry:Theory and experimental results from SIR-B[J].International Journal of Remote Sensing,1988,9(8):857-872.

[2] 范景輝,燕云鵬,葛大慶,等.全國(guó)地表形變遙感地質(zhì)(InSAR)調(diào)查技術(shù)指南[M].北京:地質(zhì)出版社,2015.

[3] Werner C,Wegmüller U,Strozzi T,et al.Interferometric point target analysis for deformation mapping[C]//Proceedings of 2003 IEEE International Geoscience and Remote Sensing Symposium.Toulouse,France:IEEE,2003,7:4362-4364.

[4] Yang H L,Peng J H.Monitoring urban subsidence with multi-master Radar interferometry based on coherent targets[J].Journal of the Indian Society of Remote Sensing,2015,43(3):529-538.

[5] Feng G C.Coseismic deformation and ionospheric variation associated with Wenchuan earthquake estimated from InSAR[D].Hong Kong:The Hong Kong Polytechnic University,2011.

[6] Lu Z,Mann D,Freymueller J T,et al.Synthetic aperture Radar interferometry of Okmok volcano,Alaska:Radar observations[J].Journal of Geophysical Research,2000,105(B5):10791-10806.

[7] Berardino P,Costantini M,Franceschetti G,et al.Use of differential SAR interferometry in monitoring and modelling large slope instability at Maratea (Basilicata,Italy)[J].Engineering Geology,2003,68(1/2):31-51.

[8] Ferretti A,Prati C,Rocca F.Permanent scatterers in SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(1):8-20.

[9] Ferretti A,Prati C,Rocca F.Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(5):2202-2212.

[10]Adam N,Kampes B M,Eineder M.Development of a scientific permanent scatterer system:Modifications for mixed ERS/Envisat time series[C]//Proceedings of the 2004 ENVISAT and ERS Symposium.Salzburg:ADS,2004.

[11]Hanssen R F.Radar Interferometry:Data Interpretation and Error Analysis[M].Netherlands:Springer,2001.

[12]陳剛,湯曉濤,余旭初,等.基于GCP的星載InSAR系統(tǒng)基線估計(jì)算法精度分析[J].測(cè)繪通報(bào),2009(6):16-18,21.

[13]楊紅磊,彭軍還,張丁軒,等.軌道誤差對(duì)InSAR數(shù)據(jù)處理的影響[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),2012,29(2):118-121,126.

[14]Zebker H A,Villasenor J.Decorrelation in interferometric Radar echoes[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(5):950-959.

[15]Cumming I G,Wong F H.合成孔徑雷達(dá)成像——算法與實(shí)現(xiàn)[M].洪文,胡東輝,譯.北京:電子技術(shù)出版社,2007.

[16]Souyris J C,Henry C,Adragna F.On the use of complex SAR image spectral analysis for target detection:Assessment of polarimetry[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(12):2725-2734.

[17]Berardino P,Fornaro G,Lanari R,et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2375-2383.

[18]安徽省地質(zhì)調(diào)查院.安徽省阜陽(yáng)市地面沉降調(diào)查報(bào)告[R].2004.

[19]安徽省地礦局第一水文地質(zhì)工程地質(zhì)隊(duì).安徽省阜陽(yáng)市水、工、環(huán)地質(zhì)綜合詳查報(bào)告(1/5萬(wàn))[R].1988—1990.

(責(zé)任編輯： 刁淑娟)

Study of land subsidence monitoring based on an improved coherent point target technology in Fuyang City

PENG Peng1, YANG Honglei2

(1.GeologicalSurveyInstitutionofAnhuiProvince,Hefei230001，China;2.ChinaUniversityofGeosciences(Beijing),Beijing100083，China)

This paper proposed an improved coherent point target InSAR technology. The value of the sub-look coherent was used to identify the coherent point targets of small dataset, and the trend filter was then applied to remove the large-scale residual phase. Using 10 Radarsat-2 Wide Model images over Fuyang City during 2012—2013, the high density coherent point targets were obtained and the reliability of the results were improved. The obtained land subsidence was in accordance with the previous measurement data and our field survey. This method provides the technical supports and application example for risk analysis and the regular monitoring of the land subsidence in Fuyang City.

coherent point target technology; land subsidence; filter; Radarsat-2

2016-03-10；

2016-05-31。

中國(guó)地質(zhì)調(diào)查局 “全國(guó)地表形變遙感地質(zhì)調(diào)查(安徽省)(編號(hào)： 1212011140037)”和安徽省國(guó)土資源廳“基于PS-InSAR的皖北地區(qū)長(zhǎng)時(shí)序地面沉降監(jiān)測(cè)方法研究(編號(hào): 2016-k-8)”項(xiàng)目共同資助。

彭鵬(1982—)，男，工程師，碩士研究生，主要從事遙感地質(zhì)調(diào)查工作。Email： 119170256@qq.com。

TP79

A

2095-8706(2016)04-0063-06

引用格式： 彭鵬,楊紅磊.基于改進(jìn)相干點(diǎn)目標(biāo)技術(shù)的阜陽(yáng)市地面沉降調(diào)查[J].中國(guó)地質(zhì)調(diào)查,2016,3(4): 63-68.