王曉輝,姜佃高
(太原理工大學(xué)測繪科學(xué)與技術(shù)系,山西太原030024)
開采沉陷預(yù)計(jì)可以定量地研究開采影響的地表移動(dòng)變形在時(shí)間和空間上的分布規(guī)律,對“三下”開采也具有重要指導(dǎo)作用。而開采沉陷預(yù)計(jì)方法中應(yīng)用最為廣泛的是概率積分法,其預(yù)計(jì)結(jié)果的精度主要取決于預(yù)計(jì)參數(shù)的可靠性[1]。因此,準(zhǔn)確確定概率積分法參數(shù)對開采沉陷理論和生產(chǎn)實(shí)踐都具有重要意義。
概率積分法的預(yù)測參數(shù)一般根據(jù)地表移動(dòng)觀測站實(shí)測資料通過最小二乘曲線擬合確定,最小二乘法有極好的配賦誤差的能力,但它易受粗差的影響,使參數(shù)估值失真[2-4]。當(dāng)觀測值中不可避免地存在粗差時(shí)[5-6],穩(wěn)健估計(jì)求參技術(shù)具有有效抵御粗差或異值對參數(shù)干擾的特性,能夠克服最小二乘法擬合求參時(shí)常出現(xiàn)的結(jié)果發(fā)散問題,保證求參結(jié)果的可靠性[7-8]。然而,不同穩(wěn)健估計(jì)方法的穩(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]研究認(rèn)為,Welsch法的穩(wěn)健性最好,其次是Tukey法,Huber法的穩(wěn)健性最差。李惠芬等[13]分析指出,Andrews法、Tukey法和IGG方案均優(yōu)于Fair法。
本文在嘉樂泉煤礦實(shí)測資料基礎(chǔ)上進(jìn)行人工異值干擾求參試驗(yàn),以觀測值中含有不同粗差值為例,比較了13種常用穩(wěn)健估計(jì)方法求取概率積分法參數(shù)時(shí)的穩(wěn)健特性。
概率積分法參數(shù)的穩(wěn)健估計(jì)模型詳細(xì)見文獻(xiàn)[14]。
用曲線擬合求取實(shí)測參數(shù)的基本算法是泰勒級數(shù)展開法,但其存在一個(gè)缺點(diǎn),即當(dāng)各參數(shù)初始值選的偏離其真值過大時(shí),各次求得的各參數(shù)值可能不逐漸趨近于真值(稱為迭代收斂),而是越來越偏離真值(稱為迭代發(fā)散)。為了克服這個(gè)缺點(diǎn),放寬對初始值的限制,采用其改進(jìn)算法,即麥夸爾特法。
首先選取參數(shù) bi(i=1,2,…,m)的初始值可從經(jīng)驗(yàn)值中選取,也可根據(jù)本礦區(qū)已有的實(shí)參數(shù)據(jù)求出,則有(Δi為第i個(gè)參數(shù)的修正值)
將f(Xk;B)附近按泰勒級數(shù)展開,并假設(shè)與bi足夠接近,使Δi足夠小,并可在展開式中略去Δi的二次及二次以上的項(xiàng),則有
式中
根據(jù)最小二乘原理VΤV=min,則有
即
寫出m個(gè)未知數(shù)Δi的線性方程組,并進(jìn)行整理、移項(xiàng),則上式可寫作
式中
上述中的d稱為阻尼因子,是一個(gè)≥0的可調(diào)整的常數(shù),其作用是適當(dāng)選擇d可使每次迭代后求得的Q'值均較迭代前的值小。d值的選擇原則是:在迭代收斂時(shí),d選取較小的值,以減少迭代次數(shù);在迭代后的Q'值大于迭代前的值時(shí),將d值放大,重新進(jìn)行迭代,直至迭代后的Q'值比迭代前的值小為止,這樣能保證經(jīng)過一次迭代Q'就減少一次,bi的值就向其真值接近一步,從而大大放寬了對接近bi的要求,最終可求得各參數(shù)值bi。對應(yīng)于概率積分求參下的參數(shù)有下沉系數(shù)q,主要影響角的正切 tanβ,主要影響傳播角 θ0和點(diǎn)偏移距 s1、s2、s3、s4等參數(shù)。
13種常用的穩(wěn)健估計(jì)方法及它們的調(diào)和系數(shù)見文獻(xiàn)[15]。
為了比較任意兩種參數(shù)估計(jì)方法在人工異值干擾試驗(yàn)下哪種方法更有效,給出參數(shù)的相對減益的概念,簡述如下:
定義:在概率積分求參中,當(dāng)觀測值中不存在粗差時(shí),通過最小二乘曲線擬合法求得的參數(shù)理論上是最優(yōu)值,而當(dāng)觀測值中存在粗差時(shí)通過各種穩(wěn)健估計(jì)方法求得同一參數(shù)的值與無粗差的最小二乘曲線擬合法求得的參數(shù)值會發(fā)生變化,二者之間的差值與無粗差的最小二乘曲線擬合法求得的參數(shù)比值的絕對值的百分?jǐn)?shù)就稱作參數(shù)的相對減益。當(dāng)用穩(wěn)健估計(jì)方法求得參數(shù)的相對減益越大,說明這種方法越不穩(wěn)健;反之,其值越小則越穩(wěn)健。
在概率積分求參中選取3個(gè)主要參數(shù):下沉系數(shù)q、主要影響角的正切tanβ和拐點(diǎn)偏移距D=(H0為平均采深)。令下沉系數(shù)q的相對減益為
主要影響角的正切tanβ的相對減益為
拐點(diǎn)偏移距D的相對減益為
上式中,當(dāng)觀測值中不存在粗差時(shí),通過最小二乘曲線擬合法求得的3個(gè)主要參數(shù)值為下沉系數(shù)q0、主要影響角的正切tanβ0和拐點(diǎn)偏移距D0;當(dāng)觀測值中存在粗差時(shí),通過各種穩(wěn)健估計(jì)方法求得的3個(gè)主要參數(shù)值為下沉系數(shù)qi、主要影響角的正切tan βi和拐點(diǎn)偏移距 Di(i=1,2,3,…,13),i為所選的第i種穩(wěn)健估計(jì)方法。
基于嘉樂泉煤礦某工作面地表移動(dòng)走向觀測線實(shí)測資料,通過人為地增加一些異值點(diǎn)的辦法,進(jìn)行人工干擾穩(wěn)健求參試驗(yàn)。該工作面采用走向長壁采煤法采煤,全部陷落法管理頂板,上覆巖層巖性綜合評價(jià)為中硬,地表沉陷規(guī)律基本符合概率積分模型,地表線A是半無限開采,其實(shí)測下沉值見表1。本文采用將走向觀測線的拐點(diǎn)和最大下沉點(diǎn)處第一次加上100 mm粗差,第二次加上200 mm的粗差,經(jīng)計(jì)算得出了13種常用穩(wěn)健估計(jì)方法與無粗差下最小二乘法3個(gè)主要參數(shù)的相對減益,其結(jié)果如圖1—圖3所示。
表1 走向觀測線地表下沉值
圖1 不同穩(wěn)健估計(jì)方法相對于無粗差的最小二乘法下沉系數(shù)的相對減益
圖2 不同穩(wěn)健估計(jì)方法相對于無粗差的最小二乘法拐點(diǎn)偏移距的相對減益
圖3 不同穩(wěn)健估計(jì)方法相對于無粗差的最小二乘法的主要影響角正切的相對減益
由圖1可看出,當(dāng)粗差取100 mm時(shí),L1法和German-McClure法下沉系數(shù)的相對減益均為0%,其次是IGGIII方案,其相對減益為1%,其余方法的結(jié)果都≥2%;當(dāng)粗差取200 mm時(shí),L1法和IGGIII方案相對減益均為6%,而German-McClure法下沉系數(shù)的相對減益為7%,其余方法的結(jié)果都≥8%。
由圖2可看出,當(dāng)粗差取100 mm時(shí),L1法和German-McClure法拐點(diǎn)偏移距的相對減益均為0%,而IGGIII方案的相對減益為1%,其余方法的結(jié)果都≥2%;當(dāng)粗差取200 mm時(shí),L1法和IGGIII方案相對減益均為4%,而German-McClure法拐點(diǎn)偏移距的相對減益為5%,其余方法的結(jié)果都≥6%。總體上看各個(gè)穩(wěn)健估計(jì)方法的拐點(diǎn)偏距的相對減益相差不是很大,都在9%以內(nèi)。
由圖3可看出,當(dāng)粗差取100 mm時(shí),L1法、L1-L2法和German-McClure法主要影響角正切的相對減益均為 0%,而是 Huber法、Andrews法、Welsch法、Tukey法、Danish法、IGG 方案、IGGIII方案和Cauchy法主要影響正切的相對減益為1%,各方法的結(jié)果都特別接近,都不大于3%;當(dāng)粗差取200 mm時(shí),L1法、German-McClure法和IGGIII方案相對減益均為4%,其余方法的結(jié)果都≥6%。
綜合圖1、圖2和圖3的試驗(yàn)結(jié)果,可知當(dāng)給走向觀測線的拐點(diǎn)和最大下沉點(diǎn)處同時(shí)加上100 mm的粗差時(shí),對于各穩(wěn)健估計(jì)方法3個(gè)主要參數(shù)的相對減益相差都不是很大,其中 L1法、German-McClure法和IGGIII方案相對更穩(wěn)健;當(dāng)給走向觀測線的拐點(diǎn)和最大下沉點(diǎn)處同時(shí)加上200 mm的粗差時(shí),各穩(wěn)健估計(jì)主要參數(shù)相對減益值差異逐漸增大,而 L1法、German-McClure法和IGGIII方案的3個(gè)主要參數(shù)相對減益值比其他穩(wěn)健估計(jì)方法的都小。對于概率積分法參數(shù)的穩(wěn)健估計(jì)方法,L1、German-McClure方法和IGGIII方案較其他穩(wěn)健估計(jì)方法相對更穩(wěn)健。
[1]何國清,楊倫,凌賡娣,等.礦山開采沉陷學(xué)[M].徐州:中國礦業(yè)大學(xué)出版社,1994.
[2]王素立,劉永.基于波動(dòng)相關(guān)性及主分量變換的多元線性回歸模型研究[J].統(tǒng)計(jì)與決策,2012(22):18-21.
[3]高山,李孝軍.關(guān)于最小二乘法的穩(wěn)健性分析[J].統(tǒng)計(jì)與決策,2006(8):125-126.
[4]何霞,劉衛(wèi)鋒.基于全最小一乘準(zhǔn)則的灰色GM(1,1)模型參數(shù)估計(jì)[J].統(tǒng)計(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]吳健平,張立.地理數(shù)據(jù)線性回歸中的穩(wěn)健估計(jì)方法[J].干旱區(qū)地理,1994,17(1):83-88.
[8]馮守平,石澤,鄒瑾.一元線性回歸模型中參數(shù)估計(jì)的幾種方法比較[J].統(tǒng)計(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]方俊濤,何楨,宋琳曦.響應(yīng)曲面建模的穩(wěn)健M-回歸方法[J].工業(yè)工程,2012,15(3):98-103.
[13]李惠芬,蔣向前,李柱.高斯濾波穩(wěn)健性能的研究與改進(jìn)[J].儀器儀表學(xué)報(bào),2004,25(5):633-637.
[14]郭廣云,汪云甲.概率積分法參數(shù)的穩(wěn)健估計(jì)模型及其應(yīng)用研究[J].測繪學(xué)報(bào),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.