抗差岭估计在概率积分法预计参数求取中的应用研究

王 友，王双亭

(河南理工大学测绘与国土信息工程学院，河南焦作454003)

基于随机介质理论的概率积分法因其所用的移动和变形预计公式中含有概率积分而得名，是目前我国用于矿山开采沉陷预计地表移动与变形的重要方法和文献［1］规定的开采沉陷预计方法之一。

在求取预计参数时，通常采用计算简便的最小二乘拟合算法，其求取的参数精度基本能满足工程需求。最小二乘估计在参数估计时具有良好的性质，当误差服从正态分布时，最小二乘估计在所有无偏估计类中具有无偏性、一致性和有效性。但是其也存在两方面的问题:一是当自变量较多，其中存在近似线性相关变量时，其参数估值与真值相差很大;二是当观测值有悖于正态分布假设，数据遭到异常污染时，最小二乘估计不具有抗干扰性，单个观测值的偏差就可能造成参数解面目全非［2－4］。

针对第1个问题，文献 ［2］提出了岭估计，岭估计是从减少均方误差的角度出发而提出的一种压缩性有偏估计，其可以改善法矩阵的病态性，稳定参数解;针对第2个问题，文献 ［4］和 ［5］提出了抗差估计，其通过选择适当的等价权，使其可以较好地克服模型偏差和异值点存在造成的求参困难。但二者都是只能改善某一个问题，不能同时对两个问题都有所改善。本文提出的抗差岭估计结合了岭估计和抗差估计的优点，可以同时起到抵抗病态法矩阵和抵抗异值或粗差对求参结果的影响，保证了求参结果的有效性和可靠性。

1 抗差岭估计

1.1 抗差岭估计原理［3，6］

设观测方程为

式中，L为n维观测向量;A为n×t阶系数矩阵; X为t维参数向量;Δ为n维误差向量，Δ:N(0，I)。

相应的误差方程为

由抗差岭估计原理［6－7］求得参数的抗差岭估计解为

由式 (3)可知，确定岭参数和等价权阵是求得抗差岭估计解的必要条件。所以抗差岭估计的关键是选择合适的岭参数和等价权阵。

1.2 岭参数的确定［2，7－10］

对岭参数的确定方法进行对比研究是当前的热点问题之一。关于岭参数的确定，目前应用较多的有岭迹法、GCV解法、双h公式法和L曲线法等等。其中岭迹法虽然应用方便，但具有太大的随意性，且没有严格的理论基础;GCV解法虽然在理论上能够选择最优的岭参数，但有时GCV函数的变化过于平缓甚至发散，此时定位它的最小值会很困难;双h公式法虽然在确定岭参数时有计算简捷、应用灵活等优点，但当法方程系数矩阵病态性较严重时，其效果并不明显;L曲线法是一种理论上非常严密的研究岭参数的方法，其具有定位准确、适用性较好等优点。因此，本文在确定岭参数时选用L曲线法。

按正则化原理，式 (1)的抗差岭估计准则为

对上式求其最大值，就可以得到最大曲率，其所对应的最大点即为所求点，然后求出该点所对应的抗差岭参数k。

1.3 等价权函数的确定［4－5，11－12］

抗差岭估计的抗差效果主要取决于等价权函数。采用的等价权函数不同，则相应的抗差估计模型就不同，其抗差效果就不同。

周江文教授提出的IGGN方案是对最小二乘参数估计解式的权阵进行改化，即以等价权代替先验权，从而可以用等价权重构抗差解式和单位权中误差。本文中采用周江文教授提出的IGGN抗差方案。等价权的形式为

式中，k0为分位点，一般取1.5;k1为淘汰点，一般取2.5。

但是由于误差绝对值大于中误差、二倍中误差、三倍中误差出现的概率分别为 31.7%，4.5%，0.3%，加之我国大多测量规范都规定以二倍中误差为极限误差，因此，我们取k0=1.0，k1= 2.0。

2 概率积分法预计参数的抗差岭估计模型

概率积分法预计参数可由多个地表任意点沿任意方向的移动变形实测资料求出。当布置的测点为常规观测站时，采用最小二乘曲线拟合法;当测点为非常规观测站即一系列散点时，采用最小二乘曲面拟合法。但是当最小二乘的法矩阵病态或观测数据中有异常值时，则易造成求参结果的失真。基于抗差岭估计的最小二乘曲面拟合法在求参时，不仅可以有效地抵御异常值和病态法矩阵的干扰，还可以适用于任意形状和残缺观测站资料的求参运算［2］。

概率积分法中所涉及的主要预计参数有5个:下沉系数q、水平移动系数b、主要影响角正切tanβ、主要影响传播角θ、拐点偏移距S1，S2，S3，S4。在进行抗差岭估计求参时，应首先利用实测的地表点下沉观测值求取q，tanβ，θ，S1，S2，S3，S4。然后，利用这7个已求得的参数和地表点的水平移动观测值迭代拟合求取b。

按概率积分法的任意点移动变形模型，任一点下沉W都可以表示成测点坐标和预计参数的函数，即

选取参数初值 q0，tanβ0，θ0，S10，S20，S30，S40，将其线性化，即

则得误差方程的一般形式为:

式中，l0i=Wi测－W0i。

式(9)的矩阵形式为

由式(5)确定岭参数k，由式(6)确定等价权，则可求出7个参数的改正数的抗差岭估计解式如式 (3)所示。求得结果后，以上次计算的结果做为下次的初值，进行迭代计算，直到最后2次迭代计算的参数估值之差满足迭代收敛精度为止。此时即可获得参数估值的抗差岭估计解。而水平移动系数b，则可以根据地表任意点水平移动预测模型求出。在此不再赘述。

3 工程应用及结果分析

工程实例采用文献［13］－［14］中开滦集团钱家营煤矿辅271地表移动观测站资料。该观测站位于辅271工作面上方。由于辅271观测站的实测资料分单一煤层和多煤层开采2种情况，而在开采单一煤层时走向和倾向观测站实测资料比较翔实，因此采用辅271观测站单一煤层开采时的资料。该观测站所对应的辅271工作面表土层厚220m，采厚3.1m，平均采深304m，煤层倾角为8°，工作面走向长830m，倾向长150m。

为了检验抗差岭估计模型的有效性，首先分别运用最小二乘法和抗差岭估计算法拟合求参，然后利用 matlab软件中的 randn命令生成标准差为10mm的随机误差加入到观测值中，再分别用最小二乘法和抗差岭估计算法拟合求参，最后对2次求参结果进行纵横向对比研究，求参结果如表1、表2所示。

表1 人为干预前数据拟合求参结果

表2 人为干预后数据拟合求参结果

对表1和表2对比分析可知:

(1)当观测值中没有人工干预时，抗差岭估计和传统最小二乘估计理论的求参结果基本一致。

(2)当人为地在观测值中加入一些异值点和随机误差时，采用传统最小二乘估计理论求参获得的参数有很大的偏差，而采用抗差岭估计理论获得的参数与没有人为干预时获得的参数相差不大。

4 结论

(1)将抗差岭估计理论应用于概率积分法预计参数求取过程中，可以自动排除病态法矩阵和异值点的干扰。

(2)通过人工干预求参试验表明，采用抗差岭估计求参理论，求得的概率积分法预计参数具有更好的有效性和可靠性。

［1］中华人民共和国煤炭工业部.建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程［M］.北京:煤炭工业出版社，2000.

［2］王明柱，郭广礼，王 磊，等.基于岭估计的概率积分法预计参数的求取［J］.煤矿开采，2012，17(2):17－19，85.

［3］隋立芬.抗差岭估计原理及其应用 ［J］.测绘通报，1994 (1):9－12.

［4］吴若飞.基于抗差估计的概率积分法的预计参数模型研究［J］.煤炭技术，2009，28(9):167－168.

［5］顾 叶，宋振柏，张胜伟.基于概率积分法的开采沉陷预计研究［J］.山东理工大学学报 (自然科学版)，2011，25 (1):33－36.

［6］归庆明，黄维彬，张建军.抗差泛岭估计［J］.测绘学报，1998(3):211－216.

［7］王 彬，高井祥，刘 超，等.L曲线法在等价权抗差岭估计模型中的应用 ［J］.大地测量与地球动力学，2012，32 (3):97－101.

［8］王振杰.用L－曲线法确定岭估计中的岭参数［J］.武汉大学学报 (信息科学版)，2004，29(3):235－238.

［9］黄海兰，牛 犇.岭参数确定的研究［J］.测绘科学，2011，36(4):31－32.

［10］田玉淼，朱建军，陶肖静.修正岭估计方法在测量数据处理中的应用研究［J］.测绘工程，2012，21(1):7－10.

［11］杨 玲，沈云中，楼立志.基于中位参数初值的等价权抗差估计方法［J］.测绘学报，2011，40(1):28－32.

［12］陈西强，黄张裕.抗差估计的选权迭代法分析与比较［J］.测绘工程，2010，19(4):8－11，15.

［13］殷作如，邹友峰，邓智毅，等.开滦矿区岩层与地表移动规律及参数［M］.北京:科学出版社，2010.

［14］陈 勇.开滦矿区深部开采地表移动规律的研究［D］.河南:河南理工大学，2010.