周文韜 張文君 楊元繼 馬旭東 冉茂瑩

1 西南科技大學(xué)環(huán)境與資源學(xué)院，四川省綿陽(yáng)市青龍大道中段59號(hào)，621010 2 國(guó)家遙感中心綿陽(yáng)科技城分部，四川省綿陽(yáng)市青龍大道中段59號(hào)，621010

礦山地表沉降機(jī)理復(fù)雜，通常單一預(yù)測(cè)模型預(yù)測(cè)數(shù)據(jù)的精度較低，不能反映礦區(qū)地表的真實(shí)情況。將單一預(yù)測(cè)模型進(jìn)行組合，提取各模型的有效信息來(lái)提高預(yù)測(cè)精度是當(dāng)前的發(fā)展趨勢(shì)[1-4]。差分整合移動(dòng)平均自回歸(autoregressive integrated moving average,ARIMA)模型和Holt-Winters模型都是基于時(shí)間序列變化的預(yù)測(cè)模型，適用于礦山地表沉降時(shí)序變化預(yù)測(cè)。本文基于誘導(dǎo)有序加權(quán)平均(IOWA)算子，以SBAS-InSAR監(jiān)測(cè)值作為實(shí)際值，將ARIMA模型和Holt-Winters模型進(jìn)行加權(quán)組合，以期為礦山地表沉降監(jiān)測(cè)預(yù)測(cè)提供新方法。

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

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

小基線集(small baseline subset, SBAS)技術(shù)可有效改善D-InSAR時(shí)空失相干、大氣延遲等缺陷，其原理是將傳統(tǒng)D-InSAR的監(jiān)測(cè)形變結(jié)果作為單一觀測(cè)值，并根據(jù)最小二乘法得到高精度的形變序列[5]。由于實(shí)驗(yàn)研究目標(biāo)為地表垂向沉降變化，因此本文忽略地表水平位移的影響，將LOS向形變轉(zhuǎn)至垂向形變。

1.2 基于IOWA算子的ARIMA和Holt-Winters模型的組合預(yù)測(cè)方法原理

ARIMA模型是針對(duì)非平穩(wěn)時(shí)間序列建模的常見(jiàn)模型。ARIMA(p,d,q)稱為差分整合移動(dòng)平均自回歸模型，其中AR為自回歸項(xiàng)，p為自回歸階數(shù)，MA為移動(dòng)平均項(xiàng)，q為移動(dòng)平均階數(shù)，d為差分階數(shù)。ARIMA(p,d,q)模型可表示為：

(1)

式中，L為滯后算子，d∈Z，d>0。當(dāng)差分階數(shù)d=0時(shí)，ARIMA模型等同于ARMA模型(即平穩(wěn)時(shí)間序列模型)[6]。本文研究數(shù)據(jù)均為平穩(wěn)時(shí)序數(shù)據(jù)，故d均為0。

指數(shù)平滑預(yù)測(cè)法是以具有某種指標(biāo)的本期監(jiān)測(cè)數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)為基礎(chǔ)，引入平滑系數(shù)，以求取平均數(shù)的一種時(shí)間序列預(yù)測(cè)法，其預(yù)測(cè)公式為：

(2)

Holt-Winters線性指數(shù)平滑法是在上述一次指數(shù)平滑的基礎(chǔ)上進(jìn)行二次指數(shù)平滑，其公式為：

(3)

誘導(dǎo)有序加權(quán)平均(IOWA)算子是根據(jù)有序加權(quán)平均(OWA)算子發(fā)展而來(lái)的新型加權(quán)方法，基于IOWA算子的組合預(yù)測(cè)模型可以有效克服傳統(tǒng)組合預(yù)測(cè)方法在某個(gè)時(shí)間點(diǎn)預(yù)測(cè)精度不同的缺點(diǎn)。IOWA算子原理如下：

假設(shè)存在一組監(jiān)測(cè)值xi(i=1,2,…,N),采用m種可行的單預(yù)測(cè)模型進(jìn)行預(yù)測(cè)，第t種方法在時(shí)刻i的預(yù)測(cè)值(或擬合值)為xti(t=1,2,…,m)。設(shè)L=(l1,l2,…,lm)為m種單項(xiàng)預(yù)測(cè)在組合預(yù)測(cè)中的加權(quán)系數(shù)，且滿足：

(4)

根據(jù)加權(quán)算術(shù)平均數(shù)計(jì)算公式，令

(5)

第t種方法在時(shí)刻i的預(yù)測(cè)精度pti可表示為：

(6)

利用m種方法獲得預(yù)測(cè)值及其預(yù)測(cè)精度，可構(gòu)成一個(gè)二維數(shù)組(p1i,x1i),(p2i,x2i),…,(pmi,xmi)，將pti按從大到小排序，設(shè)a-index(ti)是序號(hào)為t的預(yù)測(cè)精度的下標(biāo)，令

FL((p1i,x1i),(p2i,x2i),…,(pmi,xmi))=

(7)

則式(7)稱為由預(yù)測(cè)精度序列p1i,p2i,…,pmi所產(chǎn)生的IOWA組合預(yù)測(cè)值。

從式(7)可以看出，組合預(yù)測(cè)的賦權(quán)系數(shù)與單項(xiàng)預(yù)測(cè)方法無(wú)關(guān)，而與單項(xiàng)預(yù)測(cè)方法在各個(gè)時(shí)間點(diǎn)的預(yù)測(cè)精度大小密切相關(guān)[8]。

令ea-index(ti)=xi-xa-index(ti)，則以誤差平方和為準(zhǔn)則的基于IOWA的組合預(yù)測(cè)最優(yōu)化模型可表示為：

(8)

2 實(shí)驗(yàn)分析

2.1 研究區(qū)概況及數(shù)據(jù)來(lái)源

本文以金昌市龍首礦西二采區(qū)為研究區(qū)(以下簡(jiǎn)稱西二采)。西二采位于我國(guó)甘肅省金昌市，其沉降區(qū)面積約1.1 km2，如圖1所示，礦區(qū)地勢(shì)較平坦，平均海拔為1 500～1 800 m。礦區(qū)地表水系不發(fā)育，常年干枯，屬溫帶大陸性氣候。近年由于地下不斷開(kāi)采，導(dǎo)致地表塌陷，裂縫明顯，從而對(duì)鎳礦資源開(kāi)采產(chǎn)生極大的安全隱患。

圖1 研究區(qū)概況Fig.1 Overview of the study area

選取時(shí)間跨度為2019-04-03～2020-01-28的連續(xù)26景SAR數(shù)據(jù)C波段的升軌Sentinel-1A斜距單視復(fù)數(shù)(SLC)影像，周期為12 d，分辨率為5 m×20 m，實(shí)驗(yàn)參數(shù)詳見(jiàn)表1。影像預(yù)處理過(guò)程中引入AUX_POEORB精密定軌星歷數(shù)據(jù)，數(shù)據(jù)處理時(shí)選用分辨率為30 m的SRTM DEM數(shù)據(jù)去除地形相位。

表1 實(shí)驗(yàn)數(shù)據(jù)參數(shù)

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

實(shí)驗(yàn)選取覆蓋研究區(qū)的26景Sentinel-1A升軌SAR影像，基于SARscape5.2.1版本進(jìn)行SBAS數(shù)據(jù)處理。為獲取可靠的形變監(jiān)測(cè)結(jié)果，設(shè)置臨界基線閾值為2 %，時(shí)間基線為150，并引入30 m分辨率的SRTM DEM數(shù)據(jù)去除地形相位。剔除效果較差的干涉對(duì)后，最終選取172個(gè)干涉對(duì)進(jìn)行差分干涉處理；設(shè)置解纏相干系數(shù)閾值為0.2，并選取47個(gè)穩(wěn)定控制點(diǎn)進(jìn)行軌道精化和重去平；進(jìn)行形變速率反演和地理編碼，去除大氣相位等誤差影響，得到視線(LOS)向年平均沉降速率和時(shí)序累積沉降結(jié)果并轉(zhuǎn)至垂向形變，其中累積垂向沉降值見(jiàn)圖2。

圖2 累積沉降值Fig.2 Cumulative settlement

由圖2可以看出，礦區(qū)在監(jiān)測(cè)時(shí)段內(nèi)由地下開(kāi)采引起的最大累積垂向沉降值達(dá)到-94.5mm。為便于后文分析礦區(qū)地表沉降預(yù)測(cè)，將26景影像視為26期時(shí)序沉降變化，其中每12 d視為1期，并選取11個(gè)監(jiān)測(cè)點(diǎn)進(jìn)行時(shí)序沉降監(jiān)測(cè)分析，包含沿走向分布的6個(gè)監(jiān)測(cè)點(diǎn)(A1～A6)、沿傾向分布的4個(gè)監(jiān)測(cè)點(diǎn)(B1～B4)及交點(diǎn)O。各點(diǎn)沉降變化曲線如圖3所示，從圖中可以看出，1～4期各監(jiān)測(cè)點(diǎn)呈輕微抬升趨勢(shì)；5～12期各監(jiān)測(cè)點(diǎn)呈沉降趨勢(shì)，且各點(diǎn)累積沉降值相差較小；12～26期整體沉降速率加快，其中接近沉降中心的A4點(diǎn)累積沉降變化最大，A1點(diǎn)距離沉降中心最遠(yuǎn)，累積沉降變化最小。

圖3 各點(diǎn)沉降變化曲線Fig.3 Settlement curve of each point

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

由于缺少地面觀測(cè)點(diǎn)數(shù)據(jù)，故將SBAS-InSAR監(jiān)測(cè)數(shù)據(jù)作為實(shí)際值，選取前23期(每期時(shí)間間隔為12 d)監(jiān)測(cè)數(shù)據(jù)對(duì)后3期(2020-01-04～28)數(shù)據(jù)進(jìn)行預(yù)測(cè)，并將后3期SBAS監(jiān)測(cè)值與預(yù)測(cè)值進(jìn)行對(duì)比分析。為證明本文提出的組合預(yù)測(cè)模型在礦區(qū)沉降預(yù)測(cè)中的可行性，建立3種方案進(jìn)行分析：方案1采用ARIMA模型進(jìn)行預(yù)測(cè)，方案2采用Holt-Winters指數(shù)平滑模型進(jìn)行預(yù)測(cè)，方案3采用基于IOWA算子的ARIMA和Holt-Winters指數(shù)平滑組合模型進(jìn)行預(yù)測(cè)。

方案1和方案2采用IBM SPSS Statistics 24.0軟件進(jìn)行時(shí)間序列預(yù)測(cè)分析，得到各預(yù)測(cè)模型前23期數(shù)據(jù)擬合度(表2)、ARIMA模型預(yù)測(cè)值(表3)和Holt-Winters模型預(yù)測(cè)值(表4)。

表4 Holt-Winters模型預(yù)測(cè)值

從表2可以看出，ARIMA模型的均方根誤差RMSE平均值為3.44，略高于Holt-Winters模型，因此Holt-Winters模型的擬合度整體上優(yōu)于ARIMA模型。但從表3、4可以看出，2種預(yù)測(cè)模型的預(yù)測(cè)值與實(shí)際值在不同期數(shù)不同點(diǎn)的精度存在差異，例如A1點(diǎn)在ARIMA模型下3期預(yù)測(cè)值的相對(duì)誤差分別為19.23 %、0.14 %和18.64 %，在Holt-Winters模型下3期預(yù)測(cè)值的相對(duì)誤差分別為4.14 %、6.39 %和20.55 %，因此引入IOWA算子可克服單一預(yù)測(cè)方法在不同時(shí)間點(diǎn)預(yù)測(cè)精度不同的缺點(diǎn)。

表2 各模型擬合度

表3 ARIMA模型預(yù)測(cè)值

方案3基于IOWA算子，由式(11)計(jì)算IOWA組合預(yù)測(cè)值，以A1點(diǎn)為例：

fL[(p1,24,x1,24),(p2,24,x2,24)]=

l1xa-index(1,24)+l2xa-index(2,24)

fL[(p1,25,x1,25),(p2,25,x2,25)]=

l1xa-index(1,25)+l2xa-index(2,25)

fL[(p1,26,x1,26),(p2,26,x2,26)]=

l1xa-index(1,26)+l2xa-index(2,26)

得出：

fL[(0.808,-17.93),(0.959,-21.28)]=

-21.28l1-17.93l2

fL[(0.999,-20.99),(0.936,-22.30)]=

-20.99l1-22.30l2

fL[(0.814,-23.87),(0.794,-23.31)]=

-23.87l1-23.31l2

將其代入式(12)，整理得到最優(yōu)化模型：

2 811.6l1l2-3 225.42l1-3 098.15l2+1 793.00

同理可得到其他點(diǎn)的最優(yōu)化模型。利用編程或MATLAB最優(yōu)化工具箱可計(jì)算出A1點(diǎn)基于IOWA算子的組合預(yù)測(cè)模型的最優(yōu)權(quán)系數(shù)為：

l1=1,l2=0

則A1點(diǎn)的預(yù)測(cè)值分別為-21.28 mm、-20.99 mm和-23.87 mm。

從A1點(diǎn)組合預(yù)測(cè)模型的IOWA最優(yōu)權(quán)系數(shù)可以看出，該方法在A1點(diǎn)的組合預(yù)測(cè)是將2個(gè)單一模型中預(yù)測(cè)精度更高的值作為其組合預(yù)測(cè)值。表5為基于IOWA算子的組合預(yù)測(cè)模型的預(yù)測(cè)值。

表5 組合預(yù)測(cè)模型預(yù)測(cè)值

為驗(yàn)證上述組合預(yù)測(cè)模型的有效性，選取均方誤差MSE和平均絕對(duì)誤差MAE作為主要評(píng)價(jià)指標(biāo)(表6)。

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

從表6可以看出，基于IOWA算子的組合模型在各點(diǎn)的均方誤差MSE和平均絕對(duì)誤差MAE均比單一模型小，表明該組合模型的預(yù)測(cè)精度更高，具有一定的有效性。

3 結(jié) 語(yǔ)

以金昌市龍首礦西二采為研究區(qū)，采用26景SAR影像進(jìn)行SBAS地表沉降監(jiān)測(cè)，得到該地區(qū)26期時(shí)序沉降變化與累積沉降值。以SBAS監(jiān)測(cè)值為實(shí)際值，基于前23期監(jiān)測(cè)值分別采用ARIMA模型和Holt-Winters模型進(jìn)行擬合并預(yù)測(cè)后3期沉降值。結(jié)果表明，單一模型在不同點(diǎn)不同時(shí)期的預(yù)測(cè)精度存在差異，其中ARIMA模型在A1、A2、A3點(diǎn)第25期的預(yù)測(cè)精度更高，Holt-Winters模型則在B1～B4點(diǎn)的預(yù)測(cè)精度更高。根據(jù)各單一模型的預(yù)測(cè)值特點(diǎn)，本文基于IOWA 算子將ARIMA模型和Holt-Winters模型進(jìn)行組合，在最優(yōu)權(quán)系數(shù)賦值時(shí)按照單一模型的預(yù)測(cè)精度大小進(jìn)行賦值，組合模型可有效克服傳統(tǒng)單一賦權(quán)的缺陷，其精度優(yōu)于各單一模型。