王曉輝，姜佃高

(太原理工大學測繪科學與技術系，山西太原030024)

一、引 言

開采沉陷預計可以定量地研究開采影響的地表移動變形在時間和空間上的分布規(guī)律，對“三下”開采也具有重要指導作用。而開采沉陷預計方法中應用最為廣泛的是概率積分法，其預計結果的精度主要取決于預計參數的可靠性［1］。因此，準確確定概率積分法參數對開采沉陷理論和生產實踐都具有重要意義。

概率積分法的預測參數一般根據地表移動觀測站實測資料通過最小二乘曲線擬合確定，最小二乘法有極好的配賦誤差的能力，但它易受粗差的影響，使參數估值失真［2-4］。當觀測值中不可避免地存在粗差時［5-6］，穩(wěn)健估計求參技術具有有效抵御粗差或異值對參數干擾的特性，能夠克服最小二乘法擬合求參時常出現的結果發(fā)散問題，保證求參結果的可靠性［7-8］。然而，不同穩(wěn)健估計方法的穩(wěn)健特性不同。Mitra等［9］研究表明，考慮到重尾噪聲分布，Huber法和Andrews法的穩(wěn)健性要優(yōu)于L1法。Li等［10］指出Fair法的穩(wěn)健性要優(yōu)于 L1法。Pennacchi［11］通過算例表明，在迭代 100次的情況下，Cauchy法比German-McClure法、Welsch法和Tukey法效果更好。方俊濤等［12］研究認為，Welsch法的穩(wěn)健性最好，其次是Tukey法，Huber法的穩(wěn)健性最差。李惠芬等［13］分析指出，Andrews法、Tukey法和IGG方案均優(yōu)于Fair法。

本文在嘉樂泉煤礦實測資料基礎上進行人工異值干擾求參試驗，以觀測值中含有不同粗差值為例，比較了13種常用穩(wěn)健估計方法求取概率積分法參數時的穩(wěn)健特性。

二、概率積分法參數的穩(wěn)健估計

1.概率積分法參數的穩(wěn)健估計模型

概率積分法參數的穩(wěn)健估計模型詳細見文獻［14］。

2.改進算法麥夸爾特法

用曲線擬合求取實測參數的基本算法是泰勒級數展開法，但其存在一個缺點，即當各參數初始值選的偏離其真值過大時，各次求得的各參數值可能不逐漸趨近于真值(稱為迭代收斂)，而是越來越偏離真值(稱為迭代發(fā)散)。為了克服這個缺點，放寬對初始值的限制，采用其改進算法，即麥夸爾特法。

首先選取參數 bi(i=1，2，…，m)的初始值可從經驗值中選取，也可根據本礦區(qū)已有的實參數據求出，則有(Δi為第i個參數的修正值)

將f(Xk;B)附近按泰勒級數展開，并假設與bi足夠接近，使Δi足夠小，并可在展開式中略去Δi的二次及二次以上的項，則有

式中

根據最小二乘原理VΤV=min，則有

即

寫出m個未知數Δi的線性方程組，并進行整理、移項，則上式可寫作

式中

上述中的d稱為阻尼因子，是一個≥0的可調整的常數，其作用是適當選擇d可使每次迭代后求得的Q'值均較迭代前的值小。d值的選擇原則是:在迭代收斂時，d選取較小的值，以減少迭代次數;在迭代后的Q'值大于迭代前的值時，將d值放大，重新進行迭代，直至迭代后的Q'值比迭代前的值小為止，這樣能保證經過一次迭代Q'就減少一次，bi的值就向其真值接近一步，從而大大放寬了對接近bi的要求，最終可求得各參數值bi。對應于概率積分求參下的參數有下沉系數q，主要影響角的正切 tanβ，主要影響傳播角 θ0和點偏移距 s1、s2、s3、s4等參數。

三、常用穩(wěn)健估計方法和穩(wěn)健性比較

1.常用穩(wěn)健估計方法

13種常用的穩(wěn)健估計方法及它們的調和系數見文獻［15］。

2.兩種參數估計方法的比較

為了比較任意兩種參數估計方法在人工異值干擾試驗下哪種方法更有效，給出參數的相對減益的概念，簡述如下:

定義:在概率積分求參中，當觀測值中不存在粗差時，通過最小二乘曲線擬合法求得的參數理論上是最優(yōu)值，而當觀測值中存在粗差時通過各種穩(wěn)健估計方法求得同一參數的值與無粗差的最小二乘曲線擬合法求得的參數值會發(fā)生變化，二者之間的差值與無粗差的最小二乘曲線擬合法求得的參數比值的絕對值的百分數就稱作參數的相對減益。當用穩(wěn)健估計方法求得參數的相對減益越大，說明這種方法越不穩(wěn)健;反之，其值越小則越穩(wěn)健。

在概率積分求參中選取3個主要參數:下沉系數q、主要影響角的正切tanβ和拐點偏移距D=(H0為平均采深)。令下沉系數q的相對減益為

主要影響角的正切tanβ的相對減益為

拐點偏移距D的相對減益為

上式中，當觀測值中不存在粗差時，通過最小二乘曲線擬合法求得的3個主要參數值為下沉系數q0、主要影響角的正切tanβ0和拐點偏移距D0;當觀測值中存在粗差時，通過各種穩(wěn)健估計方法求得的3個主要參數值為下沉系數qi、主要影響角的正切tan βi和拐點偏移距 Di(i=1，2，3，…，13)，i為所選的第i種穩(wěn)健估計方法。

四、結果與討論

基于嘉樂泉煤礦某工作面地表移動走向觀測線實測資料，通過人為地增加一些異值點的辦法，進行人工干擾穩(wěn)健求參試驗。該工作面采用走向長壁采煤法采煤，全部陷落法管理頂板，上覆巖層巖性綜合評價為中硬，地表沉陷規(guī)律基本符合概率積分模型，地表線A是半無限開采，其實測下沉值見表1。本文采用將走向觀測線的拐點和最大下沉點處第一次加上100 mm粗差，第二次加上200 mm的粗差，經計算得出了13種常用穩(wěn)健估計方法與無粗差下最小二乘法3個主要參數的相對減益，其結果如圖1—圖3所示。

表1 走向觀測線地表下沉值

圖1 不同穩(wěn)健估計方法相對于無粗差的最小二乘法下沉系數的相對減益

圖2 不同穩(wěn)健估計方法相對于無粗差的最小二乘法拐點偏移距的相對減益

圖3 不同穩(wěn)健估計方法相對于無粗差的最小二乘法的主要影響角正切的相對減益

由圖1可看出，當粗差取100 mm時，L1法和German-McClure法下沉系數的相對減益均為0%，其次是IGGIII方案，其相對減益為1%，其余方法的結果都≥2%;當粗差取200 mm時，L1法和IGGIII方案相對減益均為6%，而German-McClure法下沉系數的相對減益為7%，其余方法的結果都≥8%。

由圖2可看出，當粗差取100 mm時，L1法和German-McClure法拐點偏移距的相對減益均為0%，而IGGIII方案的相對減益為1%，其余方法的結果都≥2%;當粗差取200 mm時，L1法和IGGIII方案相對減益均為4%，而German-McClure法拐點偏移距的相對減益為5%，其余方法的結果都≥6%?？傮w上看各個穩(wěn)健估計方法的拐點偏距的相對減益相差不是很大，都在9%以內。

由圖3可看出，當粗差取100 mm時，L1法、L1-L2法和German-McClure法主要影響角正切的相對減益均為 0%，而是 Huber法、Andrews法、Welsch法、Tukey法、Danish法、IGG 方案、IGGIII方案和Cauchy法主要影響正切的相對減益為1%，各方法的結果都特別接近，都不大于3%;當粗差取200 mm時，L1法、German-McClure法和IGGIII方案相對減益均為4%，其余方法的結果都≥6%。

五、結束語

綜合圖1、圖2和圖3的試驗結果，可知當給走向觀測線的拐點和最大下沉點處同時加上100 mm的粗差時，對于各穩(wěn)健估計方法3個主要參數的相對減益相差都不是很大，其中 L1法、German-McClure法和IGGIII方案相對更穩(wěn)健;當給走向觀測線的拐點和最大下沉點處同時加上200 mm的粗差時，各穩(wěn)健估計主要參數相對減益值差異逐漸增大，而 L1法、German-McClure法和IGGIII方案的3個主要參數相對減益值比其他穩(wěn)健估計方法的都小。對于概率積分法參數的穩(wěn)健估計方法，L1、German-McClure方法和IGGIII方案較其他穩(wěn)健估計方法相對更穩(wěn)健。

［1］何國清，楊倫，凌賡娣，等.礦山開采沉陷學［M］.徐州:中國礦業(yè)大學出版社，1994.

［2］王素立，劉永.基于波動相關性及主分量變換的多元線性回歸模型研究［J］.統計與決策，2012(22):18-21.

［3］高山，李孝軍.關于最小二乘法的穩(wěn)健性分析［J］.統計與決策，2006(8):125-126.

［4］何霞，劉衛(wèi)鋒.基于全最小一乘準則的灰色GM(1，1)模型參數估計［J］.統計與決策，2006(8):30-33.

［5］GRIEP M I，WAKELING I N，VANKEERBERGHEN P，et al.Comparison of Semi-robust and Robust Partial Least Squares Procedures［J］.Chemometrics and Intelligent Laboratory Systems，1995，29(1):37-50.

［6］GUO JF，OU J K，WANG H T.Robust Estimation for Correlated Observations:Two Local Sensitivity-based Down Weighting Strategies ［J］.Journal of Geodesy，2010，84(4):243-250.

［7］吳健平，張立.地理數據線性回歸中的穩(wěn)健估計方法［J］.干旱區(qū)地理，1994，17(1):83-88.

［8］馮守平，石澤，鄒瑾.一元線性回歸模型中參數估計的幾種方法比較［J］.統計與決策，2008(24):152-153.

［9］MITRA S，MITRA A，KUNDU D.Genetic Algorithm and M-Estimator Based Robust Sequential Estimation of Parameters of Nonlinear Sinusoidal Signals［J］.Communications in Nonlinear Science and Numerical Simulation，2010，16(7):2796-2809.

［10］LI H F，JIANG X Q，LI Z.Robust Estimation in Gaussian Filtering for Engineering Surface Characterization［J］.Precision Engineering，2004，28(2):186-193.

［11］PENNACCHI P.Robust Estimate of Excitations in Mechanical Systems Using M-Estimators-Theoretical Background and Numerical Applications［J］.Journal of Sound and Vibration，2008，310(4-5):923-946.

［12］方俊濤，何楨，宋琳曦.響應曲面建模的穩(wěn)健M-回歸方法［J］.工業(yè)工程，2012，15(3):98-103.

［13］李惠芬，蔣向前，李柱.高斯濾波穩(wěn)健性能的研究與改進［J］.儀器儀表學報，2004，25(5):633-637.

［14］郭廣云，汪云甲.概率積分法參數的穩(wěn)健估計模型及其應用研究［J］.測繪學報，2000，29(2):162-165.

［15］GE Y H，YUAN Y，JIA N N.More Efficient Methods among Commonly Used Robust Estimation Methods for GPS Coordinate Transformation ［J］.Survey Review，2013，45(330):229-234.