高分辨率TerraSAR-X時序差分干涉沉降監(jiān)測及精度驗證

于 冰， 譚青雪， 劉國祥， 劉福臻， 周志偉， 何智勇

(1.西南石油大學土木工程與測繪學院，成都 610500； 2.西南石油大學油氣空間信息工程研究所，成都 610500； 3.西南石油大學油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都 610500； 4.中國科學院精密測量科學與技術(shù)創(chuàng)新研究院，大地測量與地球動力學國家重點實驗室，武漢 430077; 5.西南交通大學地球科學與環(huán)境工程學院，成都 611756)

0 引言

地面沉降是指在一定的地面范圍內(nèi)發(fā)生的地面標高降低的現(xiàn)象，它通常是在自然和人為因素作用下引發(fā)的一種地質(zhì)災害[1-2]。地面沉降不僅會造成建筑物下沉、建筑基礎(chǔ)變形、地面塌陷，還會造成水準點高程失真等，更嚴重的是地面沉降可能會間接地引發(fā)或加劇其他自然災害的發(fā)生，如海平面上升、海水倒流、土地鹽堿化、地面裂縫等[3-6]。天津市是環(huán)渤海地區(qū)的經(jīng)濟發(fā)展中心，是第一批沿海開放的城市，位于華北平原北部，地理位置較為特殊，地面沉降災害尤為嚴重[6]，對天津市開展沉降監(jiān)測和分析具有重要的現(xiàn)實意義。

傳統(tǒng)的沉降監(jiān)測方法，如精密水準測量、全球定位系統(tǒng)(global positioning system,GPS)測量等，人工成本較大，作業(yè)時易受天氣影響，且監(jiān)測速度慢，監(jiān)測范圍較小[7-8]。合成孔徑雷達干涉(interferometric synthetic aperture Radar，InSAR)測量具有觀測范圍廣，空間覆蓋率和分辨率高、監(jiān)測效率高等優(yōu)勢[9]。尤其是時序InSAR，如永久散射體干涉(permanent scatterer InSAR，PSInSAR)、小基線集(small baseline subsets，SBAS)干涉和相干點目標分析(interferometric point target analysis, IPTA)等，在地面沉降、滑坡形變監(jiān)測等方面得到了較多的應用[10-14]。

IPTA是在PSInSAR技術(shù)基礎(chǔ)上加以改進而來，相對于PSInSAR技術(shù)來說，IPTA技術(shù)提高了相干點的空間采樣率和解算效率[10]。時序差分干涉測量(differential InSAR,DInSAR)雖有較多應用，但其精度驗證仍是目前研究中的重要組成部分，現(xiàn)有文獻大多會采用GPS或水準數(shù)據(jù)評估其精度。隨著高時空分辨率SAR數(shù)據(jù)的投入使用，往往存在檢核數(shù)據(jù)的時間采樣率不足而無法對時序DInSAR所獲取的時序沉降進行驗證。比如，實際中GPS和水準觀測周期比SAR衛(wèi)星重訪周期長，尤其是很多區(qū)域只能獲取極少甚至只有一期觀測數(shù)據(jù)，如何對時序沉降進行驗證仍是一個需要解決的問題。

本文以天津市西青區(qū)及其周邊區(qū)域為研究區(qū)域，以高分辨率的TerraSAR-X SAR影像數(shù)據(jù)為數(shù)據(jù)源，采用IPTA技術(shù)進行沉降信息的提取與分析，并提出一種最小二乘擬合的沉降時間序列驗證方法，對監(jiān)測結(jié)果進行精度驗證，進一步檢驗高分辨率時序DInSAR應用于城市沉降監(jiān)測的精度，并為缺少時序地面測量數(shù)據(jù)的情況下驗證時序DInSAR結(jié)果提供一種有效的技術(shù)思路。

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

1.1 研究區(qū)概況

天津市處于華北平原北部，向東鄰近渤海，向北依靠燕山[15]。土地總面積約為11 966 km2，地質(zhì)構(gòu)造復雜。早期，天津市的地面沉降主要發(fā)生在城市地區(qū)，但是隨著天津市經(jīng)濟與城市建設(shè)的外圍發(fā)展，并加上主城區(qū)地下水開采控制，地面沉降逐漸從城市擴展到郊區(qū)[16]。研究區(qū)域主要包括西青區(qū)及北辰區(qū)部分區(qū)域，如圖1所示，紅色框為影像范圍，底圖為數(shù)字高程模型(digital elevation model,DEM)。

圖1 研究區(qū)域示意圖Fig.1 Study area

1.2 數(shù)據(jù)源

本文選取了2009年4月7日—2010年12月14日的34幅高分辨率TerraSAR-X降軌影像數(shù)據(jù)作為時序DInSAR處理的數(shù)據(jù)源，所有影像均為單視復數(shù)格式，波長為3.1 cm，圖像成像空間分辨率約為3 m。在考慮影像獲取時間(時間居中，時間最好在春冬季，且天氣狀況較好為最適宜)和基線長度(基線差較小為最好)的條件下，選取了2009年11月13日的影像為主影像，其他33幅影像為輔影像。確定一景影像為主影像后，要先對影像做預處理： 將其他33景輔影像配準重采樣到主影像的幾何空間。影像獲取日期及輔影像相對于主影像的時空基線參數(shù)如表1所示。

表1 干涉對及時空基線Tab.1 Interferometric pairs and spatiotemporal baselines

2 IPTA基本原理及數(shù)據(jù)處理流程

IPTA的核心思路是從配準好的SAR影像集中選取具有穩(wěn)定散射特性的相干點目標，提取其相位信息并以點為單位進行DInSAR處理，以線性形變模型和高程殘差為基礎(chǔ)建立二維線性回歸模型，并迭代計算線性形變速率和高程殘差，然后從扣除線性形變和高程殘差的殘余相位中分離非線性形變、大氣延遲和軌道誤差等，最后將線性形變和非線性形變相疊加獲得時間序列[5,17]。

在本文中，基于所獲取的34幅TerraSAR-X SAR影像，采用IPTA提供的時序振幅統(tǒng)計穩(wěn)定性和光譜相關(guān)性2種方法選取相干點目標，并對二者進行合并，以此增加點目標的密度。然后，從34幅復數(shù)影像中提取相干點目標的原始相位進行干涉處理，本文采取公共主影像干涉組合模式，共獲取33幅干涉相位圖。此時，相干點上的干涉相位由參考橢球相位、地形相位、形變相位、大氣延遲相位和失相干噪聲等組成。接下來，基于所獲取的外部DEM計算參考橢球相位和地形相位并從干涉相位中去除，得到33幅差分干涉圖。此時，干涉圖上第i個相干點的差分干涉相位可表達為：

φdiff,i=φdef,i+φtop_res,i+φatm,i+φnoi,i+φorb_res,i

(1)

式中：φdiff,i為相干點上的差分干涉相位；φdef,i為形變相位；φtop_res,i和φorb_res,i分別為殘余地形及軌道誤差相位；φatm,i和φnoi,i分別為大氣延遲和噪聲相位。

在得到相干點上的差分干涉相位后，IPTA基于鄰域差分的思路建立相鄰相干點之間的相位增量(差異)模型，該模型是關(guān)于相鄰2點間線性形變速率增量和地形殘差增量的線性模型，并包含其他噪聲或者誤差項，具體可表達為：

(2)

式中: Δφdiff,i-j為2相鄰相干點間的二次差分相位； Δν為2相鄰相干點間的線性形變速率增量；Βt為差分干涉圖的時間基線；λ,R和θ分別為雷達波長、成像斜距和入射角； Δδh為2相干點間的高程改正增量；Bs為差分干涉圖的空間基線；φorb_res,i-j,φatm,i-j和φnoi,i-j分別為相應的軌道誤差增量、大氣延遲增量和噪聲增量。

在上面模型的基礎(chǔ)上，IPTA方法數(shù)據(jù)處理的流程如圖2所示。IPTA以相干點上所有的差分干涉相位為觀測值，采用線性回歸最優(yōu)化求解線性形變速率增量和高程增量，并選定參考點求解所有相干點上相對于參考點的線性形變速率和高程改正。本文中IPTA解算的參考點經(jīng)緯度為(E116.976 966 9°，N39.093 383 8°)。除此之外，IPTA提供迭代策略逐步精化這2個線性分量的解算結(jié)果。然后，從差分干涉相位中扣除線性形變和高程改正值，得到殘余相位，主要由非線性形變、大氣延遲、軌道誤差和失相干噪聲等分量組成。根據(jù)這些相位分量的時空頻率特性不同，IPTA通過空間域和時間域濾波從殘余相位中分離出非線性形變。將線性形變和非線性形變疊加即可得到形變時間序列，最后將形變時間序列轉(zhuǎn)換到垂直向得到沉降時間序列。

圖2 IPTA數(shù)據(jù)處理流程Fig.2 Flow chart of IPTA data processing

3 IPTA數(shù)據(jù)處理結(jié)果及精度驗證

3.1 沉降速率精度驗證

根據(jù)第2部分的數(shù)據(jù)處理方法和步驟，本文提取了1 231 705個相干點，獲取了研究區(qū)域內(nèi)的年沉降速率，如圖3所示。

圖3 相干點沉降速率及水準點分布Fig.3 Distribution of coherent point subsidence rate and leveling points

為了驗證本文所提取沉降結(jié)果的可靠性與精度，根據(jù)水準實測沉降數(shù)據(jù)對所得結(jié)果進行檢驗(為保證后續(xù)結(jié)果分析的可靠性，先對監(jiān)測結(jié)果進行驗證)?，F(xiàn)有研究也多是采用這種驗證形式對雷達DInSAR結(jié)果進行檢驗[2,5,8,10,12]。在水準沉降監(jiān)測中，采用二等精密水準測量，獲取不同時期的高程數(shù)據(jù)，通過高程作差得到相應的沉降量。高程測量過程中每km高差中誤差為1 mm，滿足沉降監(jiān)測規(guī)范要求，因此可以作為驗證IPTA結(jié)果精度的基準數(shù)據(jù)。在此需要說明的是，由于時序DInSAR處理方法得到的沉降速率為相對沉降速率，而水準數(shù)據(jù)的沉降速率為絕對值，所以選用其中的一個水準點為校正點，來校正其他相干點，利用校正后的結(jié)果進行對比分析。水準點的空間分布如圖3中黑色三角形所示，共計13個(1個基準點和12個檢核點)。首先，從IPTA沉降結(jié)果中提取出與這13個水準點相近的相干點(如圖3中紅色圓點所示)，并記錄沉降結(jié)果，用基準點對IPTA結(jié)果進行校正后，采用12個檢核點進行驗證。

在對二者進行對比分析時，采用差異均值和均方根誤差[1,3](root mean square error,RMSE)2個指標，其公式分別為：

(3)

(4)

式中：k=12；Μv為水準點與相干點沉降速率差異的均值；vKi表示第i個水準點的沉降速率值；vLi表示第i個相干點的沉降速率值；N為點的數(shù)目；sv表示水準點與相干點沉降速率的RMSE。

本文選取12號點為校正點，其余12個水準點作為檢核點。如圖4所示為13個水準點和相干點沉降速率值的差異對比(12號點為校正點，差異為0)。圖4中紅色折線表示相干點的沉降速率，黑色折線表示水準點的沉降速率，藍色折線表示水準點和相干點沉降速率的差異。從圖4中可看出，這13個點的沉降速率值大致在-45 mm/a和-20 mm/a之間，12個檢核點中，有8個點差值在3 mm/a內(nèi)，差異最大值為-6.62 mm/a，最小值為0.06 mm/a，從圖4中藍色折線的走勢也可看出相干點的沉降速率與水準點處的沉降速率偏差較小[2]。根據(jù)式(3)—(4)計算得到水準點和相干點沉降速率差異的均值和RMSE分別為： 1.02 mm/a和±3.15 mm/a。上述對比分析表明結(jié)合高分辨率TerraSAR-X數(shù)據(jù)的IPTA技術(shù)所提取的沉降速率與精密水準沉降監(jiān)測結(jié)果具有較好的符合度。

圖4 水準點與相干點沉降結(jié)果對比Fig.4 Comparison of subsidence between leveling points and cohevent points

3.2 基于最小二乘擬合的沉降時間序列精度驗證

前已述及，在缺少時序水準數(shù)據(jù)(尤其是僅能獲取1期水準數(shù)據(jù))時，如何對時序DInSAR所獲取的時序沉降進行驗證是需要解決的問題。在本文中，為了進一步驗證時序DInSAR處理方法的精度及可靠性，使用基于最小二乘的方法，根據(jù)沉降時間序列和觀測時間段來擬合沉降速率，并分析擬合沉降速率與水準沉降速率的差異。對于每個相干點，可建立K個觀測方程，即

(5)

式中：ΔtK表示第K景影像到第一景影像的時間間隔，d;v表示相干點上的沉降速率，每一個相干點對應于一個v;dK表示相干點在第K景影像所對應時間的沉降量;vdK為對應沉降量的改正數(shù)。將上述觀測方程轉(zhuǎn)換成矩陣形式為：

(6)

據(jù)此可解得每個相干點的沉降速率為：

v=(BΤPB)-1BTL

(7)

式中:B=[Δt1,Δt2,…,ΔtK]Τ;L=[d1,d2,…,dK]Τ;P為單位權(quán)陣。

在得到相干點上的擬合沉降速率后，采用和3.1節(jié)中相同的方案進行驗證。相干點上擬合沉降速率和水準沉降速率的折線圖如圖5所示。

圖5 擬合沉降速率驗證結(jié)果Fig.5 Result of fitted subsidence rate verification

圖5中藍色折線表示水準沉降速率，綠色折線表示擬合的沉降速率，橙色折線表示水準沉降速率與擬合沉降速率差異值，差異絕對值的最大值為7.06 mm/a。根據(jù)式(3)—(4)計算得到水準沉降速率與擬合沉降速率差異的均值和RMSE分別為-3.4 mm/a和±3.25 mm/a，上述對比分析表明通過IPTA時序DInSAR處理方法所獲取的時序沉降結(jié)果擬合所得的沉降速率與精密水準沉降監(jiān)測結(jié)果具有較好的一致性。

4 研究區(qū)域沉降及成因分析

4.1 研究區(qū)域整體沉降分析

基于上述IPTA方法對2009年4月—2010年10月天津市西青區(qū)34景高分辨率TerraSAR-X SAR數(shù)據(jù)進行分析處理，獲得了該段時間的年均沉降速率，結(jié)果如圖6所示。在數(shù)據(jù)處理過程中共獲取了1 231 705個相干點，這些相干點大多分布于建筑物、構(gòu)筑物和道路等具有較高散射特性的地物上。從沉降速率結(jié)果中可看出，在監(jiān)測期間研究區(qū)域地面沉降總體呈現(xiàn)不均勻性，其中與西青區(qū)鄰近的北辰區(qū)沉降速率較大。

圖6 沉降速率結(jié)果Fig.6 Result of subsidence rate

具體來說，由圖6可看出研究區(qū)域不同地方的沉降情況，不同的顏色等級代表不同的沉降速率大小，從顏色分布情況來看，可將沉降速率分布大致分為3個級別。如圖中暗紅色、紅色、橙色、黃色部分位于北辰區(qū)附近，沉降速率較大，在-128.41～-74.26 mm/a范圍內(nèi)，說明該地區(qū)沉降較為嚴重； 綠色、青色部分主要位于西青區(qū)附近，其沉降速率在-74.25～-33.65 mm/a范圍內(nèi)； 藍色和深藏青色部分覆蓋了北辰區(qū)部分地區(qū)、紅橋區(qū)和南開區(qū)3個區(qū)域，沉降速率在-33.64～2.45 mm/a之間，該沉降區(qū)間覆蓋的范圍較大。由沉降結(jié)果及沉降分布情況可知，該研究區(qū)主要的沉降位于北辰區(qū)和西青區(qū)，結(jié)合相關(guān)資料可知該區(qū)域范圍內(nèi)存在較多的居民地、工業(yè)園區(qū)和高大建筑物聚集地等，工業(yè)生產(chǎn)用水、居民生活用水以及建筑物外部荷載等都會造成地面沉降。

為了對形變結(jié)果的空間分布規(guī)律進行進一步分析，分別在西青區(qū)和北辰區(qū)各提取了1個局部剖面AB和CD。剖面AB、剖面CD的沉降剖面圖，如圖7所示，圖7中沉降速率是進行克里金插值后的值。如圖7中所示，紅色實線表示AB的沉降剖面，可看出在剖面上距離A點不同位置處的沉降速率，在距離A點約1.3 km處沉降速率最大。藍色實線表示CD的沉降剖面，在距C點約2.1 km處沉降速率最大，是該區(qū)段的沉降漏斗中心。

圖7 AB及CD的沉降剖面Fig.7 Subsidence profile of AB and CD

提取圖6中點M，N，O，P，Q，R，S，T的沉降時間序列如圖8所示，即根據(jù)累積沉降量與時間間隔的關(guān)系繪制的形變時間序列圖。圖8中所示，點M和點N處的累計沉降量最大，點M處的累積沉降量最大達到-170 mm，點O處最大沉降量達到-165 mm。點M和O所在區(qū)域為居民地和工業(yè)園區(qū)，年平均沉降速率在-100 mm/a左右。點N,P,Q,R,S,T沉降量在120 mm左右，所在區(qū)域為居民地、工業(yè)園區(qū)以及交通要道，年平均速率在-80 mm/a以上。

圖8 沉降時間序列Fig.8 Subsidence history

4.2 典型沉降區(qū)及交通要道沉降分析

為了更好地對沉降區(qū)域的沉降情況及成因進行探討和分析，本文從整個研究區(qū)域中提取了部分塊狀區(qū)域和線狀地物的沉降速率，并獲取了對應地物的可見光影像。圖9(a)為北辰區(qū)的沉降速率，圖9(b)為西青區(qū)的沉降速率。

(a) 北辰區(qū)沉降速率 (b) 西青區(qū)沉降速率

(c) 工業(yè)園區(qū)c (d) 居民地d (e) 居民地e (f) 工業(yè)園區(qū)f圖9 2個地區(qū)的沉降速率及相應的用地類型Fig.9 Subsidence rate and land use category analysis of the two selected areas

在圖9(a)中選取了c、d、e這3個沉降速率較大的區(qū)域，在圖9(b)中選取了沉降速率較大的f區(qū)域，圖9(c),(d),(e),(f)為其相應地區(qū)的用地類型。圖9(d)和(e)所在土地用地為居民地，(c)和(f)用地類型為工業(yè)園區(qū)。結(jié)合圖6和圖9，可知圖6中紅色和黃色區(qū)域(圖9中(c),(d),(e),(f))在實地上是工業(yè)園區(qū)、廠房以及成片的居民地等，圖9(a)中c和d區(qū)域的沉降速率最大達到了-128.41 mm/a，說明該區(qū)域每年的沉降量最多。c和d分別對應于圖9(c)和(d)，用地類型主要為工業(yè)區(qū)和居民地，該區(qū)域沉降如此明顯，推測可能與工業(yè)生產(chǎn)和居民生活對地下水資源需求較大進而導致過度開采地下水相關(guān)。

研究區(qū)域內(nèi)3條高速公路局部路段的沉降速率如圖10所示，分別為津保高速、京臺高速和榮烏高速公路。從圖10中可知，津保高速公路的沉降速率最大，沉降速率大致在-92.3～-60.73 mm/a之間，高速公路附近有聚集的居民地、工廠和農(nóng)田等，地下水資源需求較大，大量的地下水開采導致了嚴重的地面沉降。京臺高速公路和榮烏高速公路的沉降速率比較接近，沉降速率大致在-51.69～-33.65 mm/a之間，高速公路附近也有農(nóng)田和居民地，工廠較少，主要的地下水需求為灌溉和生活用水，開采量相對較小，所以沉降速率相對較小。

(a) 津保高速 (b) 京臺高速 (c) 榮烏高速

圖10 3條高速公路沉降速率Fig.10 Subsidence rate of three highways

5 結(jié)論

本文以天津市西青區(qū)及其周邊地區(qū)為研究區(qū)域，以覆蓋研究區(qū)的34幅X波段TerraSAR-X SAR影像為數(shù)據(jù)源，使用基于IPTA的時序差分干涉處理方法，提取了研究區(qū)域的地面沉降信息，并基于精密水準數(shù)據(jù)進行驗證，給出一種基于最小二乘沉降速率擬合的沉降時間序列驗證方法。最后基于驗證后的沉降結(jié)果分析了研究區(qū)的區(qū)域沉降、主要道路沉降及其成因。研究表明：

1)經(jīng)與水準數(shù)據(jù)對比，IPTA解算的沉降速率、基于最小二乘擬合的沉降速率與水準點的沉降速率均較為吻合，進一步驗證了IPTA時序差分干涉處理方法的精度和可靠性。本文所提出的基于最小二乘擬合的時序差分干涉沉降時間序列驗證方法可在缺少時序地面測量數(shù)據(jù)的情況下對監(jiān)測所得時序結(jié)果進行驗證，在一定程度上可以緩解缺少地面驗證數(shù)據(jù)的問題，并完善時序差分干涉沉降監(jiān)測結(jié)果驗證的技術(shù)方案。

2)通過對整個研究區(qū)域沉降結(jié)果分析可知，在2009—2010年間，研究區(qū)域總體呈現(xiàn)出顯著的不均勻性，區(qū)域內(nèi)沉降最嚴重的區(qū)域在北辰區(qū)一帶，最大年平均沉降速率達到了-128.41 mm/a，這與現(xiàn)有研究較為吻合。結(jié)合該研究區(qū)的人口分布特點、經(jīng)濟發(fā)展情況和地物分布情況，并對研究區(qū)域進行局部沉降分析可知，造成區(qū)域沉降不均勻的原因是不同地物類型對地下水的需求量不同，進而造成地下水開采程度有所差異，引發(fā)不均勻沉降。

綜上所述，基于高分辨率TerraSAR-X數(shù)據(jù)的IPTA時序差分干涉處理方法在城市沉降速率及時間序列監(jiān)測中具有較高的精度和可靠性，可以獲取詳細的地表沉降分布情況。同時，高分辨率SAR影像在監(jiān)測道路等線狀地物沉降方面表現(xiàn)出較好的效果。本文所提出的基于最小二乘擬合的時序沉降驗證方法可有效緩解地面驗證數(shù)據(jù)不足的問題。由于地面測量數(shù)據(jù)較難獲取，本文未能采用更多的水準監(jiān)測點進行結(jié)果驗證，為了進一步完善時序沉降驗證方案，后續(xù)研究將進一步對所提出的時序沉降驗證方法進行更多的論證。此外，采用最新的SAR數(shù)據(jù)對該區(qū)域進行沉降監(jiān)測并與歷史沉降進行對比分析，探討區(qū)域沉降的時空演變規(guī)律也是非常值得研究的問題。