朱 寶，羅孝文，吳自銀

(1.上海交通大學(xué)海洋學(xué)院，上海 200240；2.自然資源部海底科學(xué)重點實驗室，浙江 杭州 310012；3.自然資源部第二海洋研究所，浙江 杭州 310012；4.自然資源部海洋空間資源管理技術(shù)重點實驗室，浙江 杭州 310012；5.浙江省海洋科學(xué)院，浙江 杭州 310012)

0 引言

海陸過渡帶是陸地與海洋之間的交界地帶，其中包含了復(fù)雜的物理過程、化學(xué)過程、生物過程以及地質(zhì)過程等，因而這里成為了地球表面動態(tài)環(huán)境最復(fù)雜多變的區(qū)域之一。目前海陸過渡帶區(qū)域生活著數(shù)十億人，是世界上最繁榮發(fā)達(dá)的區(qū)域之一，其演變過程與人類的可持續(xù)發(fā)展直接相關(guān)，該區(qū)域受到全球各個沿海國家的高度關(guān)注。

地面沉降是指地下支撐物的移動導(dǎo)致地面標(biāo)高損失的一種地質(zhì)災(zāi)害[1]，它是自然因素與人為因素共同影響的地質(zhì)過程。其中自然因素包括構(gòu)造下沉、地震及火山活動、氣候變化等，人為因素包括開采地下水、油氣資源、礦產(chǎn)資源等。海陸過渡帶區(qū)域的長期地面沉降累積，會導(dǎo)致洪澇災(zāi)害的加劇、建筑物壽命的縮短以及交通網(wǎng)絡(luò)的損失等[2]，這些都會嚴(yán)重威脅到人們的生命財產(chǎn)安全。因而有效監(jiān)測地面沉降并給出相應(yīng)的預(yù)警信息十分必要。杭州灣地區(qū)為典型的海陸過渡帶，其圍墾歷史悠久，存在著地面沉降現(xiàn)象[3]，本文選擇了該區(qū)域進(jìn)行地面沉降預(yù)測研究。

地面沉降監(jiān)測方法主要有水準(zhǔn)測量、GPS(Global Positioning System, 全球定位系統(tǒng))測量以及InSAR(Interferometric Synthetic Aperture Radar, 合成孔徑干涉雷達(dá))測量等[4]。其中InSAR技術(shù)具有全天時、全天候、高精度、實時性、連續(xù)性等優(yōu)點，廣泛應(yīng)用于陸地地面沉降的監(jiān)測。例如，王超 等[5]和廖明生 等[6]分別利用該技術(shù)得到了蘇州市和上海市地面沉降的分布情況，精度均達(dá)到了毫米級。但針對海陸過渡帶復(fù)雜多變的情況，如何利用InSAR時序數(shù)據(jù)對地面沉降進(jìn)行預(yù)測的相關(guān)工作開展得還不夠。

地面沉降預(yù)測方法，主要有基于物理機(jī)制的方法、基于數(shù)理統(tǒng)計的方法以及基于機(jī)器學(xué)習(xí)的方法等[7]?；谖锢頇C(jī)制的方法需要大量的實地監(jiān)測數(shù)據(jù)，提出基于巖石、水文特性演化的物理模型，并對大量復(fù)雜參數(shù)進(jìn)行求解，并且受到時效性的嚴(yán)重影響；基于數(shù)理統(tǒng)計的方法則是根據(jù)大量歷史測量數(shù)據(jù)得到統(tǒng)計規(guī)律，并進(jìn)一步進(jìn)行分析；基于機(jī)器學(xué)習(xí)的方法是從更高的數(shù)據(jù)維度提取現(xiàn)有歷史數(shù)據(jù)的特征，進(jìn)而得到對現(xiàn)有數(shù)據(jù)的反饋，不會受限于物理參數(shù)的影響[7]。

上述預(yù)測方法，其本質(zhì)是對于時間序列的預(yù)測問題，主要的模型[8]包括：傳統(tǒng)的ARIMA模型(Auto Regressive Integrated Moving Average model，差分整合移動平均自回歸模型)[9]，機(jī)器學(xué)習(xí)的支持向量機(jī)[10-11]、貝葉斯網(wǎng)絡(luò)[12]、矩陣分解[13]、高斯過程[14]、深度學(xué)習(xí)等[15-23]。

由于海陸過渡帶動態(tài)環(huán)境復(fù)雜多變，在對InSAR數(shù)據(jù)得到的地面沉降時間序列進(jìn)行單一預(yù)測模型處理后，得到的預(yù)測效果并不理想。為此，本文提出了一種混合了ARIMA和深度學(xué)習(xí)中的LSTM(Long Short-Term Memory, 長短期記憶單元)的預(yù)測方法，即對InSAR得到的形變量時間序列與ARIMA方法得到的預(yù)測時間序列作差，然后利用LSTM對得到的時間序列差值進(jìn)行學(xué)習(xí)，并把LSTM模型預(yù)測的結(jié)果加在ARIMA預(yù)測數(shù)據(jù)上，即可得到原始InSAR形變量時間序列的預(yù)測結(jié)果。本文使用杭州灣部分地區(qū)2017—2019年的地面沉降InSAR監(jiān)測數(shù)據(jù)進(jìn)行時間序列建模來對此方法進(jìn)行驗證，結(jié)果表明，該混合預(yù)測方法相對單一預(yù)測方法的精度有顯著提升。

1 方法原理

本文所采用的方法流程如圖1所示。InSAR的SLC(Single Look Complex, 單視復(fù)數(shù))影像數(shù)據(jù)首先經(jīng)過歐洲航天局的SNAP(SentiNel Application Platform)軟件預(yù)處理，得到干涉圖；然后利用StaMPS(Stanford Method for Persistent Scatterers)得到特定點的時間序列；接著利用ARIMA方法進(jìn)行初步預(yù)測，與原始序列作差，并利用LSTM對該差值序列進(jìn)行學(xué)習(xí)，以改進(jìn)ARIMA方法的誤差項；最后在此基礎(chǔ)上進(jìn)行預(yù)測。

圖1 混合ARIMA和LSTM方法時序預(yù)測流程圖

1.1 基于StaMPS方法的InSAR數(shù)據(jù)處理

StaMPS方法是由HOOPER et al[24]發(fā)展而來的一種從合成孔徑雷達(dá)采集的時間序列中提取地面形變量的算法軟件包。相比于傳統(tǒng)的PS-InSAR(Persistent Scatterer Interferometric Synthetic Aperture Radar,永久散射體合成孔徑干涉雷達(dá))選取永久散射體點的方法，StaMPS方法的選取規(guī)則可以識別散射強度較低但穩(wěn)定性較高的目標(biāo)，可以改善在非城鎮(zhèn)地區(qū)傳統(tǒng)PS-InSAR方法選取PS(Persistent Scatterer,永久散射體)點不足的情況[25]，更加適用于海陸過渡帶區(qū)域。

干涉相位圖的相位組成如式(1)所示：

φ=W{φdef+φatm+Δφorb+Δφθ+φN}

(1)

式中：φdef代表地形變化相位，φatm代表大氣延遲相位，Δφorb代表殘余軌道誤差相位，Δφθ代表殘余視向角誤差相位，φN代表噪聲相位。所有這些相位糾纏在一起構(gòu)成了干涉相位圖的相位信息[24]。

StaMPS利用了三維解纏算法，在時間上計算每個PS點的相位差異，然后設(shè)置參考點用最小二乘法在空間上進(jìn)行解纏[26]，相比于傳統(tǒng)的方法精度更高。

解纏后的相位通過時間上的高通濾波，空間上的低通濾波，得到大氣延遲相位并將其減去，從而矯正了大氣誤差。最后經(jīng)過相位反演得到形變信息。

1.2 ARIMA方法

經(jīng)典的時間序列預(yù)測模型包括自回歸模型(Auto Regressive, AR)，移動平均模型(Moving Average, MA)，以及兩者相結(jié)合的自回歸移動平均模型(Auto Regressive Moving Average, ARMA)[8]，進(jìn)一步利用差分運算，就得到了ARIMA模型。AR(p)模型如式(2)所示，MA(q)模型如式(3)所示，ARMA(p,q)模型如式(4)所示[8]：

(2)

(3)

(4)

式中：Xt是t時刻的隨機(jī)變量，αi和βi是權(quán)重系數(shù)，p和q表示時間窗長度，εt則代表了t時刻的白噪聲。

上述3種模型是弱平穩(wěn)的，適用于平穩(wěn)時間序列預(yù)測，若涉及到非平穩(wěn)時間序列，則需要結(jié)合差分運算得到平穩(wěn)序列后再進(jìn)行預(yù)測，即ARIMA模型。經(jīng)過實驗比對，ARIMA模型的預(yù)測誤差穩(wěn)定性較好[27]，可以用于后續(xù)與InSAR數(shù)據(jù)作差，進(jìn)而使利用LSTM進(jìn)行殘差序列學(xué)習(xí)結(jié)果的收斂性得到保障。

1.3 LSTM方法

循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network, RNN)是一種能夠?qū)?dāng)前的數(shù)據(jù)與之前的數(shù)據(jù)建立聯(lián)系的深度神經(jīng)網(wǎng)絡(luò)。但由于非線性函數(shù)的存在，其中離現(xiàn)在數(shù)據(jù)較遠(yuǎn)的數(shù)據(jù)會受到梯度衰減的影響，進(jìn)而造成梯度消失問題；如果離當(dāng)前數(shù)據(jù)較近的數(shù)據(jù)記錄過多，會造成梯度爆炸問題。

為解決RNN的潛在問題，LSTM模型被提出[28]，它繼承了RNN的記憶性，并增加了門控制單元對信息進(jìn)行取舍，然后經(jīng)過反向傳播的梯度調(diào)整，從而完成模型的建立。其中以隱藏層最為重要，它由多個細(xì)胞單元組成，其細(xì)胞單元結(jié)構(gòu)如圖2所示。

圖2 LSTM隱藏層細(xì)胞結(jié)構(gòu)

每個細(xì)胞單元中含有3個控制門，分別是遺忘門f、輸入門i以及輸出門o。其中，遺忘門決定了上一時刻的細(xì)胞單元狀態(tài)ct-1有多少成分保存到ct中；輸入門決定了當(dāng)前時刻輸入的xt有多少成分保存到ct之中；輸出門決定了ct中有多少成分被傳遞到ht之中。從輸入到輸出的方向上的前向計算方法為

ft=σ(Wf·[ht-1,xt]+bf)

(5)

it=σ(Wi·[ht-1,xt]+bi)

(6)

ct=ftct-1+ittanh(Wc·[ht-1,xt]+bc)

(7)

ot=σ(Wo·[ht-1,xt]+bo)

(8)

ht=ottanh(ct)

(9)

式中：xt表示t時刻該細(xì)胞單元的輸入，ct表示t時刻該細(xì)胞單元的狀態(tài)，ht表示該細(xì)胞單元的隱藏層輸出，W和b分別為對應(yīng)的權(quán)重系數(shù)矩陣和偏置項，σ和tanh分別為sigmoid函數(shù)和雙曲正切激活函數(shù)。

LSTM的模型訓(xùn)練過程與循環(huán)神經(jīng)網(wǎng)絡(luò)一致，通過反向傳播方法計算誤差項δ，然后根據(jù)誤差項計算每個權(quán)重的梯度，進(jìn)而利用梯度優(yōu)化算法更新權(quán)重完成訓(xùn)練。

1.4 利用LSTM對ARIMA的殘差序列進(jìn)行學(xué)習(xí)的預(yù)測方法

假設(shè)利用StaMPS方法得到若干個PS點，每個PS點對應(yīng)的形變量時間序列記作{Yt}，其長度為T。

(10)

(11)

利用LSTM模型對上述標(biāo)準(zhǔn)化殘差序列進(jìn)行訓(xùn)練。設(shè)測試集的長度為Ttest，訓(xùn)練用時間序列長度為L。實踐表明，單步預(yù)測的誤差是最小的，因而對時間序列進(jìn)行標(biāo)注，使其成為監(jiān)督數(shù)據(jù)時，只需要最后一個作為標(biāo)簽即可。故訓(xùn)練用的樣本輸入為(L-1)個。

LSTM的主要參數(shù)包括訓(xùn)練用時間序列長度L，神經(jīng)網(wǎng)絡(luò)層數(shù)K，隱藏層神經(jīng)元個數(shù)N，訓(xùn)練次數(shù)epochs以及學(xué)習(xí)率η?？梢圆捎镁W(wǎng)格搜索法得到最優(yōu)參數(shù)，同時采用隨訓(xùn)練次數(shù)增加逐步衰減的學(xué)習(xí)率。

確定LSTM參數(shù)后，可以開始對其進(jìn)行訓(xùn)練。采用ADAM優(yōu)化算法[29]，設(shè)定損失函數(shù)loss為絕對誤差：

(12)

(13)

(14)

(15)

(16)

式中：RMSE表示均方根誤差，MAE表示平均絕對誤差，MAP表示平均實際絕對量，μ表示平均預(yù)測精度。

2 結(jié)果與分析

2.1 數(shù)據(jù)及預(yù)處理

選取杭州灣部分區(qū)域(30°00′—30°30′N, 120°45′—120°55′E)為研究區(qū)。由于水準(zhǔn)數(shù)據(jù)采用的是杭州灣南側(cè)的紹興市上虞區(qū)數(shù)據(jù)，故通過相對穩(wěn)定的嘉紹大橋溝通南北兩岸，以驗證杭州灣北側(cè)的形變數(shù)據(jù)的準(zhǔn)確度。經(jīng)過與分布于研究區(qū)的水準(zhǔn)點測量數(shù)據(jù)對比，本文所采用的InSAR方法得到的結(jié)果精度達(dá)到毫米級，可以滿足研究需求。

以2017年1月至2019年12月獲取的Sentinel-1A影像的上升軌道、VV極化方式的SLC影像為遙感數(shù)據(jù)源，利用PS-InSAR技術(shù)選取2018年7月16日的影像作為主影像，共得到配準(zhǔn)的主、副影像干涉圖83幅，主、副影像的時間基線以及空間垂直基線的關(guān)系如圖3所示。

圖3 主、副影像時間基線與空間垂直基線關(guān)系

利用StaMPS方法篩選，共計得到274 926個PS點，進(jìn)一步得到沿雷達(dá)視線方向的地形變化時間序列，其中每個點擁有83個形變記錄，其年平均形變量及其疊加在衛(wèi)星圖上的情況如 圖4 所示。從圖4可以看出，沿杭州灣兩岸地形整體呈現(xiàn)相對下降的趨勢，而處在更靠近陸地的位置則在整體上呈現(xiàn)出相對上升的趨勢。

圖4 研究區(qū)年平均形變量

對所有PS點的形變量數(shù)據(jù)、年平均形變量數(shù)據(jù)做統(tǒng)計分析。3年時間內(nèi)，該地區(qū)的形變量在-71.98～90.25 mm范圍內(nèi)，年平均形變量為-23.65～10.93 mm/a。圖5a代表了所有PS點的形變量數(shù)據(jù)分布情況，絕大多數(shù)形變量數(shù)據(jù)集中在-20～20 mm之間；圖5b 為年平均形變量直方圖分布，年平均形變量大多集中在-10～5 mm/a的范圍內(nèi)。

圖5 所有PS點的形變量(a)以及年平均形變量(b)分布直方圖

2.2 模型預(yù)測結(jié)果

針對地面沉降漏斗區(qū)域的監(jiān)測情況，選取兩處地面沉降較為明顯的PS點，分別利用ARIMA、LSTM以及ARIMA-LSTM三種方法進(jìn)行預(yù)測以及對比。在 圖6 中用深藍(lán)色框標(biāo)注了選取的兩個點，并分別命名為A和B。采用最大似然估計得到ARIMA的最優(yōu)參數(shù)，用網(wǎng)格搜索法得到LSTM的最優(yōu)參數(shù)。ARIMA方法的(p,d,q)參數(shù)最終選擇為(5,1,0)；LSTM方法則按照時間先后順序劃分訓(xùn)練集以及測試集，2017年和2018年的數(shù)據(jù)作為訓(xùn)練集，2019年的數(shù)據(jù)作為測試集，epochs次數(shù)為100，隱藏層神經(jīng)元個數(shù)為10個；ARIMA-LSTM方法以2019年的前2/3年數(shù)據(jù)作為訓(xùn)練集，后1/3年數(shù)據(jù)為測試集，epochs次數(shù)為50，隱藏層神經(jīng)元個數(shù)為12個。

圖6 選取的PS點A(a)和B(b)

圖7為在A、B兩點不同方法時序預(yù)測結(jié)果的對比，圖8為在A、B兩點不同方法相對誤差結(jié)果的對比，其RMSE、MAE以及μ指標(biāo)結(jié)果如表1所示。從圖7中可以看到3種方法均會出現(xiàn)預(yù)測結(jié)果滯后的情況，這源自于對過去信息的一種記憶方式。在對t時刻的值進(jìn)行預(yù)測的時候，預(yù)測值往往與t時刻的前幾個時刻的值有一定關(guān)系，當(dāng)信息數(shù)量很小的情況下，前一個時刻的實際值yt-1所對應(yīng)的權(quán)重較大，因而其預(yù)測結(jié)果會更偏向于yt-1，反映在圖中就會顯示為預(yù)測結(jié)果滯后。從圖8中可以更明顯地看出3種方法在同一個測試集中的表現(xiàn)，圖中展示了3種方法的預(yù)測相對誤差按照由小到大的順序排列而得到的結(jié)果。這里的相對誤差指的是3種方法得到的時序值，分別與InSAR時間序列值作差的絕對值除以InSAR時間序列值的絕對值而得到的數(shù)值。其中，ARIMA和LSTM方法的相對誤差大小比較接近，而ARIMA-LSTM方法的相對誤差明顯更小。

圖7 3種預(yù)測方法時序預(yù)測結(jié)果對比

圖8 3種預(yù)測方法的相對誤差結(jié)果對比

表1 3種預(yù)測方法的精度指標(biāo)對比

由表1的統(tǒng)計結(jié)果可知，本文提出的ARIMA-LSTM方法，相比于傳統(tǒng)單一的ARIMA方法或者LSTM方法，預(yù)測精度明顯提升。其中，RMSE至少減小了 2.23 mm，MAE至少減小了0.98 mm，μ至少提升了15.19%。

3 結(jié)論

本文建立了一種聯(lián)合ARIMA和LSTM的地面沉降監(jiān)測和預(yù)測方法，對InSAR形變量時間序列與ARIMA預(yù)測形變量時間序列之間的差值序列進(jìn)行預(yù)測，與傳統(tǒng)的單一LSTM或ARIMA模型預(yù)測相比，其預(yù)測精度有所提升，其中，RMSE至少減小了 2.23 mm，MAE至少減小了0.98 mm，μ至少提升了15.19%。

海陸過渡帶地面沉降的監(jiān)測和預(yù)測意義重大，其中地面沉降預(yù)測更是國內(nèi)外的難點問題，通過本文的研究，不僅發(fā)現(xiàn)了杭州灣地區(qū)的嚴(yán)重沉降區(qū)域，并且使用ARIMA-LSTM方法對歷史數(shù)據(jù)進(jìn)行了計算分析，建立了一種新的預(yù)測模型，可為海陸過渡帶地面沉降的預(yù)警提供技術(shù)借鑒。