北京東部平原區地面沉降時空演化特征及預測

2022-12-23 08:28:44于文宮輝力陳蓓蓓周超凡

自然資源遙感 2022年4期

于文，宮輝力，陳蓓蓓，周超凡

(1.首都師范大學水資源安全北京實驗室，北京 100048； 2.首都師范大學地面沉降機理與防控教育部重點實驗室，北京 100048； 3.首都師范大學城市環境過程與數字模擬國家重點實驗室培育基地，北京 100048； 4.首都師范大學京津冀平原地下水與地面沉降國家野外科學觀測研究站，北京 100048)

0 引言

地面沉降是在自然和人為因素作用下地面高程降低的地質現象。已成為城市發展中普遍存在的環境地質問題，同時也是制約社會經濟可持續發展的重要地質災害之一。我國的區域地面沉降不可小覷。據《全國地面沉降防治規劃(2011—2020年)》統計，已有50多個城市發生地面沉降，大于200 mm地面沉降的區域面積超過7萬km2，其中華北平原、長江三角洲和汾渭盆地是其中的“重災區”[1]。北京作為首都，城市發展迅速，近年來，地面沉降作為北京平原區的主要地質災害之一，影響著城市社會與經濟的健康平穩發展，其潛在的危害和經濟損失已經受到社會和政府的廣泛關注[2]。然而，傳統監測地面沉降的精密水準測量、分層標技術以及全球定位系統等技術雖然精度較高，但是設備造價高、空間分辨率低和難以獲取大范圍信息等問題給地面沉降監測工作帶來困難。20世紀80年代出現的合成孔徑雷達干涉測量技術(interferometric synthetic aperture Radar，InSAR)彌補了傳統測量的不足，該技術具有全天時、全天候、監測范圍廣和監測精度高的特點[3]，已被廣泛應用于各個領域，如地形測繪[4]、地震[5]和地面沉降[6]等。

InSAR技術是空間對地觀測技術的一次革命性飛躍，Gabriel等[7]1989年首次論證合成孔徑雷達差分干涉測量(differential interferometric SAR，D-InSAR)技術可用于探測亞厘米級的微小地表形變。但由于D-InSAR技術存在時空失相干[8]，相位解纏中參數值評估[9]以及大氣相位影響[10]等自身局限性,為了在一定程度上解決常規D-InSAR在空間、時間去相關、大氣誤差、軌道誤差和地形誤差去除中的局限性，多時相合成孔徑差分干涉測量技術(multi-temporal InSAR, MT-InSAR)陸續應運而生。其中，永久散射體(persistent scatters，PS)-InSAR 技術[11]僅有一幅主影像，其空間基線是將所有的從影像與唯一的主影像進行配準，空間基線相對較長，在時間序列上識別具有穩定散射特性的點，對目標點經過處理得到可靠的地表形變估算結果。

在地面沉降演化特征分析方面，已有相關研究，有關學者利用全局莫蘭指數和局部莫蘭指數對北京典型地面沉降區的格局特征進行了分析，結果表明，區域地面沉降全局空間自相關十分顯著，局部空間格局分為高低聚集區與非相關區3個亞區[12]。此外有研究利用標準差橢圓方法對北京東部地區演化特征進行分析，研究結合水文地質資料，分析地面沉降的空間演化模式。結果表明，直到2012年，地面沉降主要向西北-東南方向發展，然后向各個方向擴展[13]。有研究提出用區域重心法來分析來廣營的地面沉降時空演化特征，得到地面沉降遷移方向和距離[14]。現有研究在地面沉降矩陣場整體時空特性的研究較少，相關演化特征分析的方法仍存在不足，需要進一步探索新方法。

地面沉降作為一種緩變的城市地質災害會對城市基礎設施造成損毀，嚴重威脅著城市安全并造成巨大的經濟損失。建立一個高效的地面沉降預測模型對于地面沉降的防治和保障城市安全有著重要意義[15]。在地面沉降預測研究方面，主要分為物理模型和數學模型。物理模型主要通過水土耦合模型建立地面沉降預測模型[16]，這種方法一般需要大量的先驗知識和實測數據，在大型的工程項目中被廣泛應用。由于一些構建地面沉降物理模型的土體參數很難獲取，并且這種模型的計算效率較低，因此很大程度上限制了這種模型的使用范圍。數學模型是利用離散的時序歷史沉降數據的數學統計規律建立地面沉降預測模型[17]。相對于物理模型，該數學模型所需數據更容易獲取且運算高效。近年來，地理人工智能快速發展，一些智能算法被廣泛應用于地面沉降。有關研究利用灰色馬爾科夫模型(grey-Markov model，IGMM)建立了北京平原的地面沉降預測模型，發現該預測模型精度較高[18]。隨著研究的深入，深度學習越來越受關注。有相關學者分別利用長短期記憶(long short term memory，LSTM)模型[19]與循環神經網絡(recurrent neural network，RNN)模型[20]對滄州市和撫順市地面沉降進行預測。其中RNN運算中每個隱含細胞單元計算最后需要經過一個非線性函數，輸出[0,1]間的結果，這使多次運算后的數值不斷衰減且無法記憶較遠位置數據。LSTM是一種特殊的RNN，解決了RNN運算的部分局限性，其通過在隱藏層增加門機制來控制信息的流失，再經過反向傳播過程的動態調整，使得網絡可以學習長距離的時間序列數據。盡管該深度學習算法取得一定進展，但在地面沉降時序預測上仍面臨著挑戰。

本文首先利用PS-InSAR技術對北京東部平原區地面沉降信息進行獲取。鑒于當前對地面沉降要素場區域整體的研究有所不足，運用經驗正交函數(empirical orthogonal function，EOF)對地面沉降信息的空間特征與時序特征進行了系統分析與規律挖掘，并利用一種改進的LSTM模型對區域地面沉降進行預測，為城市健康發展提供依據。

1 研究區及其數據源

1.1 研究區概況

北京位于華北平原北部邊緣，地理位置在E 115° 20′～117° 33′，N39° 23′～41° 05′之間。研究區(圖1)位于北京平原東部，是地面沉降嚴重的地區，是典型的溫帶大陸性季風氣候，年均氣溫為11～12 ℃，四季分明。春季多風，降雨較少，易發生干旱；夏季降雨較多且溫度較高；秋季秋高氣爽，溫度適中，降雨適中；冬季盛行西北風，寒冷干燥，降雨較少。根據歷史資料記載，北京西單和東單一帶最早在1935年就發生了地面沉降，建國之后，隨著城市的不斷發展，地面沉降的范圍也在逐步擴大。從沉降發育歷史來看，北京平原區地面沉降先后經歷了形成階段(1955—1973年)、發展階段(1973—1983年)、擴展階段[21](1983—1999年)、快速發展階段[22](1999—2014年)、區域發展不平衡階段[23](2014—2016年)。其中朝陽、通州沉降區連成一片，成為北京平原區沉降最發育的地區。

圖1 研究區概況Fig.1 Overview of the study area

1.2 數據源

本研究中使用的數據集來自2顆不同的衛星，一個是加拿大的RADARSAT-2衛星，另一個是來自歐空局的Sentinel-1衛星。其中，本研究選取48景存檔RADARSAT-2降軌數據，時間跨度為2010年11月22日—2015年11月20日；選取61景存檔Sentinel-1升軌數據，時間跨度為2016年1月14日—2018年11月11日。數據集參數如表1所示。研究中利用Sarproz軟件來處理SAR數據集。在PS-InSAR處理過程中，使用空間分辨率為30 m的航天飛機雷達地形任務(shuttle Radar topography mission，SRTM)數字高程模型(digital elevation model，DEM)去除地形相位，并對干涉圖進行地理編碼。

表1 S1A雷達影像信息情況Tab.1 S1A Radar image information

2 研究方法

2.1 PS-InSAR

PS-InSAR彌補了差分干涉測量的不足，為監測城市微小形變帶來便利。該方法是由Ferretti等[11]在2000年提出。利用研究區域內的多景SAR影像，通過分析幅度與相位信息，找出不受時間、空間基線失相關等因素影響的穩定點目標，并將經過長時間的推移，仍能保持穩定散射特性的點稱為PS點。

該方法主要針對覆蓋同一個區域的N景SAR影像，選擇出最佳主影像，將其他的N-1景SAR影像與主影像進行配準，除去地形相位的影響，得到N-1幅差分干涉結果圖。則第i幅干涉圖的第j個點目標的相位公式為：

(1)

差分干涉處理后，進行PS點的選取，選取散射特性強且穩定的像素作為PS點，可以是具有二面角和散射特性強的建筑物，如道路邊緣、橋梁、裸露的巖石等。然后對提取的PS點相位信息進行濾波，去除大氣相位的影響，得到變形相位信息，從而得到地面沉降信息。

為了便于對從地面獲得的數據進行比較和分析，需要根據雷達成像幾何將視線方向(line of sight,LOS)獲取數據轉換為垂直數據du，公式為：

(2)

式中：dlos為LOS方向的變形。

2.2 EOF

EOF也稱為特征向量分析，是一種分析矩陣數據結構特征并提取原始數據特征量的方法，可以用來分析變量場的結構特征，包括空間模態和時間序列。Lorenz[25]在20世紀50年代首次將其引入氣象和氣候研究，現在被廣泛應用于地球科學、水文學和其他學科[26-27]。EOF分析為時空分解，既反映空間特征，也體現時間變化。即

X=EOFm×mPCm×n，

(3)

式中：m為PS點數；n為時間月數;PCm×n為主成分矩陣。具體算法為：

1)數據矩陣標準化預處理，得到一個數據矩陣Xm×n。

2)計算協方差矩陣，即

(4)

3)計算Am×m的特征根λ和特征向量Vm×m，兩者滿足

Am×mVm×m=Vm×mEm×m，

(5)

式中E是m×m維對角陣，即

(6)

一般特征根從大到小排列,每一個非零的特征根λ對應一列特征向量值，也稱EOF。如λ1對應的特征向量稱第一個EOF模態。

4)計算主成分。將EOF投影到原始矩陣，得到每個空間特征向量對應的時間系數，即

(7)

5)計算方差貢獻率。X的方差大小可以用特征根λ表示，第k個模態對總的方差貢獻率P為：

(8)

2.3 Attention-LSTM模型

LSTM在傳統 RNN 的基礎上引入3個門[28]：遺忘門、輸入門、輸出門，結構如圖2所示。

圖2 LSTM單元結構示意圖Fig.2 Schematic diagram of LSTM unit structure

具體表達式為：

ft=σ(Wf·[Ht-1,Xt]+Bf)，

(9)

it=σ(Wi·[Ht-1,Xt]+Bi)，

(10)

(11)

(12)

ot=σ(Wo·[Ht-1,Xt]+Bo)，

(13)

Ht=ot·tanhCt，

(14)

Attention機制是一種模擬人腦注意力的模型，可以根據輸入的每項特征對輸出的影響為神經網絡中的隱層狀態賦予不同的權重。本文將該機制引入地面沉降時序預測模型中，以選擇性地關注不同時間步的輸入對預測結果的影響，從而改善預測效果。Attention-LSTM的地面沉降預測框架如圖3所示。

圖3 Attention-LSTM的地面沉降預測框架Fig.3 Land subsidence prediction framework using Attention-LSTM

3 結果與分析

3.1 北京東部平原區地面沉降信息獲取

通過PS-InSAR方法對RARDARSAT-2與Sentinel-1數據處理，得到時間序列具有穩定后向散射信號的PS點。利用ArcGIS軟件中克里金插值方法對PS點在不同時間的年沉降量屬性信息進行插值處理。圖4為北京東部平原2011—2018年地面沉降年累計圖，可以看出該時段地面沉降在持續發生。截至2018年底，最大累計沉降量達1 030 mm。沉降區域主要分布在北京市朝陽—通州交界處。

圖4-1 研究區累計沉降量Fig.4-1 Cumulative settlement of the study area

圖4-2 研究區累計沉降量Fig.4-2 Cumulative settlement of the study area

圖5為年沉降圖，其中可以看出2011—2018年間年沉降量，最大沉降發生在朝陽金盞，其中2011—2015年平均沉降速率最大約為146 mm/a，2016—2018年平均沉降速率最大約為138 mm/a。通過時序年沉降圖，可以明顯看出2015年后，年沉降量有所減緩。為了對地面沉降結果進行驗證，其中以水準點為中心，利用ArcGIS軟件進行緩沖區分析，獲取所有水準點200 m緩沖范圍，提取緩沖區內的PS點，分別將各自緩沖范圍內的PS點取平均值，然后與相應的水準點進行驗證，結果如圖6所示，得到較高的相關性值，其中2015—2016年RADARSAT-2結果與相應時期水準信息的相關系數達0.978； 2017—2018年Sentinel-1結果與相應水準信息的相關系數達0.954。

圖5 研究區時序年沉降量Fig.5 Time series annual settlement map of the study area

(a) 2015—2016年 (b) 2017—2018年圖6 InSAR結果與水準結果驗證Fig.6 Verification of InSAR results and leveling results

3.2 基于EOF的北京東部平原區典型地面沉降時空演化特征分析

分別對研究區2011—2018年時間段內年尺度累計沉降變化與月尺度沉降變化進行EOF分析。由于PS點較密，為了提高運行效率，對點進行抽稀預處理。在整個年變化時序處理中，使用9 852個點，圖7顯示了年沉降變化的空間模態(EOF-1和EOF-2)以及相關的主成分分量(PC1和PC2)。空間模態，也稱空間特征向量，反映地面沉降場的空間分布特征；主成分對應的是時間變化，也稱時間系數。EOF分析的前2種模態捕捉了長期趨勢和季節變化，分別解釋了總方差的99.5%和0.4%。年沉降變化模態1的特征值為138 503 654，模態2的特征值為65 543。從圖7可以看出模態1的特征向量以朝陽金盞、管莊和通州八里橋為中心向外輻射，該中心點為高值中心，且是負值，該模態對應的PC1與地面沉降變化的長期趨勢和年際變化有關。結合振幅時間序列PC1，負高值中心與負時間系數變化疊加效應為整體增加。結果表明，東部平原地面沉降年尺度變化在2011—2014年期間不斷增大，在2015—2016年期間沉降變化減緩，高值中心即紅色渲染部分附近的強度較其他區域變化幅度要大。PC2與主要的年際周期特征有關，峰值在2015年，谷值分別在2011年和2018年。結合模式2的空間格局，可以看出該模態存在2個大紅渲染區域正高值中心，位于朝陽金盞—樓梓莊區域與通州臺湖村，1個負高值中心位于朝陽化工橋，正高值中心與負高值中心區域出現相反的年際周期變化特征。

(a) 2011—2018年模態1 (b) 2011—2018年模態2

(c) 模態1對應時間系數PC1 (d) 模態2對應的時間系數PC2圖7 2011—2018年北京東部平原典型沉降區年尺度沉降的特征向量分布及對應時間系數Fig.7 Distribution of subsidence fecture vectors and corresponding time coefficients of typical subsidence areas in the eastern plain of Beijing from 2011 to 2018

由于南水北調工程實施影響地面沉降演化[23],為了進一步探究該影響，分別對2012—2014年的10 669個PS點與2016—2018年的11 914個PS點這2個時間段的月尺度沉降累計變化進行EOF分析，圖8與圖9分別展示了2012—2014年與2016—2018年期間月尺度沉降變化的空間模式(EOF-1和EOF-2)以及相關的主成分分量(PC1和PC2)。2012—2014年期間，EOF-1與EOF-2分別解釋了總方差的98.8%和0.9%，特征值分別為25 768 800與250 655； 2016—2018年期間，EOF-1與EOF-2分別解釋了總方差的99.4%和0.2%，特征值分別為7 205 189與13 443。

(a) 2012—2014年模態1 (b) 2012—2014年模態2

(c) 模態1對應時間系數PC1 (d) 模態2對應的時間系數PC2圖8 2012—2014年北京東部平原典型沉降區月尺度沉降特征向量分布及對應時間系數Fig.8 Distribution of monthly-scale subsidence feature vectors and corresponding time coefficients of typical subsidence area in the eastern plain of Beijing from 2012 to 2014

(a) 2016—2018年模態1 (b) 2016—2018年模態2

(c) 模態1對應時間系數PC1 (d) 模態2對應的時間系數PC2圖9 2016—2018年北京東部平原典型沉降區月尺度沉降的特征向量分布及對應時間系數Fig.9 Distribution of monthly-scale subsidence feature vectors and corresponding time coefficients of typical subsidence area in the eastern plain of Beijing from 2016 to 2018

從圖8可以看出模態1的特征向量以朝陽金盞、管莊和通州八里橋為中心向外輻射，該中心點為高值中心，且是正值，該模態對應的PC1與地面沉降變化的長期趨勢保持一致。正高值中心與正值時間系數變化疊加效應為整體增加。結合振幅時間序列PC1，結果表明，東部平原地面沉降月尺度變化在2012—2013年7月期間不斷增大，在2013年7月—2015年底減小，變化幅度最大也是高值中心處。PC2與主要的地面沉降季節性特征有關。PC2的時間系數變化振幅近似于正弦波，峰值在每年7—9月，谷值在1—3月。結合模式2的空間格局，可以看出很少一部分地區出現明顯的季節性周期變化，大多數地區顯示微弱的季節性周期變化。

從圖9可以看出月尺度模態1的特征向量高值中心區域范圍較2012—2014年有所縮減，時間系數值從±8 000減小到±6 000以內。該高值中心為正值，該模態對應的PC1與地面沉降變化的長期趨勢有關。結合振幅時間序列PC1，結果表明，東部平原地面沉降月尺度變化在2016—2017年下旬不斷增大，在2017年下旬—2018年底減小，負值中心的變化模式相反。PC2與主要的季節性周期特征有關。PC2的振幅近似于正弦波，峰值在1—3月，谷值在7—9月。結合模式2的空間格局，可以看出正值高中心與負值高中心相間分布，均具有明顯的季節性周期變化。

3.3 基于Attention-LSTM模型的北京東部平原區地面沉降預測

分別利用傳統LSTM與Attention-LSTM方法對時序地面沉降進行預測。發現添加Attention網絡機制后，預測精度得到提高。通過對LSTM模型增加Attention機制，可以充分學習不同沉降時間序列中的非線性關聯，進而可以捕獲研究區域中的復雜沉降機理。

在模型訓練時，通過不斷調節訓練次數epoch來提高訓練精度，圖10顯示epoch與損失函數loss的關系，可以看出，隨著epoch增大，損失函數均方根誤差不斷減小，并趨于穩定。可以看出Attention-LSTM的最終損失函數均方根誤差優于LSTM，其中Attention-LSTM的損失函數均方根誤差小于0.01，達到很高的精度。

分別對空間上點(包括PS-InSAR獲取的真實沉降點、LSTM預測點、Attention-LSTM預測點)進行克里金插值處理，得到累計沉降結果如圖11所示。從渲染色帶分析，三者結果的趨勢大體相同，但局部存在差異，在空間分布上，通過顏色拉伸結果，大體可以看出Attention-LSTM預測點的分布與PS-InSAR方法獲取地面沉降有一定的吻合度，而LSTM預測點與其他兩者局部地區有明顯的色調差異。

(d) 2016—2018年真實值 (e) 2016—2018年LSTM預測 (f) 2016—2018年Attention-LSTM預測圖11-2 LSTM與Attention-LSTM的區域沉降預測結果Fig.11-2 Regional settlement prediction results of LSTM and Attention-LSTM

為了進一步對預測結果的值進行分析，在研究區域選擇了2條剖面，一條東西向剖面，一條南北向剖面，如圖11中所示。利用ArcGIS軟件 3D分析工具，對預測值柵格影像與真實沉降柵格影響剖面所在區域的值進行提取(圖12)，對剖面所處位置可以看出剖面位置點真實值與預測值的變化趨勢保持一致性。但同時看出在剖面上的絕大多數點，Attention-LSTM預測更靠近真實測量值。

(a) 2012—2015年東西向剖面真實值與預測值 (b) 2012—2015年南北向剖面真實值與預測值

4 結論與展望

地下水的過度開采導致北京平原區域內產生嚴重的地面沉降現象。隨著南水北調工程項目的實施，地面沉降演化特征將發生變化。本次研究選取北京平原區沉降最嚴重的區域，東部平原朝陽—通州區域，探討該區域地面沉降演化情況，利用EOF模型分析地面沉降要素場的空間分布特征與時間變化。并利用深度學習方法對該區地面沉降進行時間序列預測。研究結論如下：

1)2011—2018年，北京東部平原區地面沉降主要發生在朝陽與通州交界處，截至2018年底，最大累計沉降量達1 030 mm，2011—2015年平均沉降速率最大約為146 mm/a，2016—2018年平均沉降速率最大約為138 mm/a。將PS-InSAR方法獲取的地面沉降結果與水準結果進行驗證，得到較高的相關系數。

2)利用EOF方法對2011—2018年年尺度地面沉降場與月尺度地面沉降場的空間特性與時間變化特性進行分析，發現研究區域空間模態1方差貢獻率很大，幾乎代表研究區域空間的整體演化情況。對應時間系數線性趨勢明顯。模態2有一定的方差貢獻率，但占比很小，對應的時間系數季節性顯著。

3)利用Attention-LSTM對地面沉降進行預測，得到精度較高的預測模型，擴展了深度學習在地面沉降研究方面的應用，利用單一變量對地面沉降進行預測，避免了由于數據源不足的困境。

在以后的研究中，訓練模型時，將盡可能考慮影響地面沉降的機理要素，并將其考慮到時間序列預測中。

志謝：感謝Sarproz軟件的制造商以及歐空局提供的Sentinel-1數據。