成 樞 朱玉明 馬衛(wèi)驕 高秀明 劉忠偉

（1.山東科技大學(xué)測繪科學(xué)與工程學(xué)院，山東青島266590；2.山東省核工業(yè)二四八地質(zhì)大隊，山東青島266041；3.山東泰豐控股集團有限公司王家寨煤礦，山東泰安271204）

由于參數(shù)模型具有很強的假設(shè)性，需要滿足一些假設(shè)條件，比如：觀測誤差中只有偶然誤差，平差模型線性化后的微小量可忽略不計等。所以，當(dāng)參數(shù)假定與實際偏差較大時，平差模型的參數(shù)估值精度降低。由于參數(shù)模型的局限性，相關(guān)學(xué)者引入非參數(shù)模型理論。非參數(shù)模型取消了變量之間在函數(shù)關(guān)系形式上的假設(shè)，其函數(shù)模型表示任意的函數(shù)關(guān)系，因此放寬了函數(shù)形式的限制［1-3］，但是非參數(shù)模型的模型解釋能力較低且觀測數(shù)據(jù)的維度受到限制。因此，針對參數(shù)模型與非參數(shù)模型各自的局限性，相關(guān)學(xué)者將二者結(jié)合到一起構(gòu)建了半?yún)?shù)模型。

在20世紀80年代，Engle等針對氣候條件改變對電力需求變化的影響提出了半?yún)?shù)模型［4-6］。半?yún)?shù)模型作為一種重要的統(tǒng)計模型，數(shù)學(xué)模型中包含參數(shù)部分與非參數(shù)部分，克服了單一的參數(shù)模型與非參數(shù)模型的不足［7］。當(dāng)觀測值中存在無法忽略的模型誤差時，半?yún)?shù)模型在數(shù)學(xué)模型方面相對于參數(shù)模型更符合客觀實際；在數(shù)值求解方面可以分別求解參數(shù)向量、非參數(shù)向量和偶然誤差，因此被廣泛應(yīng)用于測量數(shù)據(jù)處理中。

本研究采用礦區(qū)地表沉降實測數(shù)據(jù)建立半?yún)?shù)灰色Verhulst 模型，考慮到模型誤差問題，灰色Verhulst 模型的時間響應(yīng)序列的系數(shù)求取可以采用半?yún)?shù)模型，根據(jù)半?yún)?shù)模型估值準則設(shè)計不同計算方案構(gòu)建半?yún)?shù)灰色Verhulst 模型，求解時間響應(yīng)序列中的參數(shù)，并進行沉降預(yù)計結(jié)果對比和精度分析。

1 半?yún)?shù)灰色Verhulst模型

1.1 半?yún)?shù)模型

通常影響觀測值的因素有很多，然而經(jīng)典平差模型不考慮系統(tǒng)誤差與粗差的影響。當(dāng)模型誤差無法忽略且其性態(tài)十分復(fù)雜時，可以在函數(shù)模型中加入一個待定量，用來描述觀測量與影響因素之間函數(shù)關(guān)系不明確的部分，這種平差模型即為半?yún)?shù)模型［4，7］。

半?yún)?shù)模型的函數(shù)模型和隨機模型分別為

式中，L 為n 維觀測向量；X 為t 維參數(shù)向量；B 為列滿秩設(shè)計矩陣；S 為n 維非參數(shù)向量；Δ 為n 維誤差向量；為方差因子；P為Δ的權(quán)矩陣。

半?yún)?shù)模型在函數(shù)模型表達式中添加了非參數(shù)部分，相對于參數(shù)模型而言，在模型誤差無法忽略時，半?yún)?shù)模型具有更大的優(yōu)越性：在模型建立方面，數(shù)學(xué)模型與客觀實際更為接近；在數(shù)值解算方面，可求出模型誤差與偶然誤差。

由式（1）可知，半?yún)?shù)模型的誤差方程為

在最小二乘準則VTPV=min 下構(gòu)造條件極值函數(shù)得到的法方程為

由于待定參數(shù)個數(shù)多于誤差方程個數(shù)，因此基于最小二乘準則構(gòu)建的半?yún)?shù)模型的法方程唯一解不存在，需要在VTPV=min基礎(chǔ)上附加最優(yōu)化準則，構(gòu)成補償最小二乘法則［8］：

式中，R 為正則矩陣；α 為正則化參數(shù)（或稱為平滑因子），其主要作用是在極小化過程中，在二次型VTPV和起到平衡作用［9-10］。

根據(jù)式（2）和式（4）以及構(gòu)造條件極值的拉格朗日乘數(shù)法，構(gòu)造的目標(biāo)函數(shù)為

式中，λ1為拉格朗日常數(shù)。

由式（6）可得：

結(jié)合式（2）和式（7），可得到基于補償最小二乘估計準則的半?yún)?shù)模型誤差方程的法方程：

由于矩陣R為正定矩陣，因此基于補償最小二乘估計構(gòu)建的法方程是正定可逆矩陣，法方程有唯一解，對式（8）求解可得：

式中，M為光滑矩陣；H為帽子矩陣。

正則矩陣選取根據(jù)實際情況，構(gòu)造方法一般有時間序列法、距離法和樣條函數(shù)法等；正則矩陣選取一般有L曲線法、信噪比值法、效率法等［7，12］。

當(dāng)法方程系數(shù)矩陣接近奇異時，應(yīng)當(dāng)采用半?yún)?shù)嶺估計準則：

式中，k 是在X 與V 的正則化參數(shù)，稱為嶺參數(shù)，嶺參數(shù)的構(gòu)造方法一般有L曲線法［7，13］、U曲線法［14-15］等，參數(shù)求解方法與式（6）相同，在此不做贅述；Q一般取單位矩陣。

1.2 灰色Verhulst模型

設(shè)原始監(jiān)測數(shù)據(jù)序列為

通過原始數(shù)據(jù)累加弱化原始數(shù)據(jù)序列的波動性和隨機性，生成1-AGO序列

根據(jù)GM（1，1）模型的建模方法，得到灰色Verhulst模型的一階白化非線性微分模型為

式中，a,b為灰色Verhulst模型待計算參數(shù)。

若按最小二乘法，解得待求參數(shù)估值為

將參數(shù)估值a,b 代入微分方程求解得到灰色Verhulst模型的時間響應(yīng)序列［16］：

最后還原數(shù)列得到模型的預(yù)測值：

從測量平差的角度分析，式（14）求取待求參數(shù)與參數(shù)平差的法方程求解形式相同。在只考慮偶然誤差的情況下，可以采用最小二乘準則求解灰色Verhulst 模型一階白化非線性微分方程中的參數(shù)最優(yōu)估值。若模型誤差不可忽略，可考慮采用半?yún)?shù)模型求解一階白化非線性微分方程的參數(shù)估值，進而建立半?yún)?shù)灰色Verhulst模型進行預(yù)計。

1.3 半?yún)?shù)灰色Verhulst模型檢驗

建立半?yún)?shù)灰色Verhulst 模型時應(yīng)當(dāng)進行模型檢驗，驗證所建立的模型是否符合實際情況。半?yún)?shù)灰色Verhulst 模型檢驗主要包括采用半?yún)?shù)模型求解一階白化非線性微分方程中參數(shù)估值的假設(shè)檢驗以及灰色模型的精度檢驗。其中，半?yún)?shù)模型檢驗采用線性假設(shè)法［17］，灰色模型精度檢驗結(jié)果等級劃分如表1所示［18］。

?

2 算例分析

2.1 試驗數(shù)據(jù)

以某礦區(qū)地表沉降監(jiān)測為例，跟據(jù)相關(guān)設(shè)計，礦區(qū)九采區(qū)共有5 個回采工作面，其中，7265 工作面回采走向長度D1=909 m，傾向長度D2=186 m，煤層厚度m=4.15 m，煤層傾角θ=0～12°（平均7°），傾斜面積168 334 m2，工業(yè)儲量106 萬t，工作面標(biāo)高-678～-743 m，地面高程34.1～37.5 m，2017年2月中旬開采，2018 年1 月12 日回采結(jié)束，回采量88.2 萬t，回采率85%。

依據(jù)《測量規(guī)程》和“三下”開采規(guī)范，結(jié)合九采區(qū)的覆巖巖性，選取巖移角值參數(shù)布設(shè)觀測線。

7265 工作面觀測期間，觀測站97 號測點具有最大下沉速度Vmax，最大下沉值達到Wmax=1 241 mm。選取97 號沉降監(jiān)測點的14 期觀測數(shù)據(jù)得到的累計沉降值W（表2）建立預(yù)計模型。以原始數(shù)據(jù)為基礎(chǔ)建立預(yù)計模型，并以最后4期數(shù)據(jù)進行預(yù)計值與實際值對比，作為精度檢核依據(jù)。

?

2.2 計算結(jié)果對比及分析

若根據(jù)礦區(qū)沉降實測數(shù)據(jù)建立半?yún)?shù)灰色Verhulst 模型，首先要判定系數(shù)矩陣B 的復(fù)共線性［19-20］。經(jīng)計算，系數(shù)矩陣B 存在較強的復(fù)共線性，因此半?yún)?shù)模型中估值準則只需要考慮半?yún)?shù)嶺估計準則。

基于半?yún)?shù)灰色Verhulst 模型在無粗差數(shù)據(jù)時設(shè)計了以下3種計算方案：①灰色Verhulst模型；②半?yún)?shù)灰色Verhulst 模型，其中正則矩陣采用時間序列一階差分法；③半?yún)?shù)灰色Verhulst 模型，其中正則矩陣采用時間序列二階差分法。正則化參數(shù)按Xu函數(shù)法確定取值α=2.436，嶺參數(shù)按U 曲線法確定取值k=0.334。經(jīng)過不同方案計算，并以最后4 期的實測值作為預(yù)計效果的評價標(biāo)準，得到了不同平差模型的預(yù)計結(jié)果以及殘差值如表3所示。

?

由表3可知：通過采用不同半?yún)?shù)回歸模型對97號監(jiān)測點進行累計沉降預(yù)計時，方案1預(yù)計誤差絕對值的最大值為33 mm，方案2 預(yù)計誤差絕對值的最大值為15 mm，方案3 預(yù)計誤差絕對值的最大值為34 mm，3 種方案的均方誤差經(jīng)計算分別為37.68 mm、14.82 mm 和35.92 mm。因此方案2 建立的半?yún)?shù)灰色Verhulst 模型求取時間響應(yīng)式參數(shù)的精度最高，在3 種方案中的均方誤差最小，相應(yīng)的預(yù)計下沉曲線更加符合實測數(shù)據(jù)的變化趨勢。

由表3 中方案2 和方案3 對比可知：當(dāng)正則矩陣采用時間序列法確定時，方案2得到的預(yù)計結(jié)果優(yōu)于方案3。因此，本次實測數(shù)據(jù)中采用半?yún)?shù)模型的估值準則計算灰色Verhulst 模型一階白化非線性微分方程的待求參數(shù)時，應(yīng)選擇時間序列一階差分法確定的正則矩陣。

在上述分析的基礎(chǔ)上，本研究進一步對實測數(shù)據(jù)進行半?yún)?shù)模型檢驗與灰色模型精度檢驗。由線性假設(shè)法計算可知，統(tǒng)計量F=87.642。由F分布表可知：F ＞F0.05( 14, 13 )＞F0.05( 15, 13 )=2.53，因此非參數(shù)向量S 為非零向量。同時，經(jīng)計算相對誤差q=0.032，方差比C=0.394，小誤差概率P=0.615，對照灰色模型精度檢驗表（表1）可知：灰色模型精度滿足Ⅱ級（優(yōu)）標(biāo)準。因此，礦區(qū)地表沉降實測數(shù)據(jù)適用于構(gòu)建半?yún)?shù)灰色Verhulst模型。

3 結(jié) 語

采用礦區(qū)沉降實測數(shù)據(jù)進行了半?yún)?shù)灰色Verhulst 模型的算法驗證，基于半?yún)?shù)模型理論和礦區(qū)沉降實測數(shù)據(jù)建立了半?yún)?shù)灰色Verhulst 模型。針對灰色Verhulst 模型的時間響應(yīng)序列的系數(shù)求取問題，根據(jù)正則矩陣、正則化參數(shù)以及嶺參數(shù)的選取方法，設(shè)計了多種方案計算半?yún)?shù)灰色Verhulst 模型的參數(shù)估值，并進行了沉降預(yù)計結(jié)果對比和精度分析。在上述分析的基礎(chǔ)上，采用模擬數(shù)據(jù)與礦區(qū)沉降實測數(shù)據(jù)進行了模型檢驗。檢驗結(jié)果表明：礦區(qū)沉降實測數(shù)據(jù)適用于半?yún)?shù)灰色Verhulst 模型。實際應(yīng)用中，關(guān)于地表沉陷動態(tài)預(yù)計的時間函數(shù)還有Knothe時間函數(shù)、Logistic 時間函數(shù)等，對于不同地表沉陷動態(tài)預(yù)計的時間函數(shù)可進行進一步討論。