南水北調(diào)中線區(qū)域地面沉降SBAS-InSAR監(jiān)測(cè)研究

張雙成 余 靜 宋明鑫 張彬玲 樊茜佑 司錦釗 張雅斐

1 長安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安市雁塔路126號(hào)，710054 2 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，西安市雁塔路中段1號(hào)，710054 3 西安市勘察測(cè)繪院，西安市南二環(huán)東段29號(hào)，710054

地面沉降是指在自然因素和人為因素的作用下，由于地殼表層未固結(jié)土體壓縮變形而引起的區(qū)域性地面高程降低的環(huán)境地質(zhì)現(xiàn)象[1]，具有持續(xù)時(shí)間長、影響范圍大和多種因素共同作用等特點(diǎn)[2]。南水北調(diào)中線工程是世界上最大的跨流域調(diào)水工程，意在解決華北平原區(qū)域性缺水的問題，保障干渠的安全運(yùn)營是推進(jìn)受水區(qū)人民生活與經(jīng)濟(jì)發(fā)展水平的重要措施[3]。而工程沿線部分渠段受城市地下水開采、礦區(qū)開采及區(qū)域地質(zhì)構(gòu)造活動(dòng)等多種威脅，區(qū)域地質(zhì)結(jié)構(gòu)穩(wěn)定性較差[4]。為避免沿線城市地面沉降影響到中線工程正常通水運(yùn)行，對(duì)沿線區(qū)域進(jìn)行長時(shí)間、高精度的地表形變監(jiān)測(cè)具有重要意義。

與傳統(tǒng)監(jiān)測(cè)技術(shù)相比，合成孔徑雷達(dá)干涉測(cè)量InSAR技術(shù)具有獲取數(shù)據(jù)速度快、可覆蓋范圍大、全天時(shí)全天候監(jiān)測(cè)及精度高等優(yōu)勢(shì)，被廣泛用于各個(gè)領(lǐng)域[5]。由于受到時(shí)空失相干和大氣延遲等的影響，常規(guī)差分干涉測(cè)量D-InSAR技術(shù)的監(jiān)測(cè)精度較低，為突破D-InSAR技術(shù)的限制，小基線集技術(shù)SBAS-InSAR[6]應(yīng)運(yùn)而生。該技術(shù)大幅提高了形變監(jiān)測(cè)精度，并且在探測(cè)長期累積的緩慢形變方面表現(xiàn)出極大的潛力。

本文利用SBAS-InSAR技術(shù)對(duì)南水北調(diào)中線工程沿線2015-04～2020-11的地表形變進(jìn)行研究，得到沿線區(qū)域形變特征分布及時(shí)空演變規(guī)律，并結(jié)合有關(guān)資料著重探討了北京市地面沉降在南水進(jìn)京后的時(shí)空演變特征，其結(jié)果可為相關(guān)單位在防災(zāi)減災(zāi)決策方面提供參考。

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

1.1 SBAS-InSAR基本原理

2002年，Berardino等[6]首次提出小基線集SBAS技術(shù)，該方法廣泛應(yīng)用于長時(shí)間序列的地表形變監(jiān)測(cè)中，基本原理如下[7]：若有同一研究區(qū)域的N+1景SAR 影像，則可能得到M個(gè)干涉對(duì)，M需滿足：

(1)

設(shè)初始時(shí)刻為t0，任意時(shí)刻ti(i=1,…,N)相對(duì)于初始時(shí)刻的相位差是未知參數(shù)φ(ti)，干涉處理獲取的差分干涉相位δφ(tk)(k=1,…,M)為觀測(cè)量。若不考慮其他相位誤差，則第k景差分干涉圖的像元(x，y)地表形變相位δφdef可表示為：

δφdef=φ(tB,x,y)-φ(tA,x,y)=

(2)

式中，tA、tB(tA

將所有差分干涉圖進(jìn)行自由組合，φ為N幅SAR影像上的干涉相位值組成的矩陣，δφ為M幅差分干涉圖上相位組成的矩陣。則有：

(3)

(4)

將上式寫成矩陣形式：

Aφ=δφ

(5)

式中，A為M×N維矩陣。

(6)

當(dāng)基線組L>1時(shí)，式(5)秩虧，秩虧數(shù)為N-L+1，為此對(duì)系數(shù)矩陣A進(jìn)行奇異值分解，求出累積形變相位φ最小范數(shù)意義上的最小二乘解。

1.2 數(shù)據(jù)處理

設(shè)置一定的時(shí)空基線閾值，滿足條件的SAR影像自由組合生成干涉圖集。將DEM數(shù)據(jù)與主影像進(jìn)行配準(zhǔn)，模擬相位對(duì)去平后的干涉圖進(jìn)行地形相位去除，得到差分干涉圖。對(duì)差分干涉圖進(jìn)行濾波來降低由于失相干引起的噪聲等，并對(duì)濾波后得到的干涉圖進(jìn)行相位解纏。構(gòu)建矩陣方程，基于奇異值分解(SVD)方法估算解纏相位的高程誤差和平均形變速率，通過時(shí)間域和空間域?yàn)V波對(duì)殘余相位中的失相干噪聲相位和大氣延遲相位進(jìn)行去除，得到非線性形變速率。結(jié)合線性形變速率和非線性形變速得到不同影像時(shí)間序列間的地表形變速率結(jié)果。SBAS方法的數(shù)據(jù)處理技術(shù)路線如圖1所示。

圖1 SBAS數(shù)據(jù)處理技術(shù)路線Fig.1 Data processing technique flowchart by SBAS

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

2.1 研究區(qū)域概況

南水北調(diào)中線工程輸水干渠源自丹江口水庫，流經(jīng)河南、河北，最后到達(dá)北京、天津，總長1 432.8 km。自2014年底通水后，中線工程已全面進(jìn)入運(yùn)行階段，截至2021-07-19，南水北調(diào)中線一期工程自陶岔渠首累積調(diào)水入渠水量達(dá)400億m3，分別向河南省、河北省、天津市、北京市供水135億m3、116億m3、65億m3及68億m3，沿線直接受益人口增加至7 900萬人，比2015年通水1周年時(shí)的3 800萬受益人口增加1倍多。

2.2 實(shí)驗(yàn)數(shù)據(jù)來源

Sentinel-1A衛(wèi)星是歐空局繼ERS和Envisat系列衛(wèi)星后發(fā)射的環(huán)境監(jiān)測(cè)衛(wèi)星，基于C波段的成像系統(tǒng)采用條帶成像、干涉寬幅、超寬幅及波浪模式4種成像模式，其中干涉寬幅模式(IW)采用中等分辨率5 m×20 m獲取幅寬250 km的影像，成為陸地覆蓋的默認(rèn)模式。本文收集IW模式下Sentinel-1A數(shù)據(jù)，其覆蓋情況如圖2所示，詳細(xì)數(shù)據(jù)參數(shù)如表1所示。除此之外，利用歐空局提供的POD精密軌道數(shù)據(jù)提高SAR衛(wèi)星影像軌道精度，利用美國宇航局提供的分辨率為30 m的SRTM DEM數(shù)據(jù)去除地形相位。

3 中線區(qū)域地面沉降時(shí)空演變特征分析

由于研究區(qū)域較大，需進(jìn)行數(shù)據(jù)分塊處理。將每景影像經(jīng)SBAS-InSAR技術(shù)處理后的結(jié)果轉(zhuǎn)換到地理坐標(biāo)系下，在ArcGIS軟件中依據(jù)一定的規(guī)則進(jìn)行柵格鑲嵌拼接，得到中線工程沿線區(qū)域形變結(jié)果。

圖3為中線工程沿線區(qū)域2015-04～2020-11 LOS向的平均形變速率，由圖可見，整個(gè)沿線區(qū)域存在多處地面沉降，主要分布在京津冀地區(qū)，包括北京市、保定市、天津市、廊坊市、邢臺(tái)市及邯鄲市等。大部分沉降分布在河北省東南部，幾乎接連成片，最大地面沉降速率達(dá)139 mm/a，位于邢臺(tái)市南宮市沉降中心。北京市最大地面沉降速率達(dá)133 mm/a，位于朝陽區(qū)沉降中心。天津市西南部最大地面沉降速率為81 mm/a，相對(duì)較小。位于河北省東南部的地面沉降由于距離中線干渠有一定的距離，對(duì)輸水影響較小，但天津市支線經(jīng)過2個(gè)小范圍沉降區(qū)，應(yīng)引起重視。

圖3 中線沿線區(qū)域2015～2020年LOS向平均形變速率Fig.3 Mean deformation velocity in LOS of the middle route project from 2015 to 2020

由于覆蓋各個(gè)區(qū)域的影像日期不統(tǒng)一，在計(jì)算累積形變時(shí)以2016-12為時(shí)間基準(zhǔn)。圖4為以2016-12為時(shí)間起點(diǎn)，至2017-12、2018-12、2019-12、2020-11獲取的中線沿線區(qū)域LOS向累積形變分布。在此期間，整個(gè)中線沿線的地表累積形變與形變速率具有相似的空間分布格局，且隨著時(shí)間的推進(jìn)呈持續(xù)發(fā)展趨勢(shì)。最大累積沉降量位于北京市朝陽區(qū)沉降中心，2016-12～2017-12的最大沉降量為158 mm，2016-12～2018-12的最大累積沉降量為287 mm，相比于2016-12～2017-12，2016-12～2018-12的累積沉降量增量為129 mm；2016-12～2019-12的最大累積沉降量為402 mm，比2016-12～2018-12增加了115 mm；2016-12～2020-11的最大累積沉降量為476 mm，比2016-12～2019-12增加了74 mm。

圖4 中線沿線區(qū)域2016-12～2020-11 LOS向累積形變量Fig.4 Cumulative deformation in LOS of the middle route from December 2016 to November 2020

以天津支線主干渠兩側(cè)5 km緩沖區(qū)范圍為研究對(duì)象，利用SBAS-InSAR技術(shù)獲取覆蓋天津支線的2015-09～2020-11年平均形變速率(圖5)。沿線區(qū)域的最大形變速率為-94 mm/a，整條線貫穿地面沉降區(qū)域，主要沉降區(qū)有2個(gè)：沉降區(qū)Ⅰ位于保定市與廊坊市交界處，區(qū)內(nèi)最大平均形變速率為-94 mm/a；沉降區(qū)Ⅱ位于廊坊市與天津市交界處，區(qū)內(nèi)最大平均形變速率為-81 mm/a。地面沉降的不均勻分布容易對(duì)地下箱涵造成破壞，產(chǎn)生裂縫，從而導(dǎo)致滲水。相比于干渠，沉降區(qū)Ⅰ、Ⅱ內(nèi)沉降的不均勻分布較為明顯，應(yīng)對(duì)上述典型地段進(jìn)行實(shí)地調(diào)查，檢查箱涵的健康狀況，防患于未然。

圖5 天津支線2015-07～2020-11 LOS向平均形變速率Fig.5 Mean deformation velocity in LOS of the Tianjin branch line from July 2015 to November 2020

4 重點(diǎn)區(qū)域時(shí)序特征分析

4.1 北京市沉降結(jié)果及驗(yàn)證

為探究區(qū)域地面沉降在南水北調(diào)中線工程通水后的時(shí)空演變特征，對(duì)北京市的地面沉降進(jìn)行重點(diǎn)分析。北京市2015-11～2020-11 LOS向平均形變速率如圖6所示，由圖可見，北京市地面沉降空間分布差異性很大，沉降區(qū)域主要分布在昌平區(qū)、海淀區(qū)、順義區(qū)、朝陽區(qū)、通州區(qū)及大興區(qū)。各個(gè)行政區(qū)內(nèi)形成多個(gè)沉降漏斗，分別為HD(海淀)、CP(昌平)、SY(順義)、CY(朝陽)、TZ(通州)、DX(大興)，其中CY的年均沉降速率最大，達(dá)到133 mm/a，HD、CP、SY、TZ、DX的年均沉降速率分別為89.72 mm/a、59.07 mm/a、57.01 mm/a、95.67 mm/a和69.15 mm/a。

圖6 北京市2015-11～2020-11 LOS向平均形變速率Fig.6 Mean deformation velocity in LOS of the Beijing from November 2015 to November 2020

為評(píng)定InSAR監(jiān)測(cè)結(jié)果精度，獲取2017～2019年北京市10個(gè)GPS站點(diǎn)的監(jiān)測(cè)結(jié)果，依據(jù)GPS站點(diǎn)位置選擇InSAR結(jié)果的相應(yīng)點(diǎn)進(jìn)行對(duì)比。為保證時(shí)間的一致性，獲取該時(shí)段內(nèi)LOS向InSAR形變速率分布，將GPS監(jiān)測(cè)結(jié)果轉(zhuǎn)換到LOS向與其對(duì)應(yīng)的InSAR結(jié)果進(jìn)行對(duì)比分析，如表2(單位mm/a)所示。由表可見，兩者之間具有較高的一致性，互差為0～7 mm，GPS觀測(cè)值與InSAR結(jié)果之間的RMSE為2.47 mm/a，并表現(xiàn)出明顯的線性相關(guān)，最大線性相關(guān)系數(shù)R2達(dá)0.974，說明本次InSAR監(jiān)測(cè)結(jié)果精度較高。

4.2 地面沉降時(shí)空特征分析

為揭示北京市地面沉降形變特征和演化規(guī)律，獲取北京市2015-11-03～2020-11-24 LOS向時(shí)序累積形變量。圖7為2015-11-03～2020-11-24北京市地面沉降空間特征演化過程，圖8為該時(shí)段的累積形變分布。由圖7和8可見，北京市的地面沉降處于持續(xù)加重狀態(tài)，最大累積沉降量為697 mm，位于朝陽區(qū)。為對(duì)比沉降區(qū)域的形變特征是否具有相似趨勢(shì)，選取4個(gè)點(diǎn)(點(diǎn)位見圖8)并提取其線性形變時(shí)序。由圖9可見，最大累積沉降點(diǎn)位于F4，累積沉降量為507.14 mm，最小累積沉降點(diǎn)位于F2，累積沉降量為342.11 mm。4個(gè)特征點(diǎn)呈不同速度的非線性沉降趨勢(shì)，經(jīng)歷了從快速到逐漸減緩的沉降形變過程，由此可知，北京市的地面沉降正在減緩。

表2 InSAR結(jié)果與GPS觀測(cè)值比較

圖7 北京市2015-11-03～2020-11-24地面沉降空間特征演化過程Fig.7 The evolution of spatial characteristics of ground subsidence in Beijing from November 3, 2015 to November 24, 2020

圖8 北京市2015-11-03～2020-11-24累積形變Fig.8 Cumulative deformation in Beijing from November 3, 2015 to November 24, 2020

圖9 4個(gè)點(diǎn)的形變時(shí)序Fig.9 Time series deformation of 4 points

圖10為2016～2020年北京市LOS向年均形變速率分布。由圖可見，不同時(shí)段內(nèi)的地面沉降分布相似，研究時(shí)段內(nèi)沉降速率最大的區(qū)域始終位于朝陽區(qū)，2016年最大沉降速率為167 mm/a，2017年為175 mm/a，2018年為145 mm/a，2019年最大沉降速率為136 mm/a，2020年最大沉降速率持續(xù)減緩，變?yōu)?8 mm/a。

為進(jìn)一步分析地面沉降速率隨時(shí)間的演化過程，選取地面沉降分布區(qū)域內(nèi)的典型特征點(diǎn)。特征點(diǎn)的選取基于局部沉降區(qū)域中較為嚴(yán)重的沉降中心，即圖6中的HD、CP、SY、CY、TZ、DX，分別提取各點(diǎn)年均速率進(jìn)行對(duì)比分析。由圖11可見，CP的沉降速率變化最為明顯，由134.69 mm/a變?yōu)?9.16 mm/a；HD的沉降速率從132.95 mm/a變?yōu)?9.72 mm/a； SY的沉降速率從91.9 mm/a變?yōu)?8.29 mm/a；CY的沉降速率從165.94 mm/a變?yōu)?7.76 mm/a；TZ的沉降速率從135.08 mm/a變?yōu)?4 mm/a；DX的沉降速率從74.86 mm/a變?yōu)?3.82 mm/a。雖然部分特征點(diǎn)在中間時(shí)段的沉降速率有所增加，但相比于2016年，2020年6個(gè)點(diǎn)的沉降速率均減小，一定程度上反映了北京市的地面沉降程度正在減緩。

4.3 地面沉降與地下水演化相關(guān)性分析

地下水的過量開采會(huì)使含水層孔隙度逐漸變小，導(dǎo)致有效應(yīng)力增加，從而引發(fā)地面沉降。北京市由于水量短缺及水體污染導(dǎo)致可用的地表水所剩無幾，不得不長期超采地下水，從而引發(fā)大范圍的地面沉降。眾多學(xué)者結(jié)合北京市的地面沉降現(xiàn)狀，分析了其與地下水演化的相關(guān)性[8-9]。

圖10 北京市2016～2020年LOS向形變速率Fig.10 Mean deformation velocity in LOS of the Beijing from 2016 to 2020

圖11 主要沉降區(qū)域特征點(diǎn)形變速率對(duì)比Fig.11 Comparison of deformation velocity of feature points in main subsidence areas

北京市密云水庫長期擔(dān)負(fù)著城市生活及農(nóng)業(yè)生產(chǎn)用水的重要任務(wù)，僅1999～2003年密云水庫水量就萎縮了3/4，全市超過70%的用水量只能靠抽取地下水維持，因此北京平原地區(qū)的地下水位以每年1 m的速度持續(xù)下降。盡管2003年后的10 a里，北京通過各項(xiàng)節(jié)水措施使用水量下降近七成，22%的用水也已被再生水替代，但地表水稀缺的現(xiàn)實(shí)、用水量增長的趨勢(shì)難以改變，地下水位仍在逐年下降。2014年底南水北調(diào)中線工程進(jìn)入全面輸水階段，2015年北京地下水水位自1999年以來首次上升，截至2018年，地下水水位已上升了1.2 m[10]。截至2020-08底，進(jìn)入密云水庫的南水累積超過5億m3，蓄水量達(dá)到23億m3，水庫水位超過147 m，徹底扭轉(zhuǎn)了供水量入不敷出的局面。因此，由北京地下水水位與地面沉降的對(duì)應(yīng)關(guān)系可知，南水北調(diào)工程與該區(qū)域的地面沉降緩解密切相關(guān)。

5 結(jié) 語

1)利用SBAS-InSAR技術(shù)可快速獲取研究區(qū)地面沉降信息，通過與GPS觀測(cè)值對(duì)比可知，InSAR監(jiān)測(cè)結(jié)果精度較高，可為研究區(qū)域地面沉降現(xiàn)狀及時(shí)空演化特征提供參考依據(jù)。

2)北京市地面沉降雖處于持續(xù)發(fā)展階段，但沉降速率呈逐漸減緩趨勢(shì)，南水北調(diào)工程在一定程度上解決了華北平原的缺水問題，地下水的開采量也隨之減少，從而緩解了地面沉降。

3)工程沿線的地面沉降分布在空間上極不均勻，河北省東南部的地面沉降雖幾乎連接成片，但不會(huì)影響正常輸水。而天津市支線經(jīng)過2個(gè)小沉降區(qū)，輸水存在安全隱患，為保障南水北調(diào)中線工程的正常運(yùn)行，需對(duì)其進(jìn)行深入研究。