Sentinel-1A京津冀平原區2016—2018年地面沉降InSAR監測

史 珉， 宮輝力， 陳蓓蓓， 高明亮， 張舜康

(1.首都師范大學地面沉降機理與防控教育部重點實驗室，北京 100048； 2.首都師范大學水資源安全北京實驗室，北京 100048； 3.首都師范大學北京市城市環境過程與數字模擬國家重點實驗室培育基地，北京 100048； 4.首都師范大學三維信息獲取與應用教育部重點實驗室，北京 100048)

0 引言

地面沉降是一種地表高程下降的環境地質現象，其發生發展增加了城市內澇、海水倒灌等風險，并可以誘發地面塌陷、地裂縫等系列環境災害，形成災害鏈。區域地面沉降不僅對城市基礎設施、高速鐵路等重大工程產生潛在威脅，還影響著城市地下空間資源利用效率，成為制約區域可持續發展的重大問題。

京津冀地區所處的華北平原是我國沉降速率最大、影響面積最大的地區[1]。自1923年首次在天津發現地面沉降以來[2]，華北平原地面沉降按期發展過程可分為4個階段： 局部孕育階段(1954年之前)、加速發展階段(1955—1974年)、急速擴張階段(1975—1984年)和控制減緩階段(1985年至今)[3]。截至2012年，華北平原沉降速率大于20 mm/a的區域面積約占其總面積的12.5%，其中90%的沉降區位于北京、天津、河北三地，局部地區累計沉降量超過3.4 m[4]。京津冀協同發展地質調查中將地面沉降列為京津冀一體化進程中需要關注的重大地質問題之一。對其地面沉降進行監測，認識其分布規律、演化特征，對京津冀城市和重大工程規劃建設，為實現京津冀協同發展的重大國家戰略提供地質環境保障。

合成孔徑雷達干涉測量技術(interferometric synthetic aperture Radar，InSAR)已被廣泛應用于地表形變監測中。與水準測量、分層標、全球定位系統(global positioning system，GPS)測量技術相比，InSAR技術具備大范圍、全天時全天候獲取高時空分辨率地表形變的能力。多時相InSAR(multi-temporal InSAR，MT-InSAR)技術的發展(如永久散射體干涉(persistent scatterer InSAR，PS-InSAR)、小基線 (small baseline subsets，SBAS)等)更是克服了大氣延遲、時空失相干對傳統InSAR方法的影響[5-6]。目前，MT-InSAR技術已被廣泛應用于京津冀地面沉降研究中，但這類研究多關注北京[7-11]、天津[12-15]、滄州[16-18]和廊坊[19-20]等城市區域，通過多源SAR數據獲取其地面沉降信息，而對京津冀平原地面沉降的分布、范圍及沉降程度缺乏整體性認識。針對單軌SAR數據源幅寬對區域地面沉降監測范圍的限制問題，本研究基于多軌SAR數據融合技術，獲取了2016—2018年京津冀平原地面沉降信息； 為保證InSAR監測結果的可靠性，選用同時期水準測量數據對InSAR結果進行精度驗證，并對相鄰軌道結果進行一致性檢驗； 并結合GIS空間分析、剖面分析等方法，獲取InSAR監測時段內京津冀平原區地面沉降演化特征。

1 研究區概況及其SAR數據源

京津冀地區(N36°05′～42°37′，E113°11′～119°45′)位于我國華北地區，包含北京、天津和河北省的11個地級市，屬暖溫帶大陸性季風型氣候，降水多集中在7—8月。該地區地形條件復雜，總體地勢為西北高、東南低，自西北向東南依次分布壩上高原、燕山—太行山山地和東南部平原[21]。京津冀地區存在多種環境地質問題，其中以地面沉降等為例的緩變性地質主要發生在平原區。京津冀平原區是華北平原的重要組成部分，自山麓到海岸，按成因和形態特征可分為第四紀山前沖積、洪積傾斜平原(山前平原)，中部沖積、湖積多層疊加平原(中部平原)和東部沖積、湖積夾海積的濱海平原(濱海平原)[22]。截至2015年，京津冀平原區累積沉降量大于200 mm的面積約為6.4萬km2[21]，約占平原區面積的73%，沉降速率大于50 mm/a的嚴重沉降區面積占全國嚴重沉降區面積的92%(http: //www.tianjin.cgs.gov.cn)。

Sentinel-1A衛星是歐洲空間局于2014年發射的重訪周期為12 d的C波段對地觀測衛星。該衛星使用近極地太陽同步軌道，軌道高度、傾角和周期分別為693 km，98.18°和99 min，具多種成像模式和極化方式(https: //sentinel.esa.int/)。單軌SAR影像受其幅寬限制，監測范圍有限，為獲取京津冀平原區的地面沉降信息，本次研究共選用了3軌時間跨度為2016—2018年的Sentinel-1A數據(Track 40，Track 142和Track 69)。所用3軌Sentinel-1A數據均為寬幅干涉模式(interferometric wide，IW)獲取的單視復數(single look complex，SLC)升軌存檔影像，其極化方式均為垂直極化(vertical-vertical，VV)，空間分辨率為5 m×20 m，單幅幅寬達250 km。影像覆蓋范圍及主要參數分別如圖1和表1所示。

圖1 研究區概況示意圖Fig.1 Location of study area

表1 Sentinel-1A影像信息Tab.1 Satellite information for the data used in this study

2 研究方法

2.1 MT-InSAR處理過程

通過SARPROZ軟件(https: //www.sarproz.com/)對3軌Sentinel-1A影像分別進行PS-InSAR處理，以獲取每軌SAR影像范圍內時序地面沉降信息。處理過程包括以下幾個主要步驟： 主影像選取，數據集配準及重采樣，基線建立，數字高程模型(digital elevation model,DEM)模擬，生成差分干涉圖，選取PS候選點，相位解纏，大氣相位估計和去除，重新選取PS點和估計PS點形變速率[23-24]。由前人研究成果已知，京津冀地區水平向形變值約為1～3 mm/a[25-28]，因此在本次研究中忽略水平向形變結果，通過式(1)將視線向(line-of-sight，LOS)InSAR形變結果轉為垂向，并將經過處理后的垂向結果用于與水準結果的精度驗證和后續的SAR影像結果拼接中，即

dv=dLOS/cosθ

(1)

式中：dv為垂向形變結果，mm；dLOS為LOS向形變結果，mm；θ為雷達LOS向與垂向的夾角。

2.2 多軌SAR數據拼接

為了獲取大區域的地面形變信息，需對上文所得單軌SAR影像結果進行拼接。首先對各軌SAR影像垂向形變結果進行了一致性檢驗和精度驗證，在保證數據精度可靠的情況下，對相鄰軌道數據進行融合。在相鄰軌道沉降速率融合中，Ge等[29]通過引入與參考點距離權重的關系進行整體平差; 孫赫[30]通過影像重疊區內PS點形變量差值眾數對相鄰軌道數據結果進行校正; 熊思婷等[31]比較了區塊法和插值法在求解異軌重疊區的形變差中的應用; 張永紅等[32]通過大量的水準觀測數據對SAR數據結果進行了平差控制。本次實驗中，依據重疊區內同名高相干點的形變速率一致準則[32-33]，計算平差權重。

由于受不同軌道衛星成像幾何差異，不同軌道SAR數據所得LOS向地表形變結果間存在系統性偏差。為消除不同視角投影的影響，首先依據式(1)將InSAR監測所得各軌LOS向形變值轉為垂向形變值。如圖1所示，Track 142與Track 69和40間都具有重疊區，因此在多軌SAR數據融合時，將Track 142形變結果作為基準。Track 142與Track 40影像重疊區域內，雷達LOS向與垂向的夾角分別為33.72° 和43.75°； Track 142與Track 69影像重疊區內LOS向與垂向夾角分別為43.77°和33.67°。根據上述θ值，將重疊區域內每軌影像的各PS點LOS向結果依據式(1)分別投影到垂直向上，選取的PS點的相干性均在0.9以上。之后，選用最鄰近法對重疊區內各軌PS點進行同名點匹配。最后，采用最小二乘法(least square method, LSM)構建相鄰軌道重疊區域的平差方程，并分別對Track 40和Track 69垂向形變結果進行偏移量補償。該方法是通過實際觀測值與模型計算值的誤差的平方和最小來尋找最優的匹配函數。假設相鄰軌道分別為x和y，其中x為參考基準，y為需要進行平差的數據，相鄰軌道之間存在n個重合點，基于LSM的原理，則

(2)

此時，通過對δ求極限，即可求得a和b的值。最終，獲得多軌SAR數據融合結果。

3 結果及精度驗證

3.1 京津冀2016—2018年地面沉降監測結果分析

為獲取京津冀平原區地面沉降信息，基于MT-InSAR技術分別獲取了3軌Sentinel-1A數據源地面沉降信息，后對相鄰軌道數據進行融合獲取了京津冀平原區內2016—2018年地面沉降速率(圖2)。為了顯示清晰，對結果進行了柵格化處理，像元大小為1 km。InSAR監測結果顯示， 2016—2018年間京津冀平原區內最大沉降速率達164 mm/a，分布在北京朝陽—通州沉降區和廊坊—天津沉降區內。研究區內沉降分布范圍廣，且空間分布不均，在北京、天津及河北省各地級市均有不同程度的地面沉降問題，且三地沉降區有連片趨勢。其中北京平原區沉降較為嚴重的地區主要分布在北京東南部朝陽、通州等地，且北京東南部沉降區已擴張到廊坊三河市界內。天津、廊坊兩地沉降較為嚴重的區域分布在兩地交界處。河北平原沉降較為嚴重的區域主要分布在其中東部地區，且該地區沉降區呈南北貫通趨勢。研究區內共獲得了1 426 967個PS點，其中2016—2018年沉降速率大于20, 50, 100 mm/a的點數分別占總點數的53.9%, 23.9%和3.6%。沉降速率大于20, 50, 100 mm/a的面積分別為3.3萬, 1.4萬和0.14萬km2。據統計，2016年沉降速率大于20, 50, 100 mm/a的點數分別占總點數的53.8%, 23.7%和3.6%； 到2017年沉降速率大于20和50 mm/a的點數略有增長，分別占總點數的54.4%和23.8%，沉降速率大于100 mm/a的點數占比則降至3.5%。2018年沉降速率大于20和50 mm/a的PS點數與2017年比有所增長，占比達57.2%和24.0%，沉降速率大于100 mm/a的點數與2017年相比保持穩定。

圖2 2016—2018年京津冀平原區PS-InSAR地表形變速率Fig.2 Mean displacement velocities throughout the Beijing-Tianjin-Hebei derived from the Sentinel-1A data by using PS-InSAR from 2016 to 2018

研究區內地面沉降大于50 mm的區域主要分布在北京東南部，廊坊北部，廊坊東部到天津西部，天津南部，保定東部到廊坊西部，保定西部，唐山—秦皇島沿海區域，以及河北平原東南部區域。對監測時段內研究區各子平原內沉降速率大于50 mm/a的嚴重沉降區面積進行統計發現，山前沖積洪積平原沉降量大于50 mm的面積有輕微增長(2016—2018年)，由2 937 km2變為2 982 km2，漲幅為1.5%。2016—2018年中部平原沉降量大于50 mm的面積的漲幅為0.3%，基本維持穩定。2016—2018年，嚴重沉降區的面積分別為10 547 km2，10 574 km2和10 577 km2。與2016年相比，濱海平原嚴重沉降區的面積在2018年增長較為明顯，由382 km2變為577 km2，漲幅達33.8%。其中唐山—秦皇島沿海沉降區范圍明顯擴大(如圖3所示)。為進一步分析研究時段內地面沉降發展過程，對沉降較嚴重區

(a) 2016年 (b) 2017年 (c) 2018年

圖3 2016—2018年各年度InSAR平均形變速率Fig.3 Mean displacement velocities derived by InSAR between 2016 and 2018

域沿東西、南北方向分別進行剖面分析，剖面線位置如圖3(c)中虛線所示。2016—2018年，不同地面沉降區的沉降空間演化存在顯著差異(如圖4所示)。其中位于唐山沉降區的TS剖面線處(圖4(e)和(f)所示)的沉降速率在監測時段呈增長趨勢。而位于北京、廊坊—天津、保定—滄州—衡水、衡水—邢臺、邯鄲沉降區的BJ，LT，BCH，HX和HD剖面處2016—2018年沉降速率基本保持一致，沉降中心發育較為穩定。綜上，InSAR監測時段內山前平原和中部平原地面沉降發展相對穩定，濱海平原地面沉降范圍及沉降速率在監測時段內增長明顯，尤其是唐山沿海地區，沉降速率漲幅最為明顯。

(a) BJ西—東 (b) BJ北—南

(c) LT西—東 (d) LT北—南

(e) TS西—東 (f) TS北—南

(g) BCH西—東 (h) BCH 北—南

(i) HX西—東 (j) HX北—南

(k) HD西—東 (l) HD北—南

圖4 剖面時間序列InSAR形變速率結果Fig.4 Profiles of the InSAR displacement velocities

為了獲取2016—2018年研究區內地面沉降時序特征，選取了8個不同位置上的沉降特征點，各沉降點的位置如圖3(b)所示，其對應時序特征如圖5所示。北京朝陽金盞、天津王慶坨、保定大崔營、保定高陽、邢臺南宮和邯鄲臨漳等沉降區，形變特征以線性下沉為主，最大累積沉降量超過300 mm。 在唐山豐南、樂亭沉降區，區域時序沉降具有明顯的季節性特征，在2016—2018年沉降呈下降趨勢，如圖5(g)和(h)所示，但在2016年9月—2017年1月和2017年9月—2018年1月沉降減緩甚至出現回彈現象，之后地表繼續下沉，至2018年10月地面累計沉降量達180～254 mm。

(a) 北京朝陽 (b) 天津王慶坨 (c) 保定大崔營 (d) 保定高陽

(e) 邢臺南宮 (f) 邯鄲臨漳 (g) 唐山豐南 (h) 唐山樂亭圖5 研究區內沉降特征點對應時序特征Fig.5 Accumulative time-series deformation revealed by the InSAR technique

3.2 精度驗證

為保證InSAR數據結果的準確性，選用了26個2017年9月—2018年9月的水準觀測值對Sentinel-1A數據所得垂向形變速率進行精度驗證。水準數據位置如圖1所示，位于Track 142和Track 69影像覆蓋范圍內，因此通過水準測量數據對這2軌影像所得形變結果進行驗證。在驗證過程中，以水準點為圓心，100 m為半徑建立緩沖區，將緩沖區內所有PS點的平均值視為該水準點對應的InSAR結果。驗證結果如圖6和表2所示。

(a) Track 142 (b) Track 69圖6 InSAR結果與水準測量結果精度驗證Fig.6 Comparison between the InSAR measurements and levelling data

表2 InSAR所得形變信息與水準測量結果比對Tab.2 Comparison of the mean subsidence rate between the Sentinel-1A PS and leveling data

其中Track 142所得研究區內垂向形變結果與水準測量結果相比，誤差最大值為11.6 mm/a，最小值為0.2 mm/a，均方根誤差為5.5 mm/a。對2組數據進行線性回歸分析，R達0.97。Track 69所得研究區內垂向形變結果與水準測量結果相比，誤差最大值為13.5 mm/a，最小值為1.0 mm/a，均方根誤差為7.2 mm/a，R達0.98。驗證結果表明InSAR監測結果與水準測量結果在同時間內具有較好的吻合度，表明InSAR結果的可靠性。

為進一步驗證不同軌SAR影像結果的準確性，對其進行了交叉驗證。由于Track 142的影像與其他2軌影像之間都具有重疊區，將其他2軌數據所得地表形變結果分別于Track 142所得結果進行比對。結果如圖7所示，由Track 142與Track 69所得形變結果之間的相關性達0.96，由Track 142與Track 40所得形變結果之間的相關性為0.97。多軌SAR自相關驗證結果表明了本次研究所用多軌InSAR數據結果的可靠性。

(a) Track 69與Track 142 (b) Track 40與Track 142圖7 多軌SAR影像形變速率交叉驗證結果Fig.7 Consistency between the vertical displacement rates derived from different datasets

4 結論

本研究基于多軌SAR數據融合技術，獲取了2016—2018年京津冀平原區地表形變信息。不同軌SAR數據集交叉驗證和水準驗證結果均表明本次研究InSAR測量所得地面沉降信息的可靠性，滿足地面沉降監測及其演化規律的研究。

從沉降場演變過程來看，2016—2018年京津冀平原區內沉降速率大于20 mm/a和50 mm/a的PS點數占總點數的比例均呈增長趨勢，沉降速率大于100 mm/a的點數占總點數的比例基本維持不變(3.6%，3.5%，3.5%)。InSAR監測時段內，濱海平原范圍內地面沉降發展較為迅速，其中2016—2018年，該區域地面沉降嚴重區(年沉降量>50 mm)面積增長195 km2，漲幅達33.8%。沉降區時序剖面結果顯示，唐山沉降區呈加速發展趨勢，其余沉降區發展較為穩定。InSAR所得地面沉降時序結果表明，研究時段內，京津冀平原區內地面沉降時序特征以線性下降為主，且唐山等沿海地區沉降區內地面沉降時序具有明顯的季節性波動特征。

研究成果一定程度上填補了京津冀平原區地面沉降整體演化特征研究較少的空缺，為京津冀地區地面沉降防治、調控，城市地下空間安全等提供科學依據。