基于時序InSAR和GPS技術(shù)的北京平原區(qū)地表三維形變場特征*

雷坤超 馬鳳山 陳蓓蓓 羅 勇 崔文君 劉 賀 田 芳

(①中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院頁巖氣與地質(zhì)工程重點(diǎn)實驗室，北京 100029，中國)(②北京市水文地質(zhì)工程地質(zhì)大隊(北京市地質(zhì)環(huán)境監(jiān)測總站)，北京 100195，中國)(③首都師范大學(xué)地面沉降機(jī)理與防控教育部重點(diǎn)實驗室，北京 100048，中國)

0 引 言

許多地球物理現(xiàn)象(如構(gòu)造運(yùn)動、活動斷裂、地震、地面沉降等)所引發(fā)的地表變形，在一定程度上均表現(xiàn)出三維形變的特征。準(zhǔn)確獲取這些地球物理現(xiàn)象三維形變場信息，對于查明其內(nèi)在發(fā)生機(jī)制、運(yùn)動過程和發(fā)展趨勢，具有重要的地學(xué)數(shù)據(jù)支撐作用(陳強(qiáng)等，2012)。北京是國際上為數(shù)不多的以地下水作為主要供水水源的特大型城市之一，地下水開采量占到全市供水總量的50%～70%，地下水資源長期處于超量開采狀態(tài)(1961～2018年，地下水資源累計虧損96億立方米)，由此導(dǎo)致平原地區(qū)第四紀(jì)地層出現(xiàn)了較為嚴(yán)重的地表變形問題(雷坤超等，2019)。以往研究中，大多數(shù)學(xué)者針對超量開采地下水引發(fā)的多孔介質(zhì)垂向固結(jié)壓密變形，即地面沉降問題，利用大量觀測資料，圍繞地面沉降監(jiān)測信息獲取(張勤等，2009；雷坤超等，2014；周毅等，2016)、地面沉降演化規(guī)律(李曼等，2016；張永紅等，2016；周超凡等，2017；張靜等，2018)和成因機(jī)制(雷坤超等，2016；郭海朋等，2017；Gong et al.，2018；Chen et al.，2019)等方面開展了較為深入的研究，取得了豐富的研究成果。但對于抽取地下水引起的多孔介質(zhì)固體骨架水平運(yùn)動或變形特征問題研究較少。主要原因是由于多孔介質(zhì)水平運(yùn)動不如垂直固結(jié)壓密變形或地面沉降那樣直觀、容易監(jiān)測(王慶良等，2002)。但大量證據(jù)表明，抽水引起的含水層水平運(yùn)動和變形確實存在(Helm，1994；Burbey，1999; Burbey et al.，2006；羅躍等，2018)。

目前，對于大尺度區(qū)域地表形變的監(jiān)測，通常采用全球定位系統(tǒng)(Global Positioning System，GPS)和衛(wèi)星雷達(dá)干涉測量(Interferometric Synthetic Aperture Radar，InSAR)方法。其中：GPS技術(shù)的優(yōu)勢在于其水平方向(E、N向)的地表位移具有較高的監(jiān)測精度和靈敏度，采用較長的觀測時段，能夠探測出毫米級地表水平位移(Hu et al.，2006)。但在垂直方向的形變監(jiān)測精度稍低，比水平精度低2～3倍(楊建圖等，2006)。InSAR技術(shù)的優(yōu)勢在于其基于影像面的大范圍測量，具有覆蓋范圍廣、空間分辨率高、無接觸遙感等特點(diǎn)，特別對于大尺度地表形變，可以有效獲取空間連續(xù)的地表形變場信息(Galloway et al.，1998；Ferretti et al.，2001；Hooper，2008；王晨興等，2018)。但由于SAR傳感器側(cè)視成像的特點(diǎn)，InSAR對于地表垂向形變較為敏感，而對于南北向的水平形變探測較為困難，由此使得該技術(shù)在三維形變信息獲取及相應(yīng)地球物理機(jī)制解釋方面受到制約(陳強(qiáng)等，2012)。因此，綜合GPS與InSAR兩種技術(shù)各自的優(yōu)勢，將高精度的GPS水平位移測量與高密度的InSAR垂直形變觀測相結(jié)合，對于獲取完整的地表三維形變場信息，進(jìn)一步認(rèn)識地表變形的形成機(jī)制具有重要意義。

本文以北京平原區(qū)為研究對象，采用GPS與InSAR技術(shù)相結(jié)合，獲取平原區(qū)時序地表三維形變場信息，分析其空間分布特征與演化規(guī)律。研究內(nèi)容主要包括以下3個方面：(1)基于永久散射體干涉測量(PS-InSAR)和GPS技術(shù)獲取平原區(qū)2013～2018年時序地表垂向形變場特征，分析其演化規(guī)律。(2)利用2013～2018年GPS同步觀測數(shù)據(jù)進(jìn)行各期聯(lián)網(wǎng)解算和平差處理，獲取各期GPS點(diǎn)在ITRF2005(International Terrestrial Reference Frame 2005)框架下和歐亞參考框架下精確三維坐標(biāo)，分析不同參考框架下GPS點(diǎn)水平運(yùn)動特征。(3)將歐亞參考框架下GPS點(diǎn)水平形變速率和InSAR垂向形變速率與不同含水層系統(tǒng)等水位線進(jìn)行耦合分析，查明地下水開采對地表水平形變產(chǎn)生的影響。本文結(jié)果可為后續(xù)深入研究北京平原區(qū)地表三維形變成因機(jī)制及水-土三維耦合模型提供重要的數(shù)據(jù)基礎(chǔ)與參考。

1 研究區(qū)概況

北京平原區(qū)第四紀(jì)地層是由5大水系(永定河水系、潮白河水系、北運(yùn)河水系、大清河水系、薊運(yùn)河水系)聯(lián)合作用形成的沖洪積扇群構(gòu)成。在山前和沖洪積扇頂部，第四系厚度在20～40 m左右，為單一的砂、砂礫石層或頂部覆蓋較薄的黏性土層。在沖洪積扇中下部和沖積平原地區(qū)，沉積物厚度逐漸增大，層次增多，顆粒逐漸由粗變細(xì)，巖性過渡為砂、砂礫石與黏性土層交錯出現(xiàn)，并以黏性土為主。在部分沉積凹陷中心，第四系厚度達(dá)1000余米(雷坤超等，2016)(圖1)。

圖1 研究區(qū)位置和地質(zhì)條件概況

與第四系沉積特征相對應(yīng)，北京平原區(qū)含水層系統(tǒng)在垂向上主要分為3個含水層組：第1含水層組(潛水和淺層承壓水)為第四系全新統(tǒng)(Q4)和上更新統(tǒng)(Q3)沖(洪)積物，含水層底板埋深在25 m和80～120 m左右。第2含水層組(中深層承壓水)為第四系中更新統(tǒng)(Q2)地層，多層結(jié)構(gòu)，巖性以中粗砂為主，部分地區(qū)含礫石，底板埋深在180 m左右。第3含水層組(深層承壓水)為第四系下更新統(tǒng)(Q1)地層，多層結(jié)構(gòu)，以中粗砂和礫石為主，底界為第四系基底，部分地區(qū)底界在260～300 m左右(圖2)。

與第四系含水層組劃分相對應(yīng)，北京平原區(qū)可劃分3個主要的壓縮層組：第1壓縮層組(Q4+Q3)，底板埋深小于100 m，廣泛分布于平原區(qū)。巖性主要為全新統(tǒng)和上更新統(tǒng)沖積、沖湖積的粉土和黏性土。第2壓縮層組(Q2)，底板埋深150～180 m，主要分布在沖洪積扇中下部、沖積平原地區(qū)。巖性主要為中更新統(tǒng)沖洪積、沖湖積粉土和黏性土層。第3壓縮層組(Q1)，底板為第四系基底，部分地區(qū)在260～300 m左右，主要分布在幾大沉積凹陷地區(qū)。巖性主要為下更新統(tǒng)河湖相沉積的黏性土層(賈三滿等，2007)(圖2)。

圖2 A—A1處水文地質(zhì)剖面

2 地表形變監(jiān)測方法

2.1 時序PS-InSAR監(jiān)測方法

永久散射體干涉測量(Permanent Scatterer for SAR Interferometry，PS-InSAR)方法主要是從長時間序列的SAR影像中選取那些保持高相干的像素點(diǎn)作為研究對象，利用它們的散射特性具有很好穩(wěn)定性的特點(diǎn)，獲得可靠的相位信息。通過對相位進(jìn)行時空分析，分解各相干目標(biāo)點(diǎn)上相位的組成，主要包括高程、地表形變以及由于大氣引起的相位變化，最后得到地表的位移。該方法可以有效克服DinSAR在空間和時間方面失相干以及大氣延遲的影響，從而提高地表形變的監(jiān)測精度，達(dá)到毫米級(Ferretti et al.，2001)。主要步驟包括：SAR影像配準(zhǔn)、差分干涉處理、PS點(diǎn)選取、地表形變速率提取、投影變換與精度驗證(雷坤超等，2016)。文中采用GAMMA軟件對覆蓋北京平原區(qū)54景RadarSAT-2 雷達(dá)衛(wèi)星 SAR影像(Wide strip模式，C波段，分辨率30 m×30 m，入射角為27.8°，2013/01/16～2018/12/22)，進(jìn)行時序差分干涉處理及PS點(diǎn)提取，剔除了基線距大于300 m的干涉像對，利用30 m分辨率SRTM數(shù)據(jù)為參考DEM，進(jìn)行地形相位去除，經(jīng)地理編碼及視線向(LOS)到垂直向投影變換后，獲取了北京平原區(qū)2013～2018年時序地表垂向形變速率。

2.2 GPS監(jiān)測方法

北京市自2008年至今，每年8～10月均對平原區(qū)地表形變GPS觀測點(diǎn)進(jìn)行聯(lián)網(wǎng)觀測。結(jié)合北京房山(BJFS)、北京十三陵(BJSH)和天津薊縣(JIXN)3個基巖基準(zhǔn)點(diǎn)和平原區(qū)14個CORS(Continuously Operating Reference Station)站，采用美國麻省理工學(xué)院研制的GAMIT/GLOBK軟件包進(jìn)行聯(lián)網(wǎng)解算和平差處理。每期(年)GPS數(shù)據(jù)處理均包括61個GPS與水準(zhǔn)測量一體點(diǎn)(GPS觀測墩底部同時建設(shè)水準(zhǔn)測量點(diǎn))和14個CORS站，共計75個GPS同步觀測數(shù)據(jù)。文中利用上述GPS處理方法，獲取了平原區(qū)2013～2018年共6期(年)GPS點(diǎn)在ITRF2005框架下精確三維坐標(biāo)。同時，為了更好地研究北京區(qū)域內(nèi)部GPS點(diǎn)水平運(yùn)動特征，以3個基巖基準(zhǔn)點(diǎn)為穩(wěn)定參考基準(zhǔn)，采用NNR-NUVEL1A模型，將ITRF框架下的水平速度場扣除歐亞板塊速度，得到所有GPS點(diǎn)在歐亞參考框架下的區(qū)域內(nèi)部水平速度場。

3 平原區(qū)地表垂向形變特征

3.1 區(qū)域地面沉降特征

圖3 利用InSAR和GPS技術(shù)獲取北京平原區(qū)2013～2018年地表垂向形變特征

圖4 北京平原區(qū)2013～2018年平均沉降速率、沉降區(qū)面積和體積統(tǒng)計

通過對北京平原區(qū)時間序列地面沉降分布特征及演化規(guī)律分析發(fā)現(xiàn)，2013～2018年地面沉降總體呈減緩趨勢，但各地沉降減緩時間存在差異性。平原區(qū)北部在2016年開始減緩，至2017年和2018年沉降減緩現(xiàn)象越發(fā)明顯；平原區(qū)東部在2018年出現(xiàn)明顯減緩。出現(xiàn)這種變化趨勢的原因主要包括3個方面：(1)南水北調(diào)工程產(chǎn)生的環(huán)境正效益。2014年底，南水正式進(jìn)京，每年向北京輸水超過10億立方米，有效緩解了北京市水資源緊缺形勢，為北京市生產(chǎn)、生活供水及地表水-地下水聯(lián)合調(diào)蓄提供了重要的水源基礎(chǔ)。(2)一系列地面沉降防控措施初見成效。南水進(jìn)京后，北京市先后對密懷順應(yīng)急水源地進(jìn)行地下水壓采，利用南水在沖洪積扇頂部開展地下水回補(bǔ)，同步實施城區(qū)自備井置換等，致使地下水開采量逐年下降，由2013年的21.10億立方米下降到2018年的16.20億立方米。2016～2018年3年間，部分地區(qū)含水層組水位出現(xiàn)不同程度回升，特別是沖洪積扇頂部(平原區(qū)北部)最為明顯。(3)平原區(qū)北部與東部地區(qū)在地層巖性及結(jié)構(gòu)特征方面差異明顯。北部地區(qū)主要是單一砂礫石層，對地下水位變化響應(yīng)較為敏感，水位抬升后，地層幾乎同步回彈或沉降明顯減緩。而東部地區(qū)主要由砂、砂礫石和黏性土層組成的多層結(jié)構(gòu)區(qū)，并以黏性土層為主，層次較多，顆粒較細(xì)，這種地層結(jié)構(gòu)對地下水位變化的響應(yīng)存在較為明顯的滯后性。

3.2 地面沉降中心特征

圖5 2013～2018年主要沉降中心地面沉降速率變化特征

圖6 2013年東部4個沉降中心剖面和2013～2018年平原區(qū)沉降中心速率統(tǒng)計

3.3 累計地面沉降特征

在累計地面沉降變化方面，文中利用InSAR、GPS和水準(zhǔn)測量3種手段，獲取了各主要沉降中心2013～2018年時間序列累計沉降量曲線(圖7)。由于選取的PS點(diǎn)和GPS-水準(zhǔn)一體點(diǎn)不能完全重合，故以GPS-水準(zhǔn)一體點(diǎn)為中心，以100 m為半徑，選取鄰近PS點(diǎn)進(jìn)行統(tǒng)計分析。結(jié)果顯示，隨著時間的推移，各沉降中心累計沉降量逐漸增大。2013～2018年朝陽金盞(P1)最大累計沉降量達(dá)816.77 mm，是北京市累計地面沉降量最大的地區(qū)。

圖7 2013～2018年InSAR、GPS和水準(zhǔn)測量獲取的累計沉降量對比圖

其次是黑莊戶地區(qū)(P2)最大累計沉降量792.10 mm；朝陽三間房(P3)、通州城區(qū)(P4)、昌平八仙莊(P5)和海淀上莊(P6)最大累計沉降量分別為749.00 mm、627.00 mm、456.00 mm和437.00 mm。通過InSAR、GPS分別與水準(zhǔn)測量結(jié)果對比發(fā)現(xiàn)，2013～2018年InSAR累計沉降量與水準(zhǔn)結(jié)果偏差在4～30 mm之間；GPS累計沉降量與水準(zhǔn)結(jié)果偏差在5～50 mm之間，但3種手段獲得的累計沉降量在發(fā)展趨勢上具有很好的一致性。

3.4 InSAR和GPS垂向形變精度分析

為定量評定InSAR和GPS兩種手段垂向形變監(jiān)測精度，在InSAR和GPS獲得的2013～2018年沉降速率圖上選取40個相同點(diǎn)位同期水準(zhǔn)測量數(shù)據(jù)進(jìn)行單點(diǎn)精度檢驗(表1、圖8)。需要說明的是，GPS測量高程值屬于大地高系統(tǒng)，水準(zhǔn)測量屬于正常高系統(tǒng)，兩者之間存在高程異常。但考慮到研究區(qū)較小，比較的是垂向形變量(即高程差值)，因此兩者之間高程異?？梢院雎圆挥?。同時，InSAR獲取的是視線向(LOS)形變量，在與水準(zhǔn)測量結(jié)果比較時，需將視線向(LOS)形變量投影變換到垂直方向。經(jīng)對比分析發(fā)現(xiàn)：InSAR監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果吻合程度較高，表現(xiàn)出相同的形變特征，兩者之間的互差均在2～8 mm范圍之間，互差均方差為4.33 mm，并且具有顯著的線性相關(guān)性，相關(guān)系數(shù)(R2)達(dá)到0.983。同時，將GPS測量結(jié)果與水準(zhǔn)結(jié)果對比分析發(fā)現(xiàn)，80%以上的GPS點(diǎn)與水準(zhǔn)點(diǎn)互差不大于10 mm，個別點(diǎn)位上兩者之間互差較大，最大為15.80 mm，互差均方差為7.56 mm，相關(guān)系數(shù)為0.946。由此說明，本次InSAR和GPS測量成果均與水準(zhǔn)測量具有很好的一致性，單點(diǎn)監(jiān)測精度較高，可以較為詳細(xì)地揭示出北京平原區(qū)地表垂向形變特征。

表1 InSAR、GPS垂向形變量與水準(zhǔn)測量結(jié)果對比

圖8 利用水準(zhǔn)測量結(jié)果對InSAR和GPS垂向形變量進(jìn)行精度檢驗

4 平原區(qū)地表水平形變特征

4.1 ITRF2005參考框架下GPS水平速度場

圖9 ITRF2005參考框架下北京平原區(qū)GPS點(diǎn)水平速度場

4.2 GPS水平速度場精度分析

表2 ITRF2005框架下GPS觀測點(diǎn)水平運(yùn)動速率及精度

4.3 歐亞參考框架下GPS水平速度場

圖10 歐亞參考框架下北京平原區(qū)內(nèi)部GPS點(diǎn)水平速度場

圖11 地下水位漏斗與地表三維形變響應(yīng)關(guān)系

5 結(jié) 論

本文將高精度的GPS水平位移測量與高密度的InSAR垂直形變觀測相結(jié)合，獲取了北京平原區(qū)時序地表三維形變場信息，分析了其分布特征與演化規(guī)律。研究結(jié)果表明：

(1)北京平原區(qū)在抽水引發(fā)的第四系附加應(yīng)力場作用下，地表呈現(xiàn)出顯著的三維變形特征，以垂向變形為主，并輔以水平向位移。

(4)在歐亞參考框架下，GPS點(diǎn)水平運(yùn)動速率明顯減小，各點(diǎn)之間非一致性變化較為明顯，各點(diǎn)運(yùn)動速率的大小和方向不具備整體趨勢性活動特征。特別是幾大活動斷裂交接部位的地面沉降嚴(yán)重區(qū)，往往也是GPS點(diǎn)水平運(yùn)動速率較大的地區(qū)。GPS點(diǎn)水平運(yùn)動方向總體指向地面沉降或地下水位降落漏斗中心，或由地下水位較高的地區(qū)指向地下水位較低的地區(qū)。這主要是抽取地下水導(dǎo)致第四系含水層系統(tǒng)在水平向產(chǎn)生的變形分量引起的。