基于EWT-Prophet方法的地表沉降預(yù)測

周子琪 周世健 陶 蕊

1 東華理工大學(xué)測繪工程學(xué)院，南昌市廣蘭大道418號，330013 2 南昌航空大學(xué)校長辦公室，南昌市豐和南大道696號，330063

大規(guī)模的地下礦產(chǎn)開采活動易引發(fā)地面塌陷、滑坡、泥石流等地質(zhì)災(zāi)害，破壞地質(zhì)結(jié)構(gòu)，嚴(yán)重威脅礦區(qū)周邊居民的生命和財產(chǎn)安全[1]，因此對礦區(qū)地表進(jìn)行長期的有效監(jiān)測和精準(zhǔn)預(yù)測尤為重要[2]。目前，已有學(xué)者對地表沉降預(yù)測進(jìn)行研究，并取得一定的進(jìn)展[3-4]。傳統(tǒng)的灰色類預(yù)測模型雖然使用廣泛，但適用性較差；而ARMA、ARIMA等經(jīng)典模型的預(yù)測精度不穩(wěn)定，需滾動預(yù)測才能保持預(yù)測精度；神經(jīng)網(wǎng)絡(luò)預(yù)測模型雖有較高的預(yù)測精度，但存在預(yù)測過程不穩(wěn)定、參數(shù)設(shè)置復(fù)雜、自適應(yīng)能力較差等問題。針對上述問題，本文提出將Prophet模型引入地表沉降預(yù)測研究。Prophet模型是一種自適應(yīng)分解預(yù)測、擬合插值的新模型，對原始數(shù)據(jù)有較好的自適應(yīng)性，使用靈活且無需對缺失項進(jìn)行插值，也無需任何先驗條件，擬合過程速度較快，有良好的預(yù)測精度[5]。本文在深入研究Prophet模型原理后，引入經(jīng)驗小波變換(EWT)方法，基于EWT和Prophet模型構(gòu)建一種分解-預(yù)測的自適應(yīng)預(yù)測新方法，以達(dá)到提升預(yù)測精度的目的[6]。

1 研究區(qū)概況及數(shù)據(jù)來源

本文以德興礦區(qū)為研究對象，研究區(qū)位于江西省上饒市德興境內(nèi)，屬于江南丘陵地區(qū)，地勢東南部偏高、西北部偏低，地貌整體起伏較大。實驗數(shù)據(jù)選取覆蓋德興礦區(qū)的30景C波段升軌Sentinel-1A影像數(shù)據(jù)，時間跨度為2018-12-29～2019-12-24，時間間隔為12 d，具體參數(shù)見表1。采用由美國航空航天局(NASA)提供的30 m分辨率SRTM1數(shù)據(jù)作為外部參考DEM，以消除地形誤差對實驗的影響。

表1 Sentinel-1A數(shù)據(jù)參數(shù)

2 SBAS沉降監(jiān)測及組合預(yù)測模型建立

2.1 SBAS-InSAR技術(shù)原理

SBAS-InSAR是以傳統(tǒng)D-InSAR技術(shù)為基礎(chǔ)，用來獲取工作區(qū)地表形變時間序列的技術(shù)。該技術(shù)選取N+1幅SAR影像，通過設(shè)置合適的時間基線和空間基線，生成N個小基線集差分干涉對，去除地形相位后生成差分干涉圖并進(jìn)行相位解纏，將常規(guī)D-InSAR監(jiān)測的觀測結(jié)果用作單個觀測值，采用奇異值分解(SVD)方法將每組小基線數(shù)據(jù)集進(jìn)行連接，以解決在時域上采樣過于稀疏的問題。同時結(jié)合穩(wěn)定散射體的干涉相位信息，以獲取更高的空間分辨率，最后根據(jù)最小二乘準(zhǔn)則計算高精度的沉降形變時間序列。

2.2 Prophet模型原理

Prophet預(yù)測模型[7]是一種高效分析時序信號的新模型，可以處理時序信號分析中產(chǎn)生的丟失值、異常值，對不同跨度樣本的時序信號進(jìn)行長短期預(yù)測，具備較好的預(yù)測精度、良好的自適應(yīng)能力和一定的抗差性和魯棒性。Prophet模型已經(jīng)在GNSS高程信號預(yù)測[5]、電力需求預(yù)測[8]、空氣指數(shù)預(yù)測[9]和TEC電離層預(yù)測[10]等領(lǐng)域得到廣泛應(yīng)用，并取得較好的預(yù)測效果，但幾乎沒有涉及地表沉降預(yù)測領(lǐng)域。Prophet模型的詳細(xì)原理見文獻(xiàn)[7]。

Prophet采用廣義加法模型來進(jìn)行擬合平滑和函數(shù)預(yù)測，將地表沉降時序數(shù)據(jù)自適應(yīng)分解為：

y(t)=g(t)+s(t)+εt

(1)

式中，y(t)為原始地表沉降時序數(shù)據(jù)；g(t)為沉降時序數(shù)據(jù)趨勢項，表示時間序列非線性增長(非周期項)部分的變化函數(shù)；s(t)為高程方向的震蕩周期項；εt為服從正態(tài)分布的殘差項，可表示未預(yù)測到的隨機(jī)噪聲或趨勢。SBAS-InSAR監(jiān)測數(shù)據(jù)復(fù)雜，趨勢項采用邏輯回歸函數(shù)，表示為：

(2)

式中，k為地表沉降趨勢增長率；m為位移量參數(shù)；c為趨勢上限值，隨著時間t的增加，g(t)趨于c。在地殼板塊運動中，s(t)主要包含季節(jié)性的周期震蕩和微小的趨勢變化，其擬合函數(shù)通過沉降時序數(shù)據(jù)的傅里葉級數(shù)進(jìn)行構(gòu)造：

(3)

2.3 EWT方法原理

經(jīng)驗小波變換(EWT)是一種計算量小且具備較強(qiáng)魯棒性的新變換方法，其核心思想是依據(jù)原始信號中的頻譜特征，進(jìn)行自適應(yīng)分割后構(gòu)建基于傅里葉變換的正交小波濾波器組，提取出具有緊支撐傅里葉頻譜的不同AM-FM(調(diào)幅-調(diào)頻)成分。具體過程為：

1)對信號進(jìn)行傅里葉變換，求得支撐區(qū)間(0，π)的傅里葉頻譜F(ω)；

2)對于自適應(yīng)分割傅里葉頻譜F(ω)，依據(jù)香農(nóng)定理將頻譜分解為N個頻帶與N-1個分界頻率，圖1為傅里葉頻譜分割示意圖；

圖1 傅里葉頻譜分割Fig.1 Fourier spectrum segmentation

3)根據(jù)各分界頻率構(gòu)建經(jīng)驗小波φ(ω)，確定經(jīng)驗尺度函數(shù)和經(jīng)驗小波函數(shù)；

4)對F(ω)*φ(ω)反傅里葉變換，得到不同的模態(tài)分量。

(4)

2.4 EWT-Prophet模型的沉降預(yù)測方法

圖2為EWT-Prophet組合預(yù)測模型算法的流程，圖3為EWT分解示意圖?？梢钥闯?，EWT將原始沉降數(shù)據(jù)自適應(yīng)分解為1個趨勢項信號和3個周期項信號，其中趨勢項信號可有效反映原始沉降數(shù)據(jù)的趨勢性，其余周期項信號可反映原始沉降數(shù)據(jù)的微小震蕩變化。

圖2 EWT-Prophet流程Fig.2 Flow chart of EWT-Prophet

圖3 EWT分解Fig.3 EWT decomposition

本文組合模型先對沉降時序信號進(jìn)行EWT分解，將原始時序信號分解為經(jīng)驗尺度分量F0和N個經(jīng)驗小波分量FN，再對各分量進(jìn)行Prophet預(yù)測，將得到的所有預(yù)測分量疊加重構(gòu)為最終的預(yù)測信號。該方法綜合了EWT方法自適應(yīng)能力強(qiáng)、理論性強(qiáng)和Prophet模型擬合預(yù)測效果好、預(yù)測速度快等特點。EWT-Prophet組合預(yù)測模型對地表沉降時序數(shù)據(jù)進(jìn)行預(yù)測的主要步驟為：

1)使用EWT對沉降數(shù)據(jù)進(jìn)行分解：

(5)

2)對所有有效分量進(jìn)行Prophet預(yù)測，將所有預(yù)測分量重構(gòu)為最終的預(yù)測數(shù)據(jù)，并與單一的Prophet預(yù)測方法進(jìn)行比較，驗證組合模型的可靠性。

3 實驗分析

3.1 礦區(qū)SBAS-InSAR監(jiān)測及分析

實驗選取覆蓋研究區(qū)的30景升軌Sentinel-1A數(shù)據(jù)，使用GAMMA軟件進(jìn)行SBAS-InSAR處理，選取2019-06-15的SAR影像作為公共主影像，再將輔影像配準(zhǔn)到主影像上(配準(zhǔn)精度達(dá)0.001像元)，對干涉對進(jìn)行多視處理(多視系數(shù)設(shè)置為4∶1)。實驗的時間基線和空間基線閾值分別為60 d和200 m，最終構(gòu)成112對干涉對用于時序分析，圖4為干涉對基線分布。根據(jù)干涉對時空基線分布圖對影像進(jìn)行干涉處理，利用外部DEM數(shù)據(jù)去除地形相位的影響，生成差分干涉圖，并對其進(jìn)行濾波增強(qiáng)處理，采用最小費用流方法進(jìn)行相位解纏，得到研究區(qū)相位的真實值。通過查驗每組產(chǎn)生的相干系數(shù)圖、濾波后的干涉圖和解纏圖，移除干涉質(zhì)量較差的干涉對，在具有高相干性的像素點上建立模型。通過奇異值分解方法反演估算形變速率并去除殘余地形，計算相位殘差并進(jìn)行殘差分離，最終得到研究區(qū)地表沉降速率及地表沉降分布情況。

圖4 干涉對時空基線分布Fig.4 Spatial-temporal baseline of interferometric pairs distribution

實驗結(jié)果表明，研究區(qū)有3處明顯沉降，本文將這3處沉降區(qū)域分別標(biāo)記為A、B、C，并在每個區(qū)域中分別選取3個沉降特征點作為研究對象(分別標(biāo)記為#1、#2、#3)，各特征點位置如圖5所示。圖6為各區(qū)域中特征點的沉降時序，可以發(fā)現(xiàn)，A區(qū)域沉降最為嚴(yán)重，最大年平均沉降速率為490 mm/a，其中A#3位于沉降區(qū)中心，最大累積沉降量達(dá)440 mm，A#1和A#2在沉降區(qū)邊緣，沉降量分別為184 mm和131 mm；B區(qū)域最大年平均沉降速率達(dá)390 mm/a，位于沉降中心的特征點B#3累積沉降量為358 mm，沉降區(qū)邊緣的B#1和B#2累積沉降量分別為120 mm和118 mm；C區(qū)域最大年平均沉降速率可達(dá)233 mm/a，C#1、C#2和C#3的累積沉降量分別為148 mm、98 mm和117 mm。

圖5 A、B、C區(qū)域沉降速率Fig.5 Settlement rates in regions A, B and C

圖6 A、B、C區(qū)域特征點時間序列形變Fig.6 Time series deformation of feature points in areas A, B and C

馬濤等[12]對采用SBAS-InSAR技術(shù)獲取的時序沉降數(shù)據(jù)進(jìn)行分析，并引入水準(zhǔn)監(jiān)測數(shù)據(jù)進(jìn)行對比。結(jié)果表明，監(jiān)測沉降數(shù)據(jù)的平均誤差為1.62 mm，符合《地面沉降干涉雷達(dá)數(shù)據(jù)處理技術(shù)規(guī)范》中10 mm的精度要求，表明SBAS-InSAR技術(shù)監(jiān)測沉降數(shù)據(jù)的精度具有可靠性。本文將監(jiān)測數(shù)據(jù)作為真實觀測數(shù)據(jù)進(jìn)行回溯性預(yù)測研究。

3.2 地表沉降預(yù)測及結(jié)果分析

以SBAS-InSAR技術(shù)獲得的特征點沉降時序數(shù)據(jù)作為實際值，選取前27期(時間間隔為12 d)監(jiān)測數(shù)據(jù)作為樣本數(shù)據(jù)進(jìn)行預(yù)測研究，預(yù)測后3期的沉降情況，并將預(yù)測結(jié)果與實際監(jiān)測結(jié)果進(jìn)行對比分析。為證明本文提出的組合預(yù)測模型在礦區(qū)沉降預(yù)測中的可行性，分別利用Prophet模型和EWT-Prophet組合模型進(jìn)行實驗，并對實驗結(jié)果進(jìn)行對比分析。

表2和表3為Prophet模型和EWT-Prophet組合模型的預(yù)測結(jié)果和相對誤差統(tǒng)計?？梢钥闯?，組合模型的平均相對誤差分別為1.55%、2.89%、1.57%，平均相對誤差相較于單一模型分別提升48.68%、53.24%、68.21%。由此可見，引入單一Prophet模型進(jìn)行對比不足以驗證本文組合方法的有效性，因此引入經(jīng)典ARMA預(yù)測模型進(jìn)行對比分析。

表2 Prophet模型預(yù)測值

表3 EWT-Prophet組合模型預(yù)測值

表4(單位mm)為不同方法預(yù)測結(jié)果的統(tǒng)計，表5(單位mm)為預(yù)測結(jié)果殘差絕對值統(tǒng)計?？梢钥闯觯跉埐罱^對值指標(biāo)中，EWT-Prophet組合模型預(yù)測結(jié)果的殘差絕對值最小，3期平均值分別為2.516 7 mm、5.065 6 mm、2.877 8 mm，表明本文方法優(yōu)于單一Prophet模型和ARMA模型。

表4 不同方法預(yù)測結(jié)果

表5 預(yù)測結(jié)果殘差絕對值

為分析本文預(yù)測方法的有效性，利用均方根誤差(RMSE)、平均百分比誤差(MAPE)等指標(biāo)對預(yù)測結(jié)果進(jìn)行對比分析(表6)。由表可知，EWT-Prophet組合模型的精度整體優(yōu)于單一Prophet模型和ARMA模型，平均均方根誤差和平均百分比誤差分別為3.97 mm和2.00%，相較于Prophet模型和ARMA模型，RMSE分別提升51.53%、59.03%，MAPE分別提升57.81%、64.85%。實驗結(jié)果表明，本文構(gòu)建的組合模型具有更好的預(yù)測精度，且具有有效性和適用性。

表6 預(yù)測精度比較

4 結(jié) 語

本文以江西德興礦區(qū)為研究區(qū)域，首先采用覆蓋研究區(qū)的30景影像進(jìn)行SBAS地表沉降監(jiān)測，獲取該地區(qū)30期時序沉降變化與累積沉降值。然后以SBAS監(jiān)測值為研究時序數(shù)據(jù)樣本，使用前27期監(jiān)測值作為樣本數(shù)據(jù)進(jìn)行回溯性預(yù)測研究，利用EWT-Prophet組合模型與單一Prophet模型和ARMA模型進(jìn)行預(yù)測研究。結(jié)果表明，組合模型結(jié)果整體優(yōu)于單一模型，本文方法預(yù)測精度更佳，相較于Prophet模型和ARMA模型，RMSE分別提升了51.53%、59.03%，MAPE分別提升了57.81%、64.85%，說明組合模型具有適用性和有效性。

本文將EWT與Prophet預(yù)測模型相結(jié)合并引入SBAS-InSAR沉降研究，給類似的沉降時序信號數(shù)據(jù)分析提供了一種自適應(yīng)強(qiáng)、預(yù)測效果好的組合方法，但此類分解-預(yù)測的組合方法也給預(yù)測過程帶來額外的工作量。在研究過程中，礦區(qū)的沉降變化與諸多因素有關(guān)，如何建立一種高效、自適應(yīng)強(qiáng)，并能及時反映突變情況的預(yù)測方法還有待進(jìn)一步研究。