熊 威，孫志杰，張必昌

（1.山西省交通科技研發(fā)有限公司，山西 太原 030032）

合成孔徑雷達(dá)干涉測量（InSAR）技術(shù)作為一項新型的衛(wèi)星對地觀測手段，具有觀測頻率高、精度高、耗資少等優(yōu)勢，由此發(fā)展而來的PS-InSAR（永久散射體雷達(dá)干涉測量）和SBAS-InSAR（小基線集雷達(dá)干涉測量）技術(shù)近十幾年來在各大城市地面沉降中發(fā)揮了重要的作用。針對天津地區(qū)地面沉降監(jiān)測與研究的工作也開展了很多[1-3]，但是缺乏最近2 a 的研究成果，并且使用的數(shù)據(jù)均為單一升軌或降軌。本文通過歐空局提供的2017-2019 年間的升降軌Sentinel-1A/B 衛(wèi)星影像，利用SBAS-InSAR 時序處理方法獲得天津地區(qū)地面沉降速率，并對升降軌結(jié)果進(jìn)行相互檢校和融合，利用融合后的結(jié)果分析天津地區(qū)近年來的面狀沉降特征和京滬高鐵線狀沉降特征。

1 SBAS-InSAR數(shù)據(jù)處理方法

傳統(tǒng)D-InSAR（差分雷達(dá)干涉測量）技術(shù)得到的干涉相位組成為:

式中，φdef為形變相位；φdem為DEM 誤差產(chǎn)生的干涉相位；φatm為大氣延遲產(chǎn)生的干涉相位；φord為軌道誤差相位；φnoi為系統(tǒng)噪聲相位。其面臨的最主要問題為:許多地區(qū)的 SAR 圖像在進(jìn)行干涉處理時由于受到時間空間的去相干而失敗，大氣延遲更是大大的降低了數(shù)據(jù)處理的精度。

為了解決上述問題，最早由 Berardino[4]提出SBASInSAR 處理方法，基本原理是通過對獲取到的大量SAR 數(shù)據(jù)選取短的空間基線像對進(jìn)行適當(dāng)?shù)慕M合，然后計算出這些組合的多視干涉圖，最終生成平均形變速度圖或得地表形變時間序列。SBAS 方法將所有的SAR 數(shù)據(jù)組成若干個集合，集合內(nèi)的基線距較小，集合間的 SAR 圖像基線大。對子基線集之間的組合利用矩陣的奇異值分解（SVD）方法求出未知參數(shù)在最小范數(shù)意義上的最小二乘解，最后考慮到大氣相位在空間上的高相關(guān)性及時間上的低相關(guān)性，利用時間和空間濾波估計出非線性形變和大氣相位成分，處理流程如圖1 所示。

圖1 SBAS-InSAR 處理流程

首先根據(jù)同一地區(qū)不同時刻（t1、t2、t3、…、tn）的N幅SAR 影像組成M個干涉像對，其中N/2 ≤M≤N×(N-1)/2，由公式（1）可知第k個差分干涉像對上的第x個像元的解纏相位組成部分有地形形變相位、地形相位殘差、不同時刻ti、tj的大氣相位差值、軌道相位殘差和系統(tǒng)誤差相位。在數(shù)據(jù)預(yù)處理階段初步除去差分干涉相位中的大氣和軌道相位，相位殘差均包含在隨機(jī)誤差項中，使得最終的差分相位成分簡化為地形形變、地形殘差和隨機(jī)誤差。其中地形形變有:

式中，A為M×N階矩陣；x為N個時間序列上的形變量與地形殘差組成的矩陣；由公式（2）可知矩陣A的每一行主影像所在歷元的系數(shù)為-1，輔影像所在歷元的系數(shù)為1，其他為0，則矩陣A可表示為公式（4）的形式。

由SBAS 方法的技術(shù)特點可知非零數(shù)據(jù)分散在多個子集中，使得方程（3）具有較強(qiáng)的奇異性，會產(chǎn)生無數(shù)多組解，此時可利用SVD 方法對方程（3）進(jìn)行求解，進(jìn)而求得到歷元間形變量的最小范數(shù)解，根據(jù)方程（5）可進(jìn)一步得到地形殘差。利用空間上的低通濾波和時間上高通濾波對線性模型的殘余相位進(jìn)行處理就能得到大氣相位和非線性形變相位，疊加之前的形變估計值即可得到所有像元在不同時間序列上的形變結(jié)果。

通過上述原理可知，SBAS-InSAR 處理方法相對于傳統(tǒng)的D-InSAR 方法具有以下優(yōu)點:①增加時空基線長短限制，SAR 數(shù)據(jù)集內(nèi)影像相干性較好，失相關(guān)現(xiàn)象和地形誤差影響降低，同時數(shù)據(jù)之間相互組合使得數(shù)據(jù)的約束性更強(qiáng)；②由于數(shù)據(jù)之間的相干性較好，增強(qiáng)了對大氣、 DEM 等誤差的估計分析能力，測量精度高，可測得mm 級的視線向運動速度，結(jié)果更為可靠。本文利用GMTSAR 軟件中SBAS 模塊對天津地區(qū)的SAR 數(shù)據(jù)進(jìn)行處理，該軟件具有開源、自動處理、易于繪圖等優(yōu)點，應(yīng)用效果較好[5]。

2 天津市地面沉降監(jiān)測

天津市地處華北平原東北部，東臨渤海，北接北京，西部和南部被河北環(huán)繞，是北方的經(jīng)濟(jì)中心和特大城市。天津市面積廣袤，下轄16 個區(qū)，人口眾多，水資源相對匱乏，除市六區(qū)外，其他各區(qū)都存在規(guī)模以上的工業(yè)片區(qū)、農(nóng)業(yè)片區(qū)和經(jīng)濟(jì)開發(fā)區(qū)，在建設(shè)高樓大廈的同時也需要消耗大量的地下水資源，地下水超采已成為天津地區(qū)地面沉降產(chǎn)生的最主要原因[6]。

2.1 數(shù)據(jù)獲取

本文采用的升軌遙感影像為歐空局提供的Sentinel-1A 干涉寬幅衛(wèi)星數(shù)據(jù)，降軌遙感影像為歐空局提供的Sentinel-1B 數(shù)據(jù)，兩者均為側(cè)視成像，分別從東西向兩個不同的角度獲取目標(biāo)的形變信息，空間分辨率均為5 m×20 m（距離向×方位向）。數(shù)據(jù)時間跨度為2017 年1 月至2019 年1 月，升軌影像數(shù)量為25 景，降軌影像數(shù)量為24 景（2018 年9 月缺失），時間間隔為24 d 或36 d，以此保證每月一景，影像信息見表1。為提高計算效率，對整幅影像進(jìn)行裁剪后計算。計算中輔助DEM 為美國太空總署提供的SRTM3數(shù)據(jù)，其地面分辨率為3 弧秒，高程精度約為±16 m，用于移去干涉圖中的地形相位。精密水準(zhǔn)數(shù)據(jù)采用均勻分布于天津地區(qū)的327 個一等水準(zhǔn)點的測量結(jié)果，水準(zhǔn)點位置及升降軌影像范圍如圖2 所示。

圖2 研究區(qū)域及影像覆蓋范圍

表1 影像信息表

2.2 升降軌InSAR 監(jiān)測結(jié)果

設(shè)置空間基線和時間基線閾值分別為200 m、500 d，解纏閾值為0.2，解纏方法為最小費用流（MCF），距離向方位向多視比為8∶2，像元分辨率為30 m，按照圖1 的流程分別對升降軌數(shù)據(jù)進(jìn)行處理。由于時序InSAR 得到的形變結(jié)果為衛(wèi)星視線向，考慮到華北地區(qū)的水平位移很小[7]，相比垂直位移基本可以忽略，故可利用公式ΔH=ΔLOS/cosθ（θ為雷達(dá)入射角，ΔLOS為視線向形變量，ΔH為垂向形變量）將衛(wèi)星視線向形變量轉(zhuǎn)化到垂直向上，視作地面沉降結(jié)果，得到天津地區(qū)2017 年1 月至2019 年1 間的年平均沉降速率結(jié)果如圖3 所示。

從圖3 可以看出，天津市六區(qū)地面整體呈輕微抬升趨勢，年上升量小于30 mm，沉降嚴(yán)重的區(qū)域（＞30 mm/yr）基本位于西青區(qū)南部、津南區(qū)、北辰區(qū)西部、武清區(qū)南部、靜海區(qū)、漢沽東部和寧河區(qū)東南部，最大沉降量級達(dá)到156 mm/yr，位于武清區(qū)王慶坨鎮(zhèn)?？傮w來看，天津地區(qū)沉降趨勢明顯，多個局部沉降中心分散在市區(qū)周邊，沉降區(qū)域有形成集中連片之勢。

圖3 升降軌沉降速率圖

2.3 與精密水準(zhǔn)監(jiān)測結(jié)果比較分析

利用天津地區(qū)2017-2018 年間的一等水準(zhǔn)監(jiān)測結(jié)果對升降軌InSAR 結(jié)果進(jìn)行驗證，由于水準(zhǔn)的施測時間在每年的9 月份，本文提取相應(yīng)時段的InSAR 結(jié)果進(jìn)行驗證?？紤]到水準(zhǔn)點與永久散射體點（PS）位置不重合及系統(tǒng)誤差，取水準(zhǔn)點周邊0.1 km 范圍內(nèi)PS 點的平均值作為參照結(jié)果，當(dāng)單個水準(zhǔn)點周邊的PS 點少于3 個時，放棄該水準(zhǔn)點作為對比目標(biāo)。提取后的升降軌及水準(zhǔn)結(jié)果如圖4a、4b 所示，對提取后的升降軌結(jié)果進(jìn)行融合，考慮圖像分辨率及圖像截取范圍，對于升降軌影像之間像元間隔小于15 m 的像素點當(dāng)成同名點，取形變均值，融合后的結(jié)果如圖4c 所示，圖像對比統(tǒng)計分析結(jié)果如表2 所示。

從圖4 和表2 可以看出升降軌影像的年沉降速率與一等水準(zhǔn)結(jié)果的差值基本在±15 mm 以內(nèi)，主要原因可能有:①InSAR 數(shù)據(jù)結(jié)果的時間范圍無法與水準(zhǔn)測量時間完全一致；②SBAS 方法解算時，解纏控制點的選擇包含一定的誤差。InSAR 數(shù)據(jù)與一等水準(zhǔn)數(shù)據(jù)的整體誤差在8 mm 左右，并且呈現(xiàn)出一定的區(qū)域集中性，東、北部多為正值，西、南部多為負(fù)值，可能與水準(zhǔn)測量的先后順序相關(guān)。由于升降軌影像分別從西、東2 個方向以約36°的入射角對地面進(jìn)行干涉測量，能夠全面反映地物的形變特征，故升降軌融合影像不僅能顯著增加像元密度，有助于分析線狀地物的形變特征，還能在一定程度上提高監(jiān)測結(jié)果精度。

圖4 升降軌InSAR 結(jié)果精度驗證

表2 圖像對比結(jié)果統(tǒng)計表

2.4 地面沉降結(jié)果分析

天津地區(qū)2017-2018 年升降軌融合沉降結(jié)果如圖5 所示。通過圖5 可知，天津地區(qū)沉降中心位于武清區(qū)王慶坨鎮(zhèn)，最大年沉降速率達(dá)到156 mm，與Zhang[1]等測得2012 年1 月至2014 年7 月間的最大年沉降速率基本相同，明顯小于李廣宇[2]等測得2015 年6 月至2016 年8 月間的最大年沉降速率224 mm/yr，略小于Zhang[8]等測得2016 年1 月至2017 年6 月間的最大年沉降速率163 mm/yr，其原因可能為2015-2016 年間的王慶坨鎮(zhèn)共享單車大量生產(chǎn)造成用水量大增，使得沉降加劇。隨著2016 年下半年開始，共享單車的生產(chǎn)量萎縮以及一批機(jī)井和地?zé)峋姆馓罘獯?，沉降趨勢得到緩和。漢沽區(qū)最大沉降速率達(dá)到98 mm/yr，其主要原因歸咎于水產(chǎn)養(yǎng)殖政策的支持下地下水開采不斷增加，靜海區(qū)作為天津市農(nóng)業(yè)發(fā)展程度較高的地區(qū)，對地下水的需求量一直很大，自上世紀(jì)80 年代初就一直是幾大沉降中心之一[6]，近年來也有愈發(fā)嚴(yán)重的趨勢，目前最大沉降速率已達(dá)到62 mm/yr。

圖5 天津地面沉降升降軌InSAR 融合結(jié)果

值得關(guān)注的是，京滬高鐵在天津境內(nèi)貫穿武清、靜海、北辰等區(qū)域，在此利用融合后的InSAR 監(jiān)測結(jié)果對京滬線周邊的沉降情況做進(jìn)一步的分析研究。提取天津境內(nèi)京滬高鐵沿線50 m 范圍內(nèi)的PS 點，繪制形變速率曲線如圖6 所示，數(shù)據(jù)的連續(xù)性較好，能完整反映線路的沉降特征。結(jié)果顯示有三處明顯的沉降中心A1、A2、A3，其中A1、A2 均位于靜海區(qū)內(nèi)，A1 處沉降速率達(dá)到50 mm/yr，其周邊沉降速率大于20 mm/yr 的南北向?qū)挾菳1-B2 達(dá)到14 km；A2 處沉降速率達(dá)到39 mm/yr，其沉降跨度較A1 窄；A3 位于武清區(qū)中部，其沉降速率達(dá)到32 mm/yr，周邊局部小沉降中心較多。由于位置靠近，A1、A2 的形變趨勢相似，反映出先緩后急的非線性變化，A3 的線性變化特征則比較明顯。

圖6 京滬高鐵線沉降監(jiān)測結(jié)果

3 結(jié) 語

本文利用2017 年1 月至2019 年1 月間的升降軌Sentinel-1A/B 數(shù)據(jù)對天津地區(qū)地面沉降進(jìn)行監(jiān)測，并對SBAS-InSAR 方法獲得的升降軌沉降結(jié)果進(jìn)行融合，利用市域范圍內(nèi)的一等水準(zhǔn)監(jiān)測結(jié)果對InSAR 結(jié)果進(jìn)行驗證，驗證結(jié)果表明升降軌融合不僅能顯著增加局部PS 點位密度，還能一定程度上提高精度。利用融合后的結(jié)果分析了天津地區(qū)和京滬高鐵線近年來的沉降特征，王慶坨鎮(zhèn)、漢沽東部和靜海區(qū)的沉降趨勢比較明顯，受制于地理位置和環(huán)境，天津境內(nèi)的京滬高鐵線也存在三處沉降速率大于20 mm/yr 的沉降中心，其中兩處位于靜海區(qū)內(nèi)，值得警惕和進(jìn)一步關(guān)注。研究結(jié)果表明升降軌融合結(jié)果對于分析線狀地物沉降特征有一定的優(yōu)勢，能較好的保持?jǐn)?shù)據(jù)的連續(xù)性和完整性。