楊翠玉，王彥兵，趙亞麗，郭亞航，黃世強(qiáng)

(1.首都師范大學(xué) 資源環(huán)境與旅游學(xué)院，北京 100048；2.首都師范大學(xué) 三維信息獲取與應(yīng)用教育部重點實驗室，北京100048；3.首都師范大學(xué) 城市環(huán)境過程與數(shù)字模擬國家重點實驗室培育基地，北京 100048)

0 引言

由于長期超量開采地下水、城市建筑群容積率增加等原因，北京地區(qū)地面沉降加劇。從20世紀(jì)60年代，北京平原區(qū)地面沉降呈快速發(fā)展，已形成昌平沙河—八仙莊、朝陽八里莊—大郊亭、朝陽—來廣營、順義—平各莊、大興榆垡—禮賢5大沉降區(qū)[1]。2004—2009年間，各沉降區(qū)的累積沉降量分別達(dá)到昌平沙河—八仙莊為272.8 mm，朝陽八里莊—大郊亭為460.5 mm，朝陽—來廣營為434.0 mm，順義—平各莊為281.6 mm，大興榆垡—禮賢為435.8 mm[2]。2012—2015年間，朝陽—來廣營沉降區(qū)、朝陽八里莊—大郊亭沉降區(qū)的沉降速率分別達(dá)到119 mm/a、115 mm/a[3]。

針對地面沉降的研究，目前主要圍繞遙感監(jiān)測技術(shù)、演化特征及成因分析等方面進(jìn)行。遙感監(jiān)測技術(shù)目前以差分干涉測量(differential interferometry synthetic aperture radar， D-InSAR)、永久散射體干涉測量(presistent scatterers for sar interferometry，PS-InSAR)、小基線干涉測量(small baseline subset，SBAS)為主要的地面沉降監(jiān)測方法[4-6]。針對北京地區(qū)的地面沉降演化與成因，大量學(xué)者進(jìn)行了研究[2，7-13]。張雯等[2]利用統(tǒng)計方法對北京5大沉降區(qū)進(jìn)行了時序演化特征分析，發(fā)現(xiàn)沉降區(qū)的沉降量均呈平穩(wěn)上升趨勢，沉降演化的空間特征不明顯。周超凡等[12]利用全局空間自相關(guān)指數(shù)和局部空間自相關(guān)指數(shù)對北京地面沉降格局特征進(jìn)行了分析，結(jié)果表明,嚴(yán)重沉降區(qū)和輕微沉降區(qū)在空間上均存在集聚現(xiàn)象，時間特性不明顯。Zuo等[13]利用標(biāo)準(zhǔn)差橢圓方法，研究發(fā)現(xiàn)北京東部地區(qū)的地面沉降由西北向東南方向發(fā)展的空間演化特征，發(fā)現(xiàn)其時間演化特征較不明顯。Li等[14]應(yīng)用區(qū)域重心法對比分析了北京平原區(qū)地面沉降演化特征與地下水、可壓縮層厚度之間的關(guān)系。

本文為進(jìn)一步挖掘北京地面沉降的演化規(guī)律，提出利用區(qū)域重心法定量描述沉降演化特征。區(qū)域重心法是一種描述地理現(xiàn)象空間集聚程度的一種方法，它能對地理現(xiàn)象的時空演化特征進(jìn)行定量描述，被廣泛應(yīng)用于人口、經(jīng)濟(jì)等領(lǐng)域的空間分布演化規(guī)律研究[15-17]。從地理學(xué)角度，地面沉降是一種地理現(xiàn)象，用區(qū)域重心法可定量描述其時空演化特征。

1 區(qū)域重心法

重心是物體在重力場中保持平衡的點。在地學(xué)領(lǐng)域，重心的位置可以通過地理坐標(biāo)和地理空間現(xiàn)象的組合計算獲取，可以用來定量描述地理現(xiàn)象的空間分布情況及演化規(guī)律[18-19]。本文采用區(qū)域重心法定量描述區(qū)域地面沉降演化特征。

假設(shè)研究區(qū)的沉降監(jiān)測點分布大致均勻，每個沉降監(jiān)測點的屬性描述包括沉降量S和地理坐標(biāo)(Xi，Yi)，則該區(qū)域的地面沉降的重心地理坐標(biāo)如式(1)所示。

(1)

式中：X、Y代表該區(qū)域地面沉降重心的地理坐標(biāo)；Xi、Yi表示第i個沉降監(jiān)測點的地理坐標(biāo)；Si是第i個沉降監(jiān)測點的沉降量。

根據(jù)公式(1)，當(dāng)研究區(qū)沉降監(jiān)測點空間分布均勻，且均勻沉降時，重心與幾何中心位置相同；當(dāng)沉降監(jiān)測點空間分布均勻，但沉降不均勻時，則重心偏離幾何中心。因此在后續(xù)的計算中，首先要對沉降監(jiān)測點進(jìn)行處理，使得其空間分布均勻。

為描述區(qū)域地面沉降重心的偏移情況，本文以重心移動的距離Lt2-t1和方位角θt2-t1定量化表達(dá)。

(2)

(3)

式中：t1、t2表示2個不同的時間；Lt2-t1表示由t1到t2時段的區(qū)域地面沉降重心遷移的距離；θt2-t1表示由t1到t2時段的區(qū)域地面沉降重心遷移的方位角[15]。

2 地面沉降數(shù)據(jù)處理及精度驗證

本文以北京來廣營沉降區(qū)為研究區(qū)域(圖1)，利用2004年1月至2010年8月的49景ENVISAT ASAR數(shù)據(jù)和2010年11月至2015年12月的40景Radarsat-2數(shù)據(jù)，基于PS-InSAR技術(shù)獲取地面沉降信息。

首先，在PS-InSAR處理過程中，由于SAR影像范圍受到軌道位置、視角、大氣狀況等條件的影響，需要對研究區(qū)的SAR影像進(jìn)行統(tǒng)一裁剪和配準(zhǔn)等處理；其次以輔影像和主影像構(gòu)成的干涉影像對為基礎(chǔ)，利用原始數(shù)據(jù)的共軛相乘獲取干涉相位，根據(jù)軌道數(shù)據(jù)和外部數(shù)字高程模型(digital elevation model， DEM)數(shù)據(jù)計算并去除地形和地平相位；然后通過振幅離差指數(shù)大于0.65且時間相干系數(shù)大于0.70篩選沉降監(jiān)測點(來廣營沉降區(qū)(128 km2)共獲取3 002個沉降監(jiān)測點)，并利用大氣延遲相位(atmospheric delay phase， APS)在空間上低頻和時間上高頻，而非線性形變在空間高頻時間低頻的特征，通過濾波的方法將大氣延遲相位去除，從而獲取來廣營地區(qū)2004年1月至2015年12月的地面沉降數(shù)據(jù)(圖1)；接著對獲取的沉降監(jiān)測點進(jìn)行均勻化處理；最后用區(qū)域重心法定量計算地面沉降重心的移動量與移動方向，分析該地區(qū)的地面沉降時空演化特征，并結(jié)合水文地質(zhì)剖面對成因進(jìn)行了分析。具體流程如圖2所示。

注：該圖基于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2019)3333號的標(biāo)準(zhǔn)地圖制作，底圖無修改。 圖1 2004—2015年北京地區(qū)年均沉降速率圖

圖2 數(shù)據(jù)處理流程圖

根據(jù)獲取的地面沉降數(shù)據(jù)初步發(fā)現(xiàn)，北京地區(qū)存在明顯的不均勻沉降，沉降嚴(yán)重的區(qū)域主要集中在東部的朝陽、通州以及順義區(qū)，最大年均沉降速率為109 mm/a；主城區(qū)、豐臺區(qū)和石景山區(qū)相對比較穩(wěn)定，年均沉降速率小于10 mm/a。

為驗證PS-InSAR方法獲取的地面沉降數(shù)據(jù)的精度，本文用2005—2013年間同一時段、相同點位的北京市水文地質(zhì)大隊以水準(zhǔn)測量方法監(jiān)測的沉降數(shù)據(jù)進(jìn)行對比分析。

首先，以水準(zhǔn)監(jiān)測點為圓心，選取半徑200 m范圍內(nèi)的所有PS點的年均沉降速率作為該點的PS-InSAR監(jiān)測值；然后計算研究區(qū)內(nèi)8個水準(zhǔn)監(jiān)測點與對應(yīng)PS-InSAR監(jiān)測值的相對誤差，并繪制相關(guān)關(guān)系圖(圖3)。結(jié)果表明，PS-InSAR監(jiān)測值與水準(zhǔn)監(jiān)測值相對誤差在10 mm/a以內(nèi)，相關(guān)系數(shù)達(dá)到0.997，2種監(jiān)測結(jié)果具有較好的一致性，反映了PS-InSAR監(jiān)測結(jié)果的可靠性。

圖3 2005—2013年InSAR監(jiān)測結(jié)果與水準(zhǔn)結(jié)果對比

3 地面沉降演化特征分析

3.1 地面沉降演化趨勢分析

在PS-InSAR監(jiān)測的來廣營地區(qū)地面沉降數(shù)據(jù)基礎(chǔ)上，本文分析了2004年1月至2015年12月該區(qū)域地面沉降的逐年演化趨勢(圖4)。通過分析可知，來廣營地區(qū)最大累積沉降量達(dá)1 278 mm，年均沉降速率為-98～-25 mm/a。來廣營地區(qū)地面沉降漏斗向東部發(fā)展，從2008年開始，在研究區(qū)的西北、東南2處形成2個沉降漏斗中心，以較高的沉降速率累積，并且東南的沉降漏斗中心沉降速率比西北的沉降漏斗中心高，存在明顯的不均勻沉降。

注：該圖基于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2019)3333號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖4 2004—2015年來廣營逐年地面沉降演化過程

3.2 區(qū)域重心法分析

根據(jù)圖4的地面沉降演化過程，可以發(fā)現(xiàn)地面沉降的發(fā)展趨勢，但難以對其進(jìn)行定量表達(dá)。本文采用區(qū)域重心法，以累積沉降量為數(shù)據(jù)源，進(jìn)一步定量描述來廣營地面沉降的時空演化特征。

首先，對沉降監(jiān)測數(shù)據(jù)進(jìn)行均勻化處理。根據(jù)式(1)可知，區(qū)域重心的位置受研究區(qū)沉降監(jiān)測點的空間分布和沉降監(jiān)測點的沉降量二者的共同影響，為保證計算的準(zhǔn)確性，研究的沉降監(jiān)測點要均勻分布。但是，由于PS-InSAR算法獲取的沉降監(jiān)測點通常分布在人工建筑物、路燈、橋梁、人工角反射器等相干性較高的硬目標(biāo)上，所以沉降監(jiān)測點的空間不均勻分布。本文通過對研究區(qū)的沉降監(jiān)測點空間插值，以100 m為間隔均勻提取樣本點，使得點位分布均勻。

其次，根據(jù)式(2)、式(3)計算2004—2015年來廣營地面沉降重心的遷移距離與方位角(表1)，繪制區(qū)域地面沉降重心遷移軌跡圖(圖5)。根據(jù)圖5和表1的結(jié)果發(fā)現(xiàn)來廣營沉降區(qū)具有如下的時空演化特征。

1)2004—2015年，來廣營地面沉降重心由西北向東南方向遷移，累積遷移1 516.43 m，遷移的平均方位角為110.56°。

表1 2004—2015年來廣營沉降重心遷移量統(tǒng)計表

2)來廣營地區(qū)地面沉降重心遷移速度呈現(xiàn)逐年遞減趨勢，2014年地面沉降重心遷移速度達(dá)極小值，2015年沉降重心遷移速度有迅速回升的趨勢，可能由于東南和西北地區(qū)地下水開采量發(fā)生變化導(dǎo)致。

3)來廣營地面沉降重心以東西向遷移為主，南北向遷移為輔，2004年1月至2015年12月期間地面沉降重心累積向東遷移1 388.17 m，東西向的遷移速度在59～200 mm/a間波動，2007年達(dá)極大值，2014年達(dá)極小值；累積向南遷移586.73 m，南北向遷移速度在1～138 mm/a間波動，2005年達(dá)極大值，2014年達(dá)極小值。

4)來廣營地面沉降重心東西向遷移速度的變化趨勢和整體遷移速度的變化趨勢保持一致，逐年遞減，南北向遷移速度先緩慢遞增，從2011年開始迅速減緩，到2014年南北向遷移速度僅為1 mm/a，2015年東西向和南北向的遷移速度都有回升趨勢。

進(jìn)一步基于區(qū)域重心法對來廣營地區(qū)不同年份季度之間地面沉降重心的遷移距離和方位角進(jìn)行了計算，發(fā)現(xiàn)其遷移量較小，且規(guī)律性不明顯。

注：該圖基于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2019)3333號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖5 2004—2015年來廣營地面沉降重心遷移軌跡圖

由于地面沉降受地質(zhì)環(huán)境條件、地下水開采強(qiáng)度、城市開發(fā)建設(shè)等因素的影響[20]。本文從中國地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所編著的《華北平原地下水可持續(xù)利用圖集》中獲取了昌平區(qū)-通州區(qū)的水文地質(zhì)剖面圖(圖6)，進(jìn)一步分析來廣營的地面沉降成因。

注：該圖基于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2019)3333號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖6 水文地質(zhì)剖面圖[21]

水文地質(zhì)剖面圖以線段AA′為剖面線，橫跨來廣營地區(qū)(圖6)。以AA′為剖面線同時繪制累積沉降量的剖面圖(圖7)。通過圖6 和圖7的對比分析可知：來廣營地區(qū)水文地質(zhì)條件基本一致，巖性以含水砂層和粉土為主，水文地質(zhì)剖面沿線累積沉降量與深層地下水水頭具有較好的對應(yīng)關(guān)系，與淺層地下水水位的關(guān)系不明顯。結(jié)合現(xiàn)有研究成果，北京地面沉降漏斗形成的主要原因是地下水超采[2，7，22]，且開采深度逐漸增加[23]。此外，結(jié)合地面沉降監(jiān)測數(shù)據(jù)和天地圖影像，發(fā)現(xiàn)來廣營地區(qū)地面建筑均以低矮村莊為主，沉降漏斗重心的變化與地物的相關(guān)性不明顯。因此，來廣營地區(qū)不均勻沉降可能是由于深層地下水開采差異導(dǎo)致。

圖7 水文地質(zhì)剖面沿線累積沉降量剖面圖

4 結(jié)束語

本文利用永久散射體干涉測量技術(shù)(PS-InSAR)獲取了來廣營地面沉降數(shù)據(jù)，通過區(qū)域重心法對來廣營地面沉降時空演化特征進(jìn)行了定量分析，得到如下結(jié)論。

1)2004年1月到2015年12月，北京來廣營地面沉降年均沉降速率為-98～-25 mm/a，空間上呈現(xiàn)出明顯的不均勻沉降，從2008年開始，來廣營地區(qū)西北、東南2處形成2個沉降漏斗，且東南沉降漏斗的沉降速率明顯高于西北沉降漏斗。

2)來廣營地面沉降重心逐年由西北向東南方向移動，累積遷移1 516.43 m，遷移的平均方位角為110.56°；其中以東西向遷移為主，南北向遷移為輔，期間累積向東遷移1 388.17 m，遷移速度在59～200 mm/a之間波動，累積向南遷移586.73 m，遷移速度在1～138 mm/a間波動。

3)來廣營地區(qū)不均勻沉降可能由于深層地下水開采差異導(dǎo)致。

來廣營地區(qū)覆蓋北京地鐵首都機(jī)場線、首都機(jī)場高速公路等多個重要軌道交通線，隨著不均勻沉降的加劇，需要對該沉降區(qū)內(nèi)軌道交通的安全運營進(jìn)行重點關(guān)注。另外，本次研究側(cè)重于從地理學(xué)的角度分析來廣營地區(qū)地面沉降年際的時空演化特征，年內(nèi)季度之間地面沉降重心遷移變化規(guī)律不明顯。此外，受限于水文地質(zhì)數(shù)據(jù)，對于地面沉降時空演化的成因后續(xù)須進(jìn)一步展開研究。