時序InSAR解譯2017～2020年北京地面沉降時空變化

張雙成 許 強 羅 勇 雷坤超,4 牛玉芬 龐校光

1 長安大學地質工程與測繪學院，西安市雁塔路126號，7100542 地理信息工程國家重點實驗室，西安市雁塔中路1號，7100543 北京市水文地質工程地質大隊，北京市西四環北路123號，1001954 中國科學院地質與地球物理研究所，北京市北土城西路19號，1000295 河北工程大學礦業與測繪工程學院，河北省邯鄲市太極路19號，056038

地面沉降是我國平原地區主要的區域性環境地質災害，可能對人民的生命和財產安全造成危害[1]。北京是中國地面沉降最嚴重的地區之一[2]，北京的沉降現象最早出現在20世紀50年代，并在近幾十年來迅速擴散[3]。截至2009年底，北京市最大年沉降量已達137.51 mm，最大累積沉降量為1 163 mm[4-5]。因此，監測和分析北京市地面沉降活動的時空變化，對于防災減災和城市化的可持續發展具有重要意義。

InSAR具有大范圍、高精度、全天候、全天時的特點，已逐漸成為監測城市地面沉降的主要技術手段[6]。Hu等[7]利用短基線集(small baseline subset，SBAS)方法獲取2003～2010年北京地面沉降時空演變特征，結果表明，研究區內存在3個大型沉降漏斗；周呂等[8]運用SBAS方法獲取2007～2010年北京地區的地表沉降分布，結果表明，不均勻沉降較為明顯，各個沉降漏斗逐漸連成一片，且有向東發展的趨勢；Zhang等[9]基于一種改進的多時相InSAR方法獲取1992～2014年京津冀地區3個時間跨度的沉降速率，結果表明，北京市1992～2000年間只有部分小規模沉降點，2003～2010年地面沉降范圍急劇擴大，2012～2014年間地面沉降持續擴散。南水北調工程于2014-12開始向北京輸送水源，截至2020-09，平原區地下水埋深22.49 m，較2015年底的25.75 m上升了3.26 m，地下水水位連續5 a持續回升。南水北調工程為北京市提供新的水源，改變了2015年以來北京市地面沉降演變的格局。因此，不斷更新北京地面沉降速率以預測危急情況，并在必要時采取緩解措施是很有必要的。

選取覆蓋北京地區的85景C波段的Sentinel-1數據集，采用SBAS技術獲取研究區域2017～2020年的沉降分布情況和沉降速率場，并通過對比分析2018年和2019年的年形變速率來探究北京地面沉降的時空變化規律，為城區的沉降監測和預警提供理論依據。

1 研究區域及實驗數據

1.1 研究區域

北京市位于華北平原的西北邊緣，中心位置為39°54′N、116°23′E。其西部、北部均為山區，東南部為向渤海緩慢傾斜的平原。山區面積為10 072 km2(占總面積的61.4%)，平原面積為6 388 km2(占38.6%)；山區平均海拔為1 000～1 500 m，平原地區平均海拔為20～60 m。北京屬暖溫帶半濕潤季風型大陸性氣候，夏季炎熱多雨，冬季寒冷干燥；年平均氣溫為10～12 ℃，年平均降雨量約為600 mm。

1.2 實驗數據

選取2017-06～2020-06覆蓋實驗區的Sentinel-1A序列數據集，共包含85景升軌影像。采用覆蓋研究區的30 m分辨率SRTM(shuttle radar topography mission)DEM作為參考高程數據，用于輔助去除地形相位。同時獲取研究區域內8個GPS站點數據，用于InSAR結果的對比驗證。

2 數據處理方法

SBAS方法的原理是通過設定一定的空間和時間閾值，將所有的SAR數據組成若干個小集合，集合內基線較小，集合間基線較大。最后通過集合內的最小二乘求解和集合間的奇異值分解方法，得到整個時間序列地表形變信息的聯合求解結果[10-11]。

假設獲取N+1幅SAR影像，其獲取的時間序列為(t0,…,tN)，依據干涉基線組合可生成M幅干涉圖，則在像元(x,r)處由tA、tB兩個時刻SAR影像生成的差分干涉相位可表示為：

δφj(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

φdef,j(x,r)+φtopo,j(x,r)+φatm,j(x,r)+

φnoise,j(x,r)，j=1,…,M

(1)

式中，(x,r)為方位向和距離向坐標，φ(tB,x,r)和φ(tA,x,r)分別為生成干涉相位的影像相位，φdef,j(x,r)為tA時刻至tB時刻間衛星視線向的形變相位，φtopo,j(x,r)為參考DEM不精確引起的地形相位誤差，φatm,j(x,r)為大氣相位誤差，φnoise,j(x,r)為噪聲相位。

由式(1)可以得到一個方程組，包含N個未知數的M個方程：

δφ=Aφ

(2)

式中，δφ為M幅干涉相位構成的矩陣，φ為N幅SAR影像上的待求形變相位構成的矩陣；A為M×N矩陣，其每一行對應一個干涉對。當小基線子集個數L=1時，A為列滿秩矩陣，其最小二乘解為：

(3)

當小基線子集個數L>1時，式(2)是秩虧的，秩虧數為N-L+1，可以對A進行奇異值分解，求出形變相位φ最小范數意義下的最小二乘解[12]。

3 時序InSAR解譯地面沉降及精度評定

3.1 地面沉降空間分布

利用SBAS方法對覆蓋北京市的Sentinel-1影像數據進行地表形變特征信息的提取，獲取北京平原地區沿衛星雷達視線向上的地表年平均形變速率，如圖1所示(負值表示地表正在遠離衛星，正值表示地表正在靠近衛星)。在解纏相位時，需要選擇一個像素作為解纏的起始點(即參考點)，該點位于高相干區域的局部最高處，變形速度場中各監測點的平均沉降速度均與該點有關，該點的位置如圖1所示。

圖1 北京2017-06～2020-06年平均形變速率Fig.1 Annual mean deformation velocity of Beijing from June 2017 to June 2020

由圖1可見，研究區域內不均勻沉降較為突出。2017-06～2020-06沉降主要發生在昌平(以下簡稱CP)、順義(SY)、通州(TZ)、大興(DX)以及北京和天津廊坊交界處(LF)，最大年形變速率在SY區域(為-111.3 mm/a)。CP沉降區主要的形變速率范圍為-20～-60 mm/a，平均值為-37.8 mm/a，最大形變速率達-84 mm/a；SY沉降區主要的形變速率范圍為-20～-70 mm/a，平均值為-43.1 mm/a，最大形變速率達-111.3 mm/a；TZ沉降區主要的形變速率范圍為-25～-65 mm/a，平均值為-45.3 mm/a，最大形變速率達-89.7 mm/a；DX沉降區主要的形變速率范圍為-20～-45 mm/a，平均值為-31.9 mm/a，最大形變速率達-63.6 mm/a；LF沉降區主要的形變速率范圍為-20～-40 mm/a，平均值為-30.2 mm/a，最大形變速率達-59.5 mm/a。

研究區域2017-06～2020-06期間根據C波段的Sentinel數據獲取的累積形變如圖2所示。2017～2020年期間，TZ沉降區最大累積量約為255 mm，SY沉降區最大累積量約為287 mm，CP沉降區最大累積量約為241 mm，DX和LF沉降區相對較小，最大累積量分別約為157 mm和175 mm。

圖2 北京2017-06～2020-06累積形變Fig.2 Cumulative deformation of Beijing from June 2017 to June 2020

3.2 地面沉降精度評定

對2017-06～2020-06視線向上的平均形變速率進行誤差評估，通過計算速度線性擬合的偏差獲取其對應的均方根誤差RMSE(圖3)。由圖3可見，總體上RMSE不超過4 mm/a，平均值為1.4 mm/a，這對于InSAR時間序列結果來說精度較高。北京市中心的RMSE較小，小于1 mm/a，形變區和西北區域的RMSE為2～4 mm/a。推測造成此誤差分布格局的原因是，北京市中心相干性較好；由于形變區沉降明顯不均勻及西北區處于山區，相干性較差。

圖3 平均形變速率的RMSEFig.3 RMSE of the mean velocity results

獲取2017～2019年期間8個GPS站點的時序形變監測值，總共3期數據，各期起始監測時間分別為2017-08-24、2018-08-25和2019-07-23，各期結束時間分別為2017-10-23、2018-10-31和2019-09-03。為了將InSAR結果與GPS觀測值在時間段上進行統一，采用2017-10-23～2019-09-03期間內的影像數據生成InSAR視線方向上的平均速率圖。

提取各個GPS點對應的InSAR結果與視線方向上的GPS形變量進行對比，結果如表1和圖4所示，最大絕對誤差為4.45 mm/a，最小絕對誤差為-5.48 mm/a，各站測量值與InSAR結果之間的RMSE為3.04 mm/a，表明兩者之間具有良好的一致性。GPS與InSAR結果的相關系數為0.96，說明InSAR獲取的視線方向的平均速率的精度較高。

表1 InSAR結果與GPS觀測值對比

圖4 InSAR結果與GPS觀測值對比Fig.4 Comparison between time series ofInSAR and GPS measurements

3.3 地面沉降時空變化

為了探究研究區地面沉降的時空變化規律，分別利用2018年和2019年期間影像數據組成干涉對獲取研究區2018年和2019年視線向上的形變速率，如圖5所示(負值表示地面位移方向為遠離衛星，正值表示地面位移方向為朝向衛星)。由圖可見，北京地區2018年與2019年的沉降區空間分布具有較高的一致性，但從沉降速率及沉降范圍來看，北京地區出現了減速沉降的現象。例如，2018年北京地區最大沉降速率為137 mm/a，而2019年最大沉降速率為110.1 mm/a，下降了26.9 mm/a。2018年和2019年形變速率圖柵格數量分別為9 090 995和9 177 465個，之間的差量相比于總數可以忽略不計。對整幅速率圖進行統計分析可知，2018年平均沉降速率為10.2 mm/a，而2019年平均沉降速率為5.6 mm/a，下降了4.6 mm/a。由此可知，北京市整體沉降減緩。

圖5 北京視線向平均形變速率Fig.5 Mean deformation velocity in LOS direction of Beijing

2018年和2019年沉降速率變化如圖6所示(柵格大小為46.60 m×69.64 m)，考慮到SBAS-InSAR結果的不確定性，認定沉降速率大于10 mm/a的區域為沉降區。由圖可見，在各個沉降速率范圍內，2019年形變速率圖包含的柵格數量都比2018年的少。在沉降速率為10～20 mm/a的范圍內，減少的柵格數量最多，達到18萬個，減少幅度為13.62%；其次是沉降速率為20～30 mm/a的范圍，減少16萬個；在沉降速率為60～70 mm/a、70～80 mm/a、80～90 mm/a和90～140 mm/a的范圍內，共計減少15萬個柵格，減少幅度較大，分別為74.72%、75.41%、81.96%和75.54%。由此可知，在各個沉降范圍內沉降的面積都在減小。

圖6 2018年和2019年沉降速率變化Fig.6 The deformation velocity change in 2018 and 2019

將2018年和2019年形變速率圖進行差異化分析，獲取這2 a速度變化幅度的統計數據，如圖7所示(負值表示2018～2019年沉降減緩，正值表示2018～2019年沉降加速)?？紤]到SBAS-InSAR結果的不確定性，將速度變化為-10～10 mm/a的柵格視為穩定點。圖7中處于沉降減緩的區域占總面積的10.3%，而加速區域僅占0.6%，減緩的面積是加速面積的16倍。由此可知，沉降減緩的面積遠大于沉降加快的面積。

圖7 2018年和2019年形變速率的差異Fig.7 Difference indeformation velocity between 2018 and 2019

分析CP、SY、TZ、DX、LF五處沉降區(圖5中沉降區對應的藍色長方形區域)的沉降狀況，并對局部沉降速度變化進行調查(圖8)。由圖可見， SY、TZ、DX、LF 四處沉降區2018～2019年的最大形變速率和平均形變速率都在下降，其中TZ區下降最為明顯，最大沉降速率由113.0 mm/a下降到77.8 mm/a，平均沉降速率由50.7 mm/a下降到23.7 mm/a；CP區沉降仍在加速，最大沉降速率由88.4 mm/a提高到110.0 mm/a，平均沉降速率由37.5 mm/a提高到49.3 mm/a。由此可知，SY、TZ、DX、LF四處沉降區沉降均在減緩，而CP區沉降仍在加快，其中TZ區沉降減緩的幅度大于CP區沉降加快的幅度。

圖8 5處沉降區形變速率統計Fig.8 Statistics of deformation velocity in five subsidence areas

4 結 語

基于2017-06～2020-06共85景Sentinel-1影像數據，利用SBAS-InSAR技術對北京地表形變進行監測，分析2017～2020年北京地面沉降時空變化。得到的SBAS-InSAR結果與GPS觀測值吻合良好，R2為0.96，RMSE為3.04 mm/a。結果表明：

1)北京不均勻沉降較為突出，2017-06～2020-06地表形變呈現5處沉降區，最大年沉降速率為-111.3 mm/a，最大累積量為287 mm。

2)綜合2018年和2019年形變速率圖可以看出，北京整體地面沉降情況在減緩，在各個沉降范圍內的沉降速率都在減緩，且沉降減緩的面積遠大于沉降加速的面積。

3)對局部沉降調查發現，北京市5處沉降區中除1處仍然在加速外，其他4處沉降速度均在減緩。

由于本次實驗獲取的SAR數據集時間較短，可能會對實驗結果造成一定的影響。隨著Sentinel-1衛星進入業務化運行，SAR數據的不斷累積可為更好地揭示研究區域地面沉降時空演變規律提供豐富的研究資料。地下水是維持土壤應力平衡的重要因素，研究地表沉降與地下水變化的關系，可進一步全面解譯南水進京后北京地面沉降的發展趨勢。

致謝：感謝歐空局提供Sentinel-1雷達影像數據、美國航空航天局提供SRTM DEM地形數據以及北京市水文地質工程地質大隊提供部分GPS監測數據。