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

張雙成 許 強(qiáng) 羅 勇 雷坤超,4 牛玉芬 龐校光

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安市雁塔路126號(hào)，7100542 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，西安市雁塔中路1號(hào)，7100543 北京市水文地質(zhì)工程地質(zhì)大隊(duì)，北京市西四環(huán)北路123號(hào)，1001954 中國科學(xué)院地質(zhì)與地球物理研究所，北京市北土城西路19號(hào)，1000295 河北工程大學(xué)礦業(yè)與測繪工程學(xué)院，河北省邯鄲市太極路19號(hào)，056038

地面沉降是我國平原地區(qū)主要的區(qū)域性環(huán)境地質(zhì)災(zāi)害，可能對人民的生命和財(cái)產(chǎn)安全造成危害[1]。北京是中國地面沉降最嚴(yán)重的地區(qū)之一[2]，北京的沉降現(xiàn)象最早出現(xiàn)在20世紀(jì)50年代，并在近幾十年來迅速擴(kuò)散[3]。截至2009年底，北京市最大年沉降量已達(dá)137.51 mm，最大累積沉降量為1 163 mm[4-5]。因此，監(jiān)測和分析北京市地面沉降活動(dòng)的時(shí)空變化，對于防災(zāi)減災(zāi)和城市化的可持續(xù)發(fā)展具有重要意義。

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

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

1 研究區(qū)域及實(shí)驗(yàn)數(shù)據(jù)

1.1 研究區(qū)域

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

1.2 實(shí)驗(yàn)數(shù)據(jù)

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

2 數(shù)據(jù)處理方法

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

假設(shè)獲取N+1幅SAR影像，其獲取的時(shí)間序列為(t0,…,tN)，依據(jù)干涉基線組合可生成M幅干涉圖，則在像元(x,r)處由tA、tB兩個(gè)時(shí)刻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)為方位向和距離向坐標(biāo)，φ(tB,x,r)和φ(tA,x,r)分別為生成干涉相位的影像相位，φdef,j(x,r)為tA時(shí)刻至tB時(shí)刻間衛(wèi)星視線向的形變相位，φtopo,j(x,r)為參考DEM不精確引起的地形相位誤差，φatm,j(x,r)為大氣相位誤差，φnoise,j(x,r)為噪聲相位。

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

δφ=Aφ

(2)

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

(3)

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

3 時(shí)序InSAR解譯地面沉降及精度評定

3.1 地面沉降空間分布

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

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

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

研究區(qū)域2017-06～2020-06期間根據(jù)C波段的Sentinel數(shù)據(jù)獲取的累積形變?nèi)鐖D2所示。2017～2020年期間，TZ沉降區(qū)最大累積量約為255 mm，SY沉降區(qū)最大累積量約為287 mm，CP沉降區(qū)最大累積量約為241 mm，DX和LF沉降區(qū)相對較小，最大累積量分別約為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視線向上的平均形變速率進(jìn)行誤差評估，通過計(jì)算速度線性擬合的偏差獲取其對應(yīng)的均方根誤差RMSE(圖3)。由圖3可見，總體上RMSE不超過4 mm/a，平均值為1.4 mm/a，這對于InSAR時(shí)間序列結(jié)果來說精度較高。北京市中心的RMSE較小，小于1 mm/a，形變區(qū)和西北區(qū)域的RMSE為2～4 mm/a。推測造成此誤差分布格局的原因是，北京市中心相干性較好；由于形變區(qū)沉降明顯不均勻及西北區(qū)處于山區(qū)，相干性較差。

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

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

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

表1 InSAR結(jié)果與GPS觀測值對比

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

3.3 地面沉降時(shí)空變化

為了探究研究區(qū)地面沉降的時(shí)空變化規(guī)律，分別利用2018年和2019年期間影像數(shù)據(jù)組成干涉對獲取研究區(qū)2018年和2019年視線向上的形變速率，如圖5所示(負(fù)值表示地面位移方向?yàn)檫h(yuǎn)離衛(wèi)星，正值表示地面位移方向?yàn)槌蛐l(wèi)星)。由圖可見，北京地區(qū)2018年與2019年的沉降區(qū)空間分布具有較高的一致性，但從沉降速率及沉降范圍來看，北京地區(qū)出現(xiàn)了減速沉降的現(xiàn)象。例如，2018年北京地區(qū)最大沉降速率為137 mm/a，而2019年最大沉降速率為110.1 mm/a，下降了26.9 mm/a。2018年和2019年形變速率圖柵格數(shù)量分別為9 090 995和9 177 465個(gè)，之間的差量相比于總數(shù)可以忽略不計(jì)。對整幅速率圖進(jìn)行統(tǒng)計(jì)分析可知，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結(jié)果的不確定性，認(rèn)定沉降速率大于10 mm/a的區(qū)域?yàn)槌两祬^(qū)。由圖可見，在各個(gè)沉降速率范圍內(nèi)，2019年形變速率圖包含的柵格數(shù)量都比2018年的少。在沉降速率為10～20 mm/a的范圍內(nèi)，減少的柵格數(shù)量最多，達(dá)到18萬個(gè)，減少幅度為13.62%；其次是沉降速率為20～30 mm/a的范圍，減少16萬個(gè)；在沉降速率為60～70 mm/a、70～80 mm/a、80～90 mm/a和90～140 mm/a的范圍內(nèi)，共計(jì)減少15萬個(gè)柵格，減少幅度較大，分別為74.72%、75.41%、81.96%和75.54%。由此可知，在各個(gè)沉降范圍內(nèi)沉降的面積都在減小。

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

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

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

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

圖8 5處沉降區(qū)形變速率統(tǒng)計(jì)Fig.8 Statistics of deformation velocity in five subsidence areas

4 結(jié) 語

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

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

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

3)對局部沉降調(diào)查發(fā)現(xiàn)，北京市5處沉降區(qū)中除1處仍然在加速外，其他4處沉降速度均在減緩。

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

致謝：感謝歐空局提供Sentinel-1雷達(dá)影像數(shù)據(jù)、美國航空航天局提供SRTM DEM地形數(shù)據(jù)以及北京市水文地質(zhì)工程地質(zhì)大隊(duì)提供部分GPS監(jiān)測數(shù)據(jù)。