基于殘差修正GM（1，1）模型的采空區(qū)殘余變形預(yù)測(cè)研究

薛永安，冀 哲，張文志

（1.河南理工大學(xué) 測(cè)繪與國土信息工程學(xué)院，河南 焦作 454000；2.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024；3.河南省采空區(qū)場(chǎng)地生態(tài)修復(fù)與建設(shè)技術(shù)工程研究中心，河南 焦作 454000）

地下礦產(chǎn)采出后，開采區(qū)周圍巖土體的原始應(yīng)力平衡狀態(tài)被破壞，巖土體上出現(xiàn)位移與變形[1]，對(duì)采煤沉陷區(qū)土地再利用和建（構(gòu)）筑物安全運(yùn)維帶來了安全隱患。剩余變形是采空區(qū)地表移動(dòng)變形基本穩(wěn)定后的變形，也稱為采煤沉陷區(qū)的殘余變形[2]。目前，針對(duì)殘余變形預(yù)計(jì)的研究眾多[3-9]，大多基于開采沉陷的概率積分法進(jìn)行修正[6-10]，并有殘余下沉系數(shù)的計(jì)算[11]，同時(shí)有基于非線性數(shù)學(xué)模型的預(yù)測(cè)方法[4,12]和組合預(yù)測(cè)模型[13-14]被實(shí)踐驗(yàn)證，對(duì)開展殘余變形精準(zhǔn)預(yù)測(cè)模型研究具有較好的參考價(jià)值。然而，采空區(qū)殘余沉降呈非線性特點(diǎn)，且時(shí)間跨度大、變形量小，一般很少開展長(zhǎng)時(shí)間序列化監(jiān)測(cè)，樣本數(shù)據(jù)偏少，給預(yù)測(cè)模型建立及驗(yàn)證帶來了一定的困難?；疑到y(tǒng)模型是目前解決非線性問題的重要數(shù)學(xué)方法之一[15]，因其建模時(shí)所需數(shù)據(jù)樣本較少、計(jì)算簡(jiǎn)單且預(yù)測(cè)精度較好等特點(diǎn)而被廣泛研究與應(yīng)用[16-18]，基于灰色系統(tǒng)模型開展采空區(qū)地表殘余變形預(yù)測(cè)具有方法可行性?；疑到y(tǒng)理論體系中，GM（1，1）模型屬于灰色預(yù)測(cè)的核心模型，針對(duì)GM（1，1）模型的改進(jìn)、優(yōu)化等預(yù)測(cè)應(yīng)用較多，主要從數(shù)據(jù)預(yù)處理[19]、殘差修正[20]、灰色模型背景值優(yōu)化[21]、綜合改進(jìn)[22]等角度進(jìn)行了探索，研究成果及行業(yè)應(yīng)用較為豐富，為灰色系統(tǒng)理論在殘余變形預(yù)測(cè)領(lǐng)域的深入應(yīng)用提供了理論與實(shí)踐參考?；谏鲜鲅芯楷F(xiàn)狀，從預(yù)測(cè)結(jié)果殘差修正入手，開展針對(duì)GM（1，1）模型的后改進(jìn)研究，對(duì)比一元二次多項(xiàng)式擬合修正GM（1，1）0模型和正化殘差后一元二次多項(xiàng)式擬合修正GM（1，1）1模型的預(yù)測(cè)精度，為老采空區(qū)殘余變形精準(zhǔn)預(yù)測(cè)模型構(gòu)建進(jìn)行初步探索。

1 殘差修正GM（1，1）預(yù)測(cè)模型

1.1 灰色GM（1，1）模型

灰色系統(tǒng)模型相比較其它預(yù)測(cè)模型所需建模信息量小，但預(yù)測(cè)精度較高，且越靠近預(yù)測(cè)時(shí)間段時(shí)預(yù)測(cè)精度越高。目前，依據(jù)不同目的所建立的灰色系統(tǒng)預(yù)測(cè)模型較多，GM（1，1）是其中較常使用的預(yù)測(cè)模型之一，主要針對(duì)1 個(gè)變量建模進(jìn)行一階灰色預(yù)測(cè)[23]。

設(shè)非負(fù)離散數(shù)據(jù)序列x(0)為：

式中：n 為序列長(zhǎng)度。

對(duì)x(0)進(jìn)行一次累加生成，即可得到1 個(gè)生成序列x(1)＝｛x(1)（1），x(1)（2），…，x(1)（n）｝，對(duì)此累加生成序列建立一階微分方程（GM（1，1）模型）：

式中：a、b 為灰參數(shù)；t 為時(shí)間。

求出u 后代入式（2），并做累減生成，可得到還原后的原始數(shù)據(jù)序列預(yù)測(cè)值x(0)（k）：

當(dāng)k＜n 時(shí)，x(0)（k）為模型擬合值；當(dāng)k=n 時(shí)，x(0)（k）為模型濾波值；當(dāng)k＞n 時(shí)，為x(0)（k）模型預(yù)測(cè)值[23]。

1.2 監(jiān)測(cè)數(shù)據(jù)序列插值

GM（1，1）模型是針對(duì)等間距數(shù)據(jù)序列的預(yù)測(cè)模型，然而受各種因素影響，實(shí)際獲取到的變形監(jiān)測(cè)數(shù)據(jù)序列一般都是不等間距[24]，限制了GM（1，1）模型的應(yīng)用場(chǎng)景及預(yù)測(cè)有效性。針對(duì)數(shù)據(jù)不等間距的問題，學(xué)者們從不同的角度開展了研究與模型改進(jìn)[25-27]，其中，插值法通過對(duì)原始非等間距監(jiān)測(cè)數(shù)據(jù)進(jìn)行等距處理使其等間距化，然后按照等間距數(shù)據(jù)序列進(jìn)行建模，最后將數(shù)據(jù)還原[28]。常用的插值方法包括加權(quán)平均法、三次樣條法、拉格朗日法、切比雪夫多項(xiàng)式法和牛頓法，經(jīng)過對(duì)比實(shí)驗(yàn)，綜合考慮插值精度和構(gòu)造簡(jiǎn)便性等因素，選用牛頓插值法作為等間距監(jiān)測(cè)數(shù)據(jù)序列的構(gòu)造方法。

設(shè)x0，x1，…，xn為n+1 個(gè)插值節(jié)點(diǎn)，且x≠xi（i＝0，1，2，…，n），牛頓插值多項(xiàng)式的表達(dá)式f（x）為[29]：

式中：Nn（x）為牛頓插值公式；Rn（x）為牛頓插值多項(xiàng)式余項(xiàng)。

Nn（x）與Rn（x）的表達(dá)式分別為：

1.3 殘差非線性擬合修正

任何1 個(gè)復(fù)雜的曲線都可以采用1 個(gè)多項(xiàng)式來進(jìn)行描述，也就是說多項(xiàng)式可以逼近任意函數(shù)。基于GM（1，1）預(yù)測(cè)模型所建立的殘差曲線，通過預(yù)測(cè)值與觀測(cè)值確定的殘差序列，使殘差序列的平方和最小，此時(shí)的曲線就是利用多項(xiàng)式所建立的殘差非線性擬合曲線，利用該殘差曲線修正所有的預(yù)測(cè)值，從而達(dá)到提升預(yù)測(cè)精度的目的。

針對(duì)非線性特征明顯的采空區(qū)地表殘余變形GM（1，1）模型預(yù)測(cè)值殘差序列，采用一元二次多項(xiàng)式擬合方式對(duì)殘差值進(jìn)行擬合修正，改善GM（1，1）模型的預(yù)測(cè)精度。

設(shè)GM（1，1）模型預(yù)測(cè)值殘差序列為：

式中：ε(0)（k）為觀測(cè)值x(0)（k）與GM（1，1）模型預(yù)測(cè)值x^(0)（k）之間的殘差；k=2，3，…，n。

對(duì)殘差序列的一元m 次多項(xiàng)式擬合函數(shù)如下[22]：

式中：k=1，2，…，n-1；m∈［2，n/3］，cp為一元m次多項(xiàng)式系數(shù)。

經(jīng)對(duì)比實(shí)驗(yàn)后取m=2，采用最小二乘原理通過迭代法求解cp值，將得到系數(shù)值代入式（9），則殘差修正預(yù)測(cè)值的一元二次多項(xiàng)式如下：

式中：c0、c1、c2為系數(shù)。

利用式（10）求得的ε^(0)（k）對(duì)ε(0)（k）進(jìn)行修正，修正公式如下：

1.4 正化殘差非線性擬合修正

預(yù)測(cè)值殘差序列的符號(hào)很難保持一致性，而是在正負(fù)之間跳躍，使得對(duì)符號(hào)一致性有要求的數(shù)學(xué)模型難以適用。為了提升殘差修正效果，先構(gòu)建原始?xì)埐钪档姆?hào)矩陣B，記錄殘差值的正負(fù)號(hào)。

其中：

將殘差序列進(jìn)行正化（符號(hào)一致為正）處理，構(gòu)建改進(jìn)正化殘差序列ε+(0)（k）：

再對(duì)正化后的殘差序列利用公式（10）進(jìn)行擬合，得到殘差修正值ε^+(0)（k）。

將得到的殘差修正值根據(jù)原始?xì)埐畹姆?hào)進(jìn)行正負(fù)恢復(fù)，令：

其中：

最后用ε^(1)（k）補(bǔ)償殘差ε(0)（k），預(yù)測(cè)值x^ε(0)（k）計(jì)算公式如下：

2 模型等級(jí)與預(yù)測(cè)精度評(píng)價(jià)

2.1 模型等級(jí)評(píng)價(jià)

GM（1，1）模型檢驗(yàn)通常以C、P 值判定，C、P 值的計(jì)算方法見文獻(xiàn)［23］，根據(jù)C、P 值的計(jì)算結(jié)果判定GM（1，1）模型精度等級(jí)，GM（1，1）模型精度等級(jí)評(píng)價(jià)標(biāo)準(zhǔn)見表1[23]，總體評(píng)價(jià)時(shí)一般取C 值與P 值評(píng)價(jià)結(jié)果中最高等級(jí)為模型等級(jí)。以C、P 值組合評(píng)價(jià)GM（1，1）擬合模型的精度等級(jí)。

表1 GM（1，1）模型精度等級(jí)評(píng)價(jià)標(biāo)準(zhǔn)Table 1 Accuracy grade evaluation standard of GM（1，1）model

2.2 預(yù)測(cè)精度評(píng)價(jià)

一般常用預(yù)測(cè)誤差平方和、均方誤差、平均絕對(duì)誤差、平均絕對(duì)百分比誤差和均方百分比誤差5 種指標(biāo)進(jìn)行預(yù)測(cè)模型的預(yù)測(cè)精度評(píng)價(jià)[13]，選擇均方誤差和平均絕對(duì)百分比誤差組合評(píng)價(jià)GM（1，1）模型的預(yù)測(cè)精度，具體計(jì)算方法如下[13]：

1）均方誤差（MSE）。以預(yù)測(cè)誤差值平方和取平方根的平均值MSE 作為評(píng)價(jià)指標(biāo)，計(jì)算公式如下：

式中：N 為觀測(cè)次數(shù)；xt（t=1，2，…，N）為觀測(cè)值；x^t為xt的預(yù)測(cè)值。

2）平均絕對(duì)百分比誤差（MAPE）。MAPE 以預(yù)測(cè)值相對(duì)誤差絕對(duì)值的平均值作為評(píng)價(jià)指標(biāo)，計(jì)算公式如下：

均方誤差和平均絕對(duì)百分比誤差的值越大，預(yù)測(cè)精度越低，預(yù)測(cè)模型也越差，反之亦然。

3 試?yán)?yàn)證

3.1 試驗(yàn)數(shù)據(jù)

3.1.1 試驗(yàn)1（等間距數(shù)據(jù)序列）

收集到陜西省北部某煤礦52406 工作面開采后地表沉降監(jiān)測(cè)數(shù)據(jù)集，觀測(cè)方法為三角高程測(cè)量，選取地表沉降趨于穩(wěn)定后的15 期實(shí)測(cè)累計(jì)沉降觀測(cè)值為數(shù)據(jù)源，以A09 號(hào)監(jiān)測(cè)點(diǎn)開展等間距數(shù)據(jù)序列殘差修正GM（1，1）模型的采空區(qū)地表殘余變形預(yù)測(cè)及結(jié)果驗(yàn)證。A09 號(hào)監(jiān)測(cè)點(diǎn)等間距累計(jì)下沉觀測(cè)值見表2。

表2 A09 號(hào)監(jiān)測(cè)點(diǎn)等間距累計(jì)下沉觀測(cè)值Table 2 Accumulated observation values of A09 with equal time distance

A09 號(hào)點(diǎn)共有15 期近似等間距沉降觀測(cè)數(shù)據(jù)，選取該數(shù)據(jù)中的前11 期數(shù)據(jù)進(jìn)行GM（1，1）模型擬合，對(duì)后4 期數(shù)據(jù)進(jìn)行預(yù)測(cè)，并對(duì)比觀測(cè)值分析預(yù)測(cè)精度。

3.1.2 試驗(yàn)2（非等間距數(shù)據(jù)序列）

收集到山西省北部某大型輸水線路下伏老采空區(qū)地表沉降監(jiān)測(cè)數(shù)據(jù)集，觀測(cè)方法為二等水準(zhǔn)測(cè)量，選取19 期實(shí)測(cè)累計(jì)沉降觀測(cè)值為數(shù)據(jù)源，以ⅡR1—1 號(hào)監(jiān)測(cè)點(diǎn)開展非等間距數(shù)據(jù)序列殘差修正GM（1，1）模型的采空區(qū)地表殘余變形預(yù)測(cè)及結(jié)果的驗(yàn)證。ⅡR1-1 號(hào)監(jiān)測(cè)點(diǎn)非等間距累計(jì)下沉觀測(cè)值見表3。

表3 ⅡR1—1 號(hào)監(jiān)測(cè)點(diǎn)非等間距累計(jì)下沉觀測(cè)值Table 3 Accumulated observation values of ⅡR1—1 with non- equal time distance

ⅡR1—1 號(hào)點(diǎn)共有19 期非等間距沉降觀測(cè)數(shù)據(jù)，首先對(duì)非等間距監(jiān)測(cè)數(shù)據(jù)按月進(jìn)行牛頓插值法內(nèi)插，獲取到等間距監(jiān)測(cè)數(shù)據(jù)共53 期，然后選取該數(shù)據(jù)序列中的前47 期數(shù)據(jù)進(jìn)行GM（1，1）模型擬合，對(duì)后6 期數(shù)據(jù)進(jìn)行預(yù)測(cè)，并對(duì)比內(nèi)插等間距監(jiān)測(cè)數(shù)據(jù)分析預(yù)測(cè)精度。

3.2 殘差修正GM（1，1）模型預(yù)測(cè)結(jié)果

記傳統(tǒng)GM（1，1）預(yù)測(cè)模型為GM（1，1），采用非線性改正進(jìn)行GM（1，1）模型預(yù)測(cè)結(jié)果殘差修正的模型為GM（1，1）0，采用正化殘差非線性改正進(jìn)行GM（1，1）模型預(yù)測(cè)結(jié)果殘差修正的模型為GM（1，1）1。

分別采用GM（1，1）、GM（1，1）0和GM（1，1）1模型對(duì)試驗(yàn)1 與試驗(yàn)2，數(shù)據(jù)進(jìn)行預(yù)測(cè)3 種模型在A09、ⅡR1—1 號(hào)監(jiān)測(cè)點(diǎn)的擬合與預(yù)測(cè)結(jié)果對(duì)比見表4，3 種模型預(yù)測(cè)殘差對(duì)比如圖1。

圖1 3 種模型預(yù)測(cè)值殘差對(duì)比圖Fig.1 Residuals comparison of predicted values of the three models

表4 3 種模型預(yù)測(cè)結(jié)果對(duì)比表Table 4 Comparison of prediction results of three models

由表4 可知：GM（1，1）模型在A09 點(diǎn)與ⅡR1—1 點(diǎn)的擬合平均相對(duì)誤差分別為0.07%和44.86%，GM（1，1）0模型分別為0.06%和42.58%，GM（1，1）1模型分別為0.03%和23.26%；GM（1，1）、GM（1，1）0、GM（1，1）1模型在A09 點(diǎn)的擬合值殘差中誤差分別為1.62、1.36、0.71，而在ⅡR1—1 點(diǎn)的預(yù)測(cè)值殘差中誤差分別為0.81、0.80、0.42。顯然，通過對(duì)GM（1，1）模型的擬合結(jié)果進(jìn)行殘差修正的GM（1，1）0模型一定程度上提升了模型精度，而正化殘差修正的GM（1，1）1模型精度則較GM（1，1）0模型進(jìn)一步提升。

由表4 可知：GM（1，1）模型在A09 點(diǎn)與ⅡR1—1 點(diǎn)的預(yù)測(cè)平均相對(duì)誤差分別為0.20%和15.38%，GM（1，1）0模型分別為0.53%和37.99%，GM（1，1）1模型分別為0.16%和6.80%；GM（1，1）、GM（1，1）0、GM（1，1）1模型在A09 點(diǎn)的預(yù)測(cè)值殘差中誤差分別為4.95、11.95、3.85，而在ⅡR1—1 點(diǎn)的預(yù)測(cè)值殘差中誤差分別為0.32、0.75、0.15。顯然，試驗(yàn)1 與試驗(yàn)2 中，GM（1，1）1模型的預(yù)測(cè)平均相對(duì)誤差和殘差中誤差總體較GM（1，1）0和GM（1，1）模型都小，預(yù)測(cè)效果改善明顯。

圖1 顯示：試驗(yàn)1 與試驗(yàn)2 中3 種模型的預(yù)測(cè)值殘差趨勢(shì)一致，GM（1，1）1模型的預(yù)測(cè)值殘差均最小，GM（1，1）模型次之，GM（1，1）0模型的預(yù)測(cè)值殘差均最大。綜合表4 中擬合平均相對(duì)誤差、擬合殘差中誤差、預(yù)測(cè)平均相對(duì)誤差和預(yù)測(cè)殘差中誤差，可以看出，采用一元二次多項(xiàng)式擬合方式對(duì)殘差值進(jìn)行擬合修正可以提升GM（1，1）模型擬合的精度，但是模型的預(yù)測(cè)精度沒有提升，且出現(xiàn)了大幅降低，而采用正化殘差非線性擬合修正則顯著提升了GM（1，1）模型在擬合與預(yù)測(cè)中的精度。

3.3 模型檢驗(yàn)與預(yù)測(cè)精度

計(jì)算得到各監(jiān)測(cè)點(diǎn)在3 種預(yù)測(cè)模型中擬合值的C、P 值和預(yù)測(cè)值的MSE、MAPE 值，模型等級(jí)與預(yù)測(cè)精度評(píng)價(jià)見表5。

表5 模型等級(jí)與預(yù)測(cè)精度評(píng)價(jià)Table 5 Evaluation of model grade and prediction accuracy

由表5 可知：GM（1，1）模型和GM（1，1）0模型在A09 點(diǎn)、ⅡR1—1 點(diǎn)的模型等級(jí)均相同，分別為二級(jí)、四級(jí)，而GM（1，1）1模型在A09 點(diǎn)、ⅡR1—1點(diǎn)的模型等級(jí)分別為一級(jí)和二級(jí)；對(duì)比C、P 值可以看出，試驗(yàn)1 中，模型等級(jí)由高到低的順序依次是GM（1，1）1模型＞GM（1，1）0模型＞GM（1，1）模型；試驗(yàn)2 中，GM（1，1）1模型的模型等級(jí)最高（二級(jí)），GM（1，1）模型和GM（1，1）0模型則基本一致，均為四級(jí)（不合格），通過正化殘差非線性擬合修正實(shí)現(xiàn)了GM（1，1）模型的合格化，且模型等級(jí)提升為二級(jí)（合格）。

由MSE、MAPE 值可以看出，GM（1，1）1模型在A09 點(diǎn)和ⅡR1—1 點(diǎn)的預(yù)測(cè)精度顯著高于GM（1，1）0和GM（1，1）模型。

3.4 討 論

1）由表4 可以看出，試驗(yàn)1 與試驗(yàn)2 中，GM（1，1）0、GM（1，1）1模型的擬合精度均較GM（1，1）模型有所提升，但GM（1，1）1模型擬合精度增益更為明顯。其中，GM（1，1）1模型在A09 點(diǎn)擬合結(jié)果殘差的中誤差值較GM（1，1）0、GM（1，1）模型分別降低了48.25%、56.73%，而在ⅡR1—1 點(diǎn)則分別降低了51.85%、52.44%，降幅均達(dá)50%左右，對(duì)提升GM（1，1）模型預(yù)測(cè)能力提供了基礎(chǔ)保障。顯然，GM（1，1）1模型相比其它2 種模型更適用于煤礦采空區(qū)地表殘余變形預(yù)測(cè)。

2）由表5 可以看出，GM（1，1）0模型較GM（1，1）模型預(yù)測(cè)精度提升不明顯，但GM（1，1）1模型較GM（1，1）模型預(yù)測(cè)精度則提升較大。其中，GM（1，1）1模型在A09 點(diǎn)預(yù)測(cè)精度的MSE 值較GM（1，1）0、GM（1，1）模型分別降低了48.72%、57.45%，MAPE 值分別降低了60.00%、66.67%。而在ⅡR1—1 點(diǎn)預(yù)測(cè)精度的MSE 值分別降低了50.00%、50.00%，MAPE 值分別降低了45.38%、48.17%。由此表明，對(duì)擬合結(jié)果殘差先進(jìn)行正化處理后再進(jìn)行非線性擬合修正比直接進(jìn)行非線性擬合修正后預(yù)測(cè)的效果要好。

3）殘差修正模型大多對(duì)殘差序列有符號(hào)要求，目前，對(duì)GM（1，1）模型的殘差修正多采用GM（1，1）模型進(jìn)行，這種修正方式顯然未考慮殘差符號(hào)各異性特征。本研究采用非線性多項(xiàng)式擬合作為殘差修正模型，對(duì)GM（1，1）模型預(yù)測(cè)結(jié)果進(jìn)行正化殘差非線性修正，預(yù)測(cè)精度顯著提高，且該方法對(duì)殘差序列無特殊要求，適用范圍更廣，為變形預(yù)測(cè)模型后改進(jìn)研究提供了參考。

4 結(jié) 語

1）以殘差修正方式對(duì)GM（1，1）模型進(jìn)行后改進(jìn)處理，試?yán)?yàn)證的平均相對(duì)誤差、殘差中誤差、C、P、MSE 和MAPE 值統(tǒng)計(jì)結(jié)果一致表明，殘差修正對(duì)提升GM（1，1）模型預(yù)測(cè)精度具有明顯的增益，是后改進(jìn)GM（1，1）變形預(yù)測(cè)模型的有效方式。

2）采用正化殘差多項(xiàng)式擬合修正的GM（1，1）1模型對(duì)提升GM（1，1）模型預(yù)測(cè)精度具有顯著作用，且該模型不受殘差符號(hào)統(tǒng)一性的限制，可直接針對(duì)預(yù)測(cè)值殘差序列進(jìn)行擬合修正，預(yù)測(cè)精度進(jìn)一步提升，模型適用范圍更廣，較GM（1，1）模型和直接采用多項(xiàng)式擬合修正殘差的GM（1，1）0模型更適用于煤礦采空區(qū)地表殘余變形預(yù)測(cè)。

3）C、P、MSE 和MAPE 值的評(píng)價(jià)結(jié)果一致，模型等級(jí)高的預(yù)測(cè)模型其預(yù)測(cè)精度高，反之亦然，這4 種指標(biāo)互為補(bǔ)充與支撐，是評(píng)價(jià)灰色系統(tǒng)模型預(yù)測(cè)能力的合理評(píng)價(jià)指標(biāo)組合。

本研究的工作對(duì)老采空區(qū)殘余變形預(yù)測(cè)模型研究具有一定的參考價(jià)值，但本研究?jī)H采用了后改進(jìn)處理方式，未對(duì)觀測(cè)數(shù)據(jù)序列進(jìn)行前改進(jìn)處理，如卡爾曼濾波、粗差探測(cè)與剔除等，今后應(yīng)就此問題繼續(xù)開展研究，推進(jìn)綜合“前+后”改進(jìn)的GM（1，1）模型在采空區(qū)地表殘余變形預(yù)測(cè)中的應(yīng)用，推動(dòng)非線性系統(tǒng)理論變形預(yù)測(cè)方法體系發(fā)展。