088-基于岭估计的概率积分法预计参数的求取

U（ x，y，φ） = U（ x） Cycosφ + U（ y） Cxsinφ （ 13）
［关键词］ 岭估计； 概率积分法； 求参； 开采沉陷； L 曲线法 ［中图分类号］ TD325. 2 ［文献标识码］ A ［文章编号］ 1006-6225 （ 2012） 02-0017-03
Prediction Parameters Calculation of Probability Integration Method Based on Ridge Estimate
王明柱1，2 ，郭广礼1，2 ，王 磊1，2 ，王 彬1，2
（ 1. 中国矿业大学 江苏省资源环境信息工程重点实验室，江苏 徐州 221116； 2. 中国矿业大学 国土环境与灾害监测国家测绘局重点实验室，江苏 徐州 221116）
［摘 要］ 采用最小二乘拟合求取概率积分法预计参数时，经常遇到发散问题，致使求参结果 失真。而求参过程中法矩阵的病态是造成结果发散的主要原因之一。将岭估计应用于预计参数的求取 中，能够有效地克服传统最小二乘拟合求参时法矩阵的病态。在给出基于岭估计求参模型的同时，使 用 Matlab 实现了参数的自动求取。最后通过模拟算例验证了该求参模型克服法矩阵病态的有效性及求 参结果的可靠性。
概率积分法是进行开采沉陷预计的主要方法， 预计过程中预计参数的准确性决定了地表移动变形 预计的精度。参数的影响因素主要有： 覆岩的岩性 和地层结构、采深、采厚等［1］。我国 大部分矿区 都根据所设的地表移动观测站反演出了适合自己矿 区的概率积分预计参数。
预计参数的求取通常采用最小二乘拟合的方 法，该方法计算简单，且求参精度能满足工程需 求。而最小二乘拟合求参过程中法矩阵一旦病态， 就致使迭代结果发散、失真，造成求参困难。为解 决由于法方程系数阵病态而导致最小二乘估计不稳 定问题，统计学家们提出了岭估计、广义岭估计、 主成分估计等有偏估计来改善这种情况下的最小二 乘估计。岭估计法从减小均方误差的角度出发，最 初由 A. E. Hoerl 于 1962 年提出，并于 1970 年由他 和 R. W. Kennard 对该方法做了系统的发展。该方 法成为目前使用最多的一种有偏估计［2］。
数经第 1 次改正后变为 Z0 +Δ Z^ 1 ，记为Z^1 。再将Z^1
作为第 2 次计算的初值，求得经第 2 次改正后的参
数Z^2 。如此 迭 代，直 至 max | ^ ^ ΔZi+1 / Zi | ＜ Eps 停
止。Eps 为限制误差，一般可取 0. 0001 ～ 0. 001。
求出 q，tanβ0 ，θ，S1 ，S2 ，S3 ， S4 后， 认 为 这些参数已知。则任一点水平移动预计公式中只剩 下水平移动系数 b 未知。其预计公式为［1］
ρ，选择不同的 k 值将得到许多 （ ρ，η） 点。将这 些点经过曲线拟合得到一条曲线，选择该曲线曲率
最大的点所对应的岭参数即为所求。利用这条曲线
来选择岭参数的方法称为 L 曲线法。
令 η^ = log‖X^ ‖2 ，ρ^ = log‖AX^ －L‖2
则 η^ = 2log‖X^ ‖ = 2η，ρ^ = 2log‖AX^ －L‖ = 2ρ
靠。基于岭估计的曲面拟合法不仅可以有效地抵御
法矩阵病态的干扰，还可适用于非正规观测站和残
缺观测站资料的求参。
概率积分法的预计参数有下沉系数 q、主要影
响角正 切 tanβ、主 要 影 响 传 播 角 θ、拐 点 偏 移 距
S1，S2，S3，S4 及水平移动系数 b。其任一点移动 变形预计具体计算公式可参见文献 ［1］。进行基
采用岭估计可有效克服病态矩阵干扰，使求参
结果更可靠，确保概率积分法开采沉陷预计精度。
1 岭估计
1. 1 岭估计基本原理 设有观测方程：
AX = L+ Δ
（ 1）
式中，L 是观测值； A 是系数阵； X 是待估计参数 的真值； Δ 是噪声，Δ ～ N （ 0，δ20 I） 。则上式常规 的最小二乘解为
Wi = f（ xi ，yi ，q，tanβ，S1 ，S2 ，S3 ，S4 ，cotθ） （ 6） 首先选 取 参 数 初 值 q0 ，tanβ0 ，S10 ，S20 ，S30 ， S40 和 cotθ0 ，初值的向量形式记为 Z0 。其值可根据 文献 ［9］ 上的经验值进行选取，也可以参照临近
矿区已有的实测参数。将 （ 6） 式在 （ xi，yi，q0， tanβ0 ，S10 ，S20 ，S30 ，S40 ，cotθ0 ） 处按泰勒级数展 开，展开式中略去二次及以上项，则有
WANG Ming-zhu1，2 ，GUO Guang-li1，2 ，WANG Lei1，2 ，WANG Bin1，2
（ 1. Key Laboratory of Resources Environment Information Engineering of Jiangsu Province，China University of Mining ＆ Technology，Xuzhou 221116， China； 2. Key Laboratory of Land Environment and Disaster Monitoring of State Bureau of Surveying and Mapping，China University of Mining ＆ Technology，Xuzhou 221116，China）
^X = N－1 ATPL
（ 2）
若法矩阵 N = ATPA 条件数很大，即为病态矩
阵时，求逆就会出现不稳定，致使最小二乘估计的 方法不可 能 得 到 未 知 参 数 的 精 确 估 值［3］。 岭 估 计 法可解决 这 类 问 题［4］。 其 基 本 思 想 是 利 用 原 最 小
二乘估计的数学模型，在其法方程系数矩阵 N 的
‖AX^ －L‖2 + kΩ （ X^ ） = ‖AX^ －L‖2 + k X^TX^
= ‖AX^ －L‖2 +k‖X^ ‖2 = min
（ 4）
式中，‖·‖为欧氏 2－范数； Ω （ X^ ） 为稳定泛函。
‖AX^ －L‖和‖X^ ‖都是岭参数 k 的函数。以
log‖X^ ‖作为纵坐标 η，log‖AX^ －L‖作为横坐标
Wi = f（ xi ，yi ，q0 ，tanβ0 ，S10 ，S20 ，S30 ，S40 ，cotθ0 ） +
qf Δq
+
f tanβ
Δtanβ
+
f S1
ΔS1
+
f S2
ΔS2
+
f S3
ΔS3
+
f S4
ΔS4
+
f cotθ
Δcotθ
（ 7）
考虑模型误差和测量误差，式 （ 7） 可写为
Wi测
+ Vi
= f Δq q
+ f Δtanβ tanβ
+
f S1
ΔS1
+
f S2
ΔS2
+
f S3
ΔS3
+
f S4
ΔS4
+
Байду номын сангаас
f Δcotθ cotθ
+
W0i
（ 8）
式中，
W0i = f（ xi ，yi ，q0 ，tanβ0 ，S10 ，S20 ，S30 ，S40 ，cotθ0 ） （ 9） 则误差方程的一般形式为
关于岭参数的确定，目前已经有不少的方法，
如岭迹法、L 曲线法、两步解法及 GCV 法等，但
是到目前为止还没有一种公认的好方法。相比较而
言，L 曲线法易于确定岭参数，是一种确定岭参数 的良好方法［5］。本文使用 L 曲线法来确定岭参数， 在此进行简述，其他方法可参见文献 ［5－7］。
根据正则化原理，（ 1） 式的岭估计准则为
Vi
=
f Δq q
+
f Δtanβ tanβ
+
f S1
ΔS1
+
f S2
ΔS2
+
f S3
ΔS3
+
f S4
ΔS4
+
f cotθ
Δcotθ
－
l0i
（ 10）
式中，lwi = Wi测 － W0i ，7 个一阶偏导式详见文献
［10］。
上式的矩阵形式为
^ Vn×1 = An×7 Δ Z7×1 －ln×1
文献 ［8］。
2 概率积分法预计参数的岭估计模型
2. 1 求参的岭估计模型 概率积分法预计参数可由地表点实测移动变形
资料求出。当布置的测点为正规观测站时，可采用 最小二乘曲线拟合法； 当测点为非正规观测站即一 系列散点时，可采用最小二乘曲面拟合法。若最小 18
二乘的法矩阵病态，则易造成求参结果失真、不可
（ 11）
由（ 4） ，（ 5） 式确定岭参数 k，则 7 个参数改正
数的岭估计解为
王明柱等： 基于岭估计的概率积分法预计参数的求取
2012 年第 2 期
^Δ Z7 ×1
=（
N7 ×7
+
kI）
A P l －1 T 7×n n×n n×1
（ 12）
即求出第 1 次的参数改正数 Δ Z^1 ，则所求参
Abstract： Divergence problem will often result in parameters distortion in solve prediction parameters of probability integration method with least square method. Morbidity of normal matrix is one of main divergence reasons. Ridge estimate was applied in solving prediction parameters，which could effectively overcome distortion of normal matrix. At the same time，Matlab was used to realize parameters＇ automatic solving. A simulation example was present to verify the model. Key words： ridge estimate； probability integration method； parameter solve； mining subsidence； L curve method