朱 寶，羅孝文，吳自銀

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

0 引言

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

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

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

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

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

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

1 方法原理

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

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

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

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

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

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

(1)

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

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

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

1.2 ARIMA方法

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

(2)

(3)

(4)

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

上述3種模型是弱平穩(wěn)的，適用于平穩(wěn)時(shí)間序列預(yù)測(cè)，若涉及到非平穩(wěn)時(shí)間序列，則需要結(jié)合差分運(yùn)算得到平穩(wěn)序列后再進(jìn)行預(yù)測(cè)，即ARIMA模型。經(jīng)過(guò)實(shí)驗(yàn)比對(duì)，ARIMA模型的預(yù)測(cè)誤差穩(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ù)會(huì)受到梯度衰減的影響，進(jìn)而造成梯度消失問(wèn)題；如果離當(dāng)前數(shù)據(jù)較近的數(shù)據(jù)記錄過(guò)多，會(huì)造成梯度爆炸問(wèn)題。

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

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

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

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

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

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

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

(10)

(11)

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

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

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

(12)

(13)

(14)

(15)

(16)

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

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ù)，故通過(guò)相對(duì)穩(wěn)定的嘉紹大橋溝通南北兩岸，以驗(yàn)證杭州灣北側(cè)的形變數(shù)據(jù)的準(zhǔn)確度。經(jīng)過(guò)與分布于研究區(qū)的水準(zhǔn)點(diǎn)測(cè)量數(shù)據(jù)對(duì)比，本文所采用的InSAR方法得到的結(jié)果精度達(dá)到毫米級(jí)，可以滿足研究需求。

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

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

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

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

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

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

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

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

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

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

圖7 3種預(yù)測(cè)方法時(shí)序預(yù)測(cè)結(jié)果對(duì)比

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

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

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

3 結(jié)論

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

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