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

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

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

0 引言

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

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

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

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

1 研究區概況及其數據源

1.1 研究區概況

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

圖1 研究區域示意圖Fig.1 Study area

1.2 數據源

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

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

2 IPTA基本原理及數據處理流程

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

在本文中，基于所獲取的34幅TerraSAR-X SAR影像，采用IPTA提供的時序振幅統計穩定性和光譜相關性2種方法選取相干點目標，并對二者進行合并，以此增加點目標的密度。然后，從34幅復數影像中提取相干點目標的原始相位進行干涉處理，本文采取公共主影像干涉組合模式，共獲取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基于鄰域差分的思路建立相鄰相干點之間的相位增量(差異)模型，該模型是關于相鄰2點間線性形變速率增量和地形殘差增量的線性模型，并包含其他噪聲或者誤差項，具體可表達為：

(2)

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

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

圖2 IPTA數據處理流程Fig.2 Flow chart of IPTA data processing

3 IPTA數據處理結果及精度驗證

3.1 沉降速率精度驗證

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

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

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

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

(3)

(4)

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

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

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

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

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

(5)

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

(6)

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

v=(BΤPB)-1BTL

(7)

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

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

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

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

4 研究區域沉降及成因分析

4.1 研究區域整體沉降分析

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

圖6 沉降速率結果Fig.6 Result of subsidence rate

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

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

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

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

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

4.2 典型沉降區及交通要道沉降分析

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

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

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

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

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

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

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

5 結論

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

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

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

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