劉青豪，張永紅，鄧 敏，吳宏安，康永輝，魏鉅杰

1. 中南大學地球科學與信息物理學院，湖南 長沙 410083； 2. 中國測繪科學研究院，北京 100830

地面沉降是指在自然和人為因素作用下，由于地殼表層土體壓縮而導致區(qū)域性地面標高降低的一種環(huán)境現(xiàn)象，又稱地面下沉或地沉[1]。由于城市化的快速發(fā)展，人類對自然資源的需求不斷增加，煤、石油、天然氣、地下水等自然資源的不斷開采使得地面沉降問題變得越發(fā)嚴重，地面沉降問題已成為世界性的地質問題。因此，開展大范圍地面沉降預測預報分析，對國土空間規(guī)劃、防災減災等具有重要的研究價值和現(xiàn)實意義。

在大范圍地表沉降預測方面,目前的研究并不豐富。文獻檢索發(fā)現(xiàn)，已有的地表沉降預測方法主要可歸納為三大類：基于物理機制的方法、基于數(shù)理統(tǒng)計的方法以及基于機器學習的方法?；谖锢頇C制的方法從沉降內部的物理演化過程出發(fā)，通過實地檢測和試驗的方式獲取包含巖性特征、水文特征在內的一系列復雜的物理參數(shù)，進而對地面沉降進行模擬與預測[2-5]，常見的模型包括采礦沉降模型、地鐵沉降模型、地下水耦合模型等。然而，由于各種參數(shù)的存在，很多情況下需要對該類方法做出嚴格的假設才能應用，而這些假設有時可能會失效，如，由于地表沉降現(xiàn)象存在復雜的成因而難以確定的趨勢性及周期性等?；跀?shù)理統(tǒng)計的方法通過分析大量歷史監(jiān)測數(shù)據(jù)的內在關系及發(fā)展規(guī)律，采用包含回歸分析、灰色預測等在內的數(shù)學模型來實現(xiàn)地面沉降的模擬。然而，由于未考慮地下巖土介質的本構關系，該類方法通常難以推廣。采用此類方法的研究主要有文獻[5—7]，能否正確選擇參數(shù)模型在很大程度上決定了該類方法預測結果的準確率?；跈C器學習的方法是在計算機的輔助下進行沉降特征的學習，不受研究區(qū)域地質及水文等復雜物理參數(shù)的限制。其原理中的回歸思想與時間序列預測有著緊密的關系[8]，經典的支持向量回歸、人工神經網絡、貝葉斯網絡、小波分析等方法在時間序列預測方面均取得了不錯的效果[9-12]。該類方法依賴歷史數(shù)據(jù)的驅動?，F(xiàn)有的機器學習方法在地表沉降預測方面尚存在兩個困難：首先，高維的數(shù)據(jù)不僅影響了機器學習的速度，而且可能包含錯誤的噪聲信息，如何從高維數(shù)據(jù)中選擇關鍵的特征進行學習直接影響著沉降預測的精度；其次，由人類活動導致的地下水位升降等相關數(shù)據(jù)與沉降觀測點的粒度難以匹配，生硬的插值必然降低訓練樣本的質量。以上兩個難題使得機器學習方法因特征選擇困難或樣本粒度的差異而難以得到較高的預測精度，并且模型的泛化能力較弱，難以對存在地質條件差異的區(qū)域進行沉降預測。因此，如何利用歷史數(shù)據(jù)建立有效的預測模型，對未來地表沉降進行模擬，是一項十分重要的工作。

綜合上述分析可以發(fā)現(xiàn)：在大范圍地表沉降預測中，現(xiàn)有沉降預測方法因模型參數(shù)難以獲取或相關數(shù)據(jù)的缺乏而難以得到可靠的預測結果[2-12]，通過分段線性逼近來實現(xiàn)沉降時間序列的非線性擬合，在本質上忽略了時間序列中的非線性關聯(lián)[12]。本文提出一種基于深度學習的地表沉降時序預測方法，并通過滄州地區(qū)2012—2016年地表沉降InSAR監(jiān)測數(shù)據(jù)進行建模驗證。結果表明，該方法具有很好的預測效果。

1 基于LSTM的地表沉降時序預測方法

近年來，深度學習作為一種新的機器學習方法受到很多行業(yè)的高度關注[13-16]。其中，循環(huán)神經網絡(recurrent neural network,RNN)能夠使得網絡當前的運行不僅跟當前的輸入數(shù)據(jù)有關，并且可以與之前的數(shù)據(jù)建立聯(lián)系。然而，每個隱含細胞單元計算的最后都要經過一個非線性函數(shù)，使其輸出為[0,1]之間的結果，這使得十幾次運算后數(shù)值會衰減到很小而無法記憶較遠位置的數(shù)據(jù)。此外，梯度爆炸也是普通RNN模型容易面臨的困難。為解決這兩種問題，文獻[16—17]提出了長短期記憶(long short-term memory,LSTM)網絡模型。該模型是一種特殊的RNN，通過在隱藏層增加門機制來控制信息的流失，再經過反向傳播過程的動態(tài)調整，使得網絡可以學習長距離的時間序列數(shù)據(jù)。本文正是利用了LSTM模型的這一優(yōu)點，提出基于LSTM模型的大范圍地表沉降時序預測方法。

1.1 長短期記憶網絡模型

隱藏層是整個LSTM網絡的核心，由多個LSTM細胞單元組成，細胞單元的結構如圖1所示。

圖1 LSTM隱藏層細胞結構Fig.1 Cell structure of LSTM hidden layer

圖1中共顯示了一層隱藏層中的3個細胞單元。xt表示樣本時間序列中第t時刻的輸入，ht表示相應細胞單元的隱含狀態(tài)輸出。每個細胞單元中依次進行數(shù)據(jù)輸入、信息遺忘、細胞狀態(tài)更新以及隱含狀態(tài)輸出，其前向計算方法可以表示為

it=σ(Wxixt+Whiht-1+Wcict-1+bi)

(1)

ft=σ(Wxfxt+Whfht-1+Wcfct-1+bf)

(2)

ct=ftct-1+ittanh(Wxcxt+Whcht-1+bc)

(3)

ot=σ(Wxoxt+Whoht-1+Wcoct+bo)

(4)

ht=ottanh(ct)

(5)

式中，i、f、c、o分別表示輸入門、遺忘門、細胞狀態(tài)、輸出門；W和b分別為對應的權重系數(shù)矩陣和偏置項；σ和tanh分別為sigmoid和雙曲正切激活函數(shù)。LSTM模型訓練過程采用的是與經典的反向傳播算法原理類似的BPTT算法[18]，主要包括4個步驟:①按照前向計算方法計算LSTM細胞的輸出值，即式(5)；②反向計算每個LSTM細胞的誤差項，包括按時間和網絡層級兩個反向傳播方向；③根據(jù)相應的誤差項，計算每個權重的梯度；④應用基于梯度的優(yōu)化算法更新權重。

1.2 基于LSTM的地表沉降時序預測模型

基于前述LSTM預測原理，結合地表沉降中單變量短時序的數(shù)據(jù)特征，構建了基于LSTM的地表沉降時序預測模型。該模型的整體框架如圖2所示，包括輸入層、隱藏層、輸出層、網絡訓練以及網絡預測5個功能模塊。輸入層負責對原始地表沉降數(shù)據(jù)進行預處理以滿足網絡輸入要求，隱藏層采用圖2表示的LSTM細胞搭建循環(huán)神經網絡，輸出層提供預測結果，網絡訓練采用ADAM優(yōu)化算法，網絡預測采用迭代的方法逐點進行預測。

通過長短期記憶網絡動態(tài)的調整，模型可以充分學習不同沉降時間序列中的非線性關聯(lián)，進而可以捕獲研究區(qū)域中的復雜沉降機理。該方法不僅可以降低預測過程對歷時數(shù)據(jù)質量的要求，而且從本質上可以提高沉降預測精度與可解釋性。

圖2 基于LSTM的地表沉降時序預測框架Fig.2 Prediction framework of land subsidence time series based on LSTM

1.2.1 數(shù)據(jù)預處理

時間序列分析法是根據(jù)序列過去的變化趨勢來預測未來的發(fā)展變化，因此事物的發(fā)展需具備一定的持續(xù)性[19]。為確定最優(yōu)的序列長度，獲取高質量的訓練樣本，對累計或差分沉降量時間序列(下稱“沉降時序”)進行預處理。首先將每個觀測點的沉降時序記為Dm={d1,d2,…,dm}，則每個沉降時序中可提取出一個長度為可控參數(shù)L的訓練樣本，記為Dtrain={dm-L,dm-L+1…,dm-1}，其中，樣本的后Y個數(shù)值用作樣本標簽，前(L-Y)個數(shù)值用作樣本輸入，且滿足約束條件2≤L

(6)

(7)

式中，m-L≤t≤m-1,t∈N。

圖3 單個沉降時序中的樣本劃分Fig.3 Sample division in single settlement time series

1.2.2 網絡訓練

網絡訓練主要以隱藏層為研究對象。從圖2可以看出，隱藏層包含(L-1)個依次連接的同構LSTM細胞單元，經過每個細胞單元的模型輸出可以表示為

Ot=LSTMforward(Xt,Ht-1,ct-1)

(8)

(9)

本文選用文獻[20]提出的適應性動量估計(adaptive moment estimation，ADAM)算法作為優(yōu)化器。該算法融合了AdaGrad[21](adaptive gradient algorithm)和RMS Prop[22](root mean square Prop)算法的優(yōu)勢，能夠對不同參數(shù)計算適應性學習率并且占用較少的存儲資源。與其他隨機優(yōu)化方法相比，ADAM算法在實際應用中整體表現(xiàn)更優(yōu)[20]。設定損失函數(shù)最小為優(yōu)化目標，給定樣本分割長度L、網絡層數(shù)K、隱藏層節(jié)點數(shù)S、學習率l以及訓練次數(shù)epochs，不斷更新網絡權重，進而得到訓練好的預測模型。

在構建上述LSTM預測模型中，涉及眾多網絡參數(shù)，其中以樣本劃分長度L、網絡層數(shù)K和每層LSTM隱藏層的特征數(shù)量S最為關鍵[17]。為了達到更好的預測效果，本文采用多層網格搜索的方法對這3個參數(shù)進行優(yōu)選。以平均預測精度最高作為參數(shù)優(yōu)選的依據(jù)，即預測樣本的沉降量預測誤差ε(Y,Ypred)最小，目標函數(shù)可以表示為

minε(Y,Ypred)=|Ypred-Y|

(10)

式中,stepL、stepK和stepS分別為對應參數(shù)的網格搜索步長；L、K和S這3個參數(shù)構成了一個三維搜索空間，可以通過多層網格搜索算法獲取最優(yōu)參數(shù)組合。搜索過程主要包括3層，從內到外分別對L、K和S進行網格搜索。固定上述參數(shù)以外的其他超參數(shù)設置，包括迭代次數(shù)、批數(shù)等，根據(jù)式(10)預設3個參數(shù)的取值范圍，其中，網格搜索范圍i與Smax需要根據(jù)預測效果由人為經驗設定，m表示沉降時間序列長度。此外，學習率l也是模型訓練過程中影響較大的參數(shù)，當學習率較高時，損失函數(shù)收斂較快，反之，損失收斂較慢。研究表明，隨著損失函數(shù)收斂情況做動態(tài)調整的學習率可以使模型得到更好的預測結果，因此，本文試驗中采用衰減學習率的方式，將學習率按照實際效果每迭代20次乘以一個0.1的衰減系數(shù)。

1.2.3 網絡預測

Ypred{y1,y2,…,yn}

(11)

式中，Ypred表示模型對不同高相干點預測結果的集合。最后，通過計算網絡訓練及網絡預測過程中輸出結果Ypred與真實樣本標簽Y的偏差，定量地給出模型訓練和預測的精度。具體表示為

(12)

(13)

(14)

式中，μ表示模型平均預測精度；MAE表示平均絕對誤差；MAP表示平均實際沉降量。

1.3 對比預測模型

為驗證本文方法與現(xiàn)有時間序列預測方法相比所具有的優(yōu)勢，將LSTM模型與傳統(tǒng)的時間序列預測模型、機器學習模型進行試驗對比，選擇模型建立時間與預測誤差作為評價指標。為降低時間成本，試驗數(shù)據(jù)采用統(tǒng)一的數(shù)據(jù)源(見后文)。

對于傳統(tǒng)的時間序列預測模型，選擇二次指數(shù)平滑(second exponential smoothing method,SES)、自回歸移動平均(autoregressive integrated moving average,ARIMA)、多元線性回歸(multiple linear regression,MLR)模型進行試驗對比。其中，二次指數(shù)平滑法是一種能夠處理含有趨勢性成分的時間序列分析方法[23]，平滑指數(shù)α的取值大小決定了序列前期實際沉降值與預期值對本期預測值的影響；自回歸移動平均法是時間序列分析中應用最廣泛的模型之一[24]，模型可以表示為ARIMA(p,d,q)，其中p、d、q分別為自回歸項數(shù)、差分次數(shù)、移動平均項數(shù)[25]。這3個參數(shù)通?？梢酝ㄟ^觀測自相關函數(shù)和偏自相關函數(shù)確定。該方法不直接考慮其他相關隨機變量的變化。此外，多元線性回歸也是一種廣泛應用于預測任務的多因素分析方法[26]。在實際問題中，一個因變量往往受到多個自變量的影響，因此，多元的回歸分析成了時間序列預測中受歡迎的方法。

對于機器學習模型，本文選擇經典的支持向量回歸(support vector regression,SVR)模型作為代表。該方法的思想是通過一個非線性核函數(shù)將多維輸入映射到更高維度的特征空間后執(zhí)行回歸運算，進而得到與輸出指標對應的非線性映射關系。在回歸模型中，不同的核函數(shù)對擬合的結果會有較大的影響。本文選擇常用的高斯徑向基函數(shù)作為非線性核函數(shù)，本文采用Matlab2016中的支持向量回歸工具箱完成該方法的預測試驗，模型中懲罰參數(shù)C與核函數(shù)參數(shù)g通過交叉驗證法求取最優(yōu)設置。

2 滄州地區(qū)地表沉降預測結果與分析

2.1 數(shù)據(jù)源

選取近年來地表沉降現(xiàn)象嚴重的滄州地區(qū)[27]為研究區(qū)域。以2017年1月至2019年12月獲取的80期Sentinel-1A影像為遙感數(shù)據(jù)源，采用自主研發(fā)的多主影像相干目標小基線(multiple master-image coherent target small-baseline interferometric SAR,MCTSB-InSAR)技術[28]開展前期數(shù)據(jù)準備工作。最終反演得到430 476個高相干點的沉降時序，其中每個觀測點位擁有80個形變記錄。為驗證InSAR測量結果的可靠性和精確性，采用滄州地區(qū)國家二等水準數(shù)據(jù)進行對比驗證。結果表明，本文數(shù)據(jù)源的反演結果精度為7.2 mm，足以支持本項研究[29]。Sentinel-1A影像時空基線參數(shù)如表1所示。對時間基線序列的平穩(wěn)性檢驗結果表明，一階差分后的時間間隔是平穩(wěn)的，因此忽略了非等時間間隔對預測誤差的影響。

表1 Sentinel-1A影像時空基線參數(shù)

2.2 試驗結果

采用上述形變反演結果，對每個觀測點的沉降時序進行了預處理，共得到容量為430 476的訓練樣本。根據(jù)實際訓練情況，本文中i與Smax分別設定為7與100時足以發(fā)現(xiàn)最優(yōu)的網絡參數(shù)。首先通過網格搜索算法分析了樣本輸入、輸出長度對預測效果的影響。圖4為不同樣本劃分方式下的預測結果，其中，網絡層數(shù)設置為3，隱藏層神經元數(shù)量設置為30。平均絕對誤差熱力圖表明：隨著輸出序列長度的增加，不同樣本輸入長度下的模型誤差均有所升高，其中，每個網格數(shù)據(jù)表示對應樣本劃分形式下多步預測的平均預測誤差。圖5進一步表示了對多步預測結果的誤差分析。篇幅所限，僅以輸入序列長度為60、預測序列長度為6示例。折線圖表明：隨著預測序列時間軸的推移，相應時刻的各項誤差均有所降低，但對于同一時刻的預測結果，相同樣本情況下單步預測的誤差更低。在Radarset 2影像數(shù)據(jù)中發(fā)現(xiàn)了相同的規(guī)律，因此，受折線圖與熱力圖中的顯著趨勢的啟發(fā)，本文在網絡構建中優(yōu)先采用單步預測方式。

圖4 樣本劃分長度與平均預測誤差的網格搜索結果(K=30,S=3)Fig.4 Grid search results of sample parameters and mean absolute error(K=30,S=3)

圖5 單步預測與多步預測中不同時刻的各項誤差表現(xiàn)(L=66，Y1=6，Y2=1)Fig.5 Error performance of different time in single step prediction and multi-step prediction(L=66，Y1=6，Y2=1)

本文進一步分析了單步預測中不同網絡參數(shù)下的預測效果。由于LSTM網絡可以學習長距離的時間序列數(shù)據(jù)，因此，盡可能取較大的樣本長度以獲取更好的預測效果(L=79)。分析不同參數(shù)設置下的結果可以發(fā)現(xiàn)，模型預測精度并非一直與模型復雜程度正相關，而是在局部存在較優(yōu)的組合，當網格搜索范圍固定時，最優(yōu)參數(shù)的分布趨向于一致。為充分挖掘樣本長度對模型預測效果的影響，本文進而用累計沉降及差分沉降兩種數(shù)據(jù)對不同樣本長度下的單步預測結果進行了分析。從圖6可知，隨著訓練樣本長度的增加，網格搜索局部最小預測誤差均逐漸降低。此外，兩種沉降數(shù)據(jù)之間的預測誤差差異反映了LSTM模型更擅長對平穩(wěn)沉降時序進行預測。值得注意的是，差分沉降數(shù)據(jù)在輸入長度大于42以后預測誤差逐漸提高。這表明并非輸入時序長度越長效果越好，很可能較遠時刻的沉降信息對當前時刻的預測存在可忽略的關聯(lián)，一味地輸入較長時序反而引入了一定的噪聲誤差。

圖6 訓練樣本長度與平均預測誤差的關系(K=3,S=40)Fig.6 Relationship between training sample length and average prediction error (K= 3,S= 40)

基于上述試驗分析，最終確定LSTM預測模型參數(shù)組合如表2所示。

表2 最優(yōu)模型參數(shù)組合搜索結果

在上述參數(shù)設置下，對滄州地區(qū)地表沉降的預測試驗取得了非常好的效果。差分(累計)形變量平均絕對誤差可控制在0.3 mm以內，完全滿足現(xiàn)有InSAR數(shù)據(jù)處理規(guī)范的精度要求[29]。

圖7從整幅影像的尺度展示了對滄州地區(qū)的沉降預測結果。結合真實形變的大尺度表達，可以發(fā)現(xiàn)本文的預測結果與真實形變高度吻合，沉降漏斗顯示清晰且大部分高相干的預測誤差在±0.5 mm以內。算法運行結果顯示，430 476個觀測點中，最大差分(累計)形變量預測誤差為8.81 mm，平均預測精度達到67.4%，相關指標的計算方法表示為式(12)。因此，本文提出的模型預測結果與實際沉降量相比偏差較小，能夠反映地表沉降隨時間變化的基本規(guī)律，預測結果具有足夠的可信度。

圖7 滄州地區(qū)差分形變的預測結果Fig.7 Prediction results of differential deformation in cangzhou area

2.3 方法對比

對于本文并列分析的其余4種預測方法，均采用網格搜索法尋找了最優(yōu)的參數(shù)設置，并詳細記錄了模型的運行時間以及預測誤差，具體參數(shù)設置及預測效果見表3。

從表3可以看出，相比前3種基于數(shù)理模型的預測方法，在預測效果方面，LSTM模型的(差分)累計沉降量平均預測誤差在0.3 mm以內，預測標準差及均方誤差也一致低于傳統(tǒng)時序預測方法，表明本文基于深度學習的預測模型精度高且具有較好的穩(wěn)健性。其中，ARIMA模型需要對每一個高相干點的沉降時序進行一次建模，在參數(shù)優(yōu)化時由于需要大量的時間成本(表3中該方法的建模時間為抽樣預估值)，并不適合應用到大范圍的地表沉降預測中。而相比作為機器學習代表的SVR模型，本文方法在建模時間上降低了89.6%，平均絕對誤差降低了67.4%，在各項評價指標上均具有更好的表現(xiàn)。分析上述試驗結果可知，基于LSTM的大區(qū)域地表沉降預測方法精度較高，時間成本較低，是一種更為滿足生產安全需求的方法。

表3 不同預測模型的試驗結果對比

2.4 空間分析

InSAR時間序列分析對于揭示宏觀的地表形變格局、分異性及其演化規(guī)律提供了重要參考，而單個高相干點的預測結果意義并不大[30]。為探索本文時序預測模型的有效性，利用ArcGIS等工具對預測結果及歷史沉降信息的空間統(tǒng)計特性進行對比分析。對于時序預測結果間存在的空間自相關、分層異質、空間格局等特性，分別采用經典的Moran’s I、Getis-Ord Gi*、q-統(tǒng)計、最鄰近指數(shù)NNI等指標來衡量。結果表明，本文方法的預測結果并未顯著改變歷史沉降的空間統(tǒng)計特性，即仍然保持了原始InSAR反演結果的形變格局，這從空間上確保了深度學習方法在地表沉降預測上的可行性。詳細的空間分析結果如表4所示。

表4 時序預測結果的空間特性變化分析

3 結 論

本文通過構建LSTM深度循環(huán)神經網絡，將深度學習引入到地表沉降預測領域，實現(xiàn)了高精度、高時效的沉降預測。試驗結果表明：

(1) 本文提出的沉降預測模型累計沉降量預測誤差在0.3 mm以內，能夠在138 s以內完成40余萬觀測點的單步預測。相比傳統(tǒng)回歸預測方法，本文方法具有精度高、用時少的顯著優(yōu)勢。

(2) 基于LSTM的沉降預測方法能夠有效保持地表形變空間格局，對空間相關性、分層異質、空間點格局等方面的分析結果驗證了本文方法的有效性與可行性。

(3) 采用單一變量預測地表沉降避免了數(shù)據(jù)源信息不充分的困難，相比其他機器學習方法，本文方法能夠在沉降驅動因素數(shù)據(jù)短缺的情況下得到更高精度的預測結果。

總的來說，本文驗證了LSTM模型在地表沉降預測領域中的適用性，擴展了深度學習技術的應用范疇。研究結果表明，基于LSTM的地表沉降預測模型得到了優(yōu)異的結果，對于生產安全防范具有重要意義?；谀壳暗墓ぷ?，后續(xù)將展開進一步研究，如考慮空間相關性對地表沉降預測的影響，得到解釋能力更強的預測結果。