岩移沉降观测中数据缺失的一种估算方法

赖礼泉

(福建省闽东南地质大队，福建 泉州 362000)

矿区地表变形监测在施工过程中至关重要，通常为持续性地对地表进行观测，并通过建立数学预测模型对地表形变趋势预测，从而科学施工，避免事故[1-5]。但在实际变形监测过程中，通常会出现观测数据缺失的问题，主要是由于观测点处于裂缝或塌陷区、被积水淹没、遭到人为破坏[6]。当变形监测数据出现缺失时通常利用填补法进行弥补，即利用某种原则构造填补值代替缺失数据，补全数据序列[7]。常用的数据填补方法为回归分析法，计算过程简单，但由于选用的影响因子具有不可测性，导致回归分析法受到限制。本文引入最大期望算法(Expectation Maximization Model,EM)，通过观测数据的边际分布对缺失数据进行几大似然估计，当E步和M步的迭代之间参数变化小于给定阈值时得到缺失数据填补值。

1 模型的建立

1.1 回归分析模型的建立

假设有一组已知的观测数据x1,x2,…,xk，已知观测数据不同的观测期数记为xk1,xk2,…,xki，由此构建的回归模型为：

yi=β0+β1x1i+β2x2i+…+βkxki+εi

(1)

式中，yi为y第i次的观测值；β0，β1，…，βk为未知回归参数；εi为随机误差，εi～N(0,σ2)，且协方差σεi·εi-τ=0(τ≠0)。

所以误差方程为：

vi=β0+β1x1i+β2x2i+…+βkxki-yi

(2)

化为矩阵形式为：

V=Aβ-Y

(3)

式中，V=(v1,v2,…,vn)T；β=(β0,β1,…,βn)T；Y=(y1,y2,…,yn)T；

从而得到参数的最小二乘解为：

(4)

解得回归方程为：

(5)

进行岩移观测，当出现数据缺失时，通过上述过程利用缺失数据附近的点建立回归模型，最终可以得到缺失数据的近似值。

1.2 最大期望算法的建立

最大期望算法是由Dempster提出的，是利用不完全数据求取几大似然估计的算法[8-13]。假设x(i)为不完全观测数据，z(i)为数据的缺失部分，x(i)与z(i)共同构成完整的数据序列，设未知参数θ=[u,σ]T，其中，u表示z(i)的均值，σ表示z(i)的方差，两者均未知。

EM算法流程如下：

随机初始化未知参数θ；循环重复E步和M步直到收敛。

E(Expectation)步：计算期望，利用对隐藏变量的现有估计值，计算其最大似然估计值；对于每一个i，根据上一次迭代的模型参数来计算出z(i)的后验概率，作为z(i)的现估计值，公式如下：

Qi(z(i))=p(z(i)|x(i);θ)

(6)

M(Maximization)步：最大化，最大化在E步上求得的最大似然值来计算参数的值，M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行。求使Q函数获得极大时的未知参数θ取值，将似然函数最大化以获得新的未知参数θ，公式如下：

(7)

Qi(z(i))求出来代入到θ，θ求出来又反代回Qi(z(i))，如此不断的迭代，就可以得到使似然函数最大化的参数θ了。

其实，E步就是求给定X下的条件期望，即后验期望，对Jensen不等式中小的那一端进行放大，使其等于大的那一端，这是一次放大；M步最大化联合分布，通过0梯度，拉格朗日法等方法求极值点，又是一次放大。似然函数是有界的，只要M步中的0梯度点是极大值点，一直放大下去就能找到最终所求。

具体计算过程如下：

如果θ(i)就是第i+1次迭代开始时的参数估计值，则第i+1次迭代时的E步与M步如下：

E步：Qi(z(i))=p(z(i)|x(i);θ(i))

(8)

M步：θ(i+1)=argmax∑i∑z(i)Qi(z(i))

(9)

通过上述步骤完成了一次迭代过程θ(i)→θ(i+1)，继续进行迭代，过程如下：

E步：Qi(z(i))=p(z(i)|x(i);θ(i+1))

(10)

M步：θ(i+2)=argmax∑i∑z(i)Qi(z(i))

(11)

通过上述步骤完成了一次迭代过程θ(i+1)→θ(i+2)，继续进行迭代，过程如下：

E步：Qi(z(i))=p(z(i)|x(i);θ(i+2))

(12)

M步：θ(i+3)=argmax∑i∑z(i)Qi(z(i))

(13)

通过上述步骤完成了一次迭代过程θ(i+2)→θ(i+3)，将E步和M步不断进行迭代，直到‖θ(i+1)-θ(i)‖充分小，说明残差已经满足要求，即停止迭代。

2 算例分析

选取某矿区730采区测量原始数据，730采区位于矿井东南部，-980 m延深水平的南部，南至铁路保护煤柱(矿井南部边界)，东南至矿井东南部边界，东至南部轨道集中巷及其东南方向延长线，北至-980 m水平一节轨道及二节胶带暗斜井，西至支断层。采区南北长1 469～2 009 m，东西宽318～1 518 m，面积约为2.47 km2。测线布置重点区域采用点间距25 m，根据现场条件做一定调整，测线布置图如图1所示，其中，加粗线为陀螺边。

图1 测线布置图

为了比较回归分析法和EM法对缺失的下沉数据预测效果，在完备数据中先剔除部分观测值，选取A11点2015年12月10日～2016年11月24日共10期观测数据进行试验，图2为A11点利用回归分析法和EM法对10期数据的估算结果与原始观测数据对比结果，其中，纵轴为沉降值，单位为m。

图2 A11点利用回归分析法和EM法的估算结果与原始观测数据对比结果

由图2可知，回归分析法和EM法对缺失数据的估算结果都在原始观测数据附近，说明这两种方法都可以对缺失沉降观测数据进行估算；EM法缺失数据估算结果走势与观测数据更加接近，表明EM法对缺失数据估算结果优于回归分析法。为了更加直观的对比两种估算方法的精度，把残差的绝对值和残差平方和作为对比指标。残差为实际观测数据与估算填补数据的差值。实际观测数据与两种方法估算填补数据残差绝对值对比结果如表1所示。

表1 实际观测数据与两种方法估算填补数据残差绝对值对比结果/m

由表1可知，回归分析法估算填补数据残差绝对值在0.643～0.986 m之间，其中，最大为第12期数据的0.986 m，最小为第14期数据的0.643 m。EM法估算填补数据残差绝对值在0.103～0.360 m之间，其中，最大值为第14期数据的0.360 m，最小为第20期数据的0.103 m。从估算填补数据的残差绝对值来看EM法效果更好。

表2 两种方法估算填补数据中误差对比结果/m

由表2可知，回归分析法填补数据中误差为0.851 m，EM法填补数据中误差为0.219 m，EM法的中误差仅为回归分析法的25.73%，说明算法中EM法的迭代计算能够有效获取原始数据的规律，更为准确地补充缺失的数据；实际计算时应根据需求的精度，调整阈值‖θ(i+1)-θ(i)‖，在达到精度时停止式(13)的计算。

3 结 论

针对岩移沉降观测中数据缺失的问题，本文介绍了回归分析法和EM法对缺失数据估算填补，利用某矿区岩移沉降观测数据进行实验，对实验结果对比分析发现：

(1)EM法估算数据走势更接近实际观测数据，回归分析法虽然误差较大，但与实测数据偏离不大；

(2)虽然EM法精度较高，但计算复杂，需要重复迭代，时间较长，若实际工程中对岩移缺失补充精度要求不高，也可以使用回归分析法计算其概略值；

(3)回归分析法填补数据中误差为0.851 m，EM法填补数据中误差为0.219 m，可知EM多次迭代后精度相对回归分析法有较大幅度提升，实际工程中若对缺失数据精度要求高，可以考虑增加迭代次数，牺牲时间以获得更好的岩移沉降观测缺失数据补充；

(4)需要注意的是，随机初始化未知参数θ的选取对最终迭代结果的准确性有着决定性的影响，若θ选取错误，最终迭代结果可能会出现不收敛或精度偏低的情况，如何准确地选取有利于迭代计算的θ值是下一步的重点研究内容。