加权递推最小二乘在绝对重力仪数据处理中的应用

罗棋，庞聪，李查玮

（中国地震局地震研究所（地震大地测量重点实验室），湖北武汉430071）

加权递推最小二乘在绝对重力仪数据处理中的应用

罗棋，庞聪，李查玮

（中国地震局地震研究所（地震大地测量重点实验室），湖北武汉430071）

针对绝对重力仪高精度的需求，在抗差参数估计的基础上提出了加权递推算法及其改进算法，加权递推和改进算法改良抗差权重因子并且避免求复杂多维的逆矩阵。通过仿真表明，该算法对含有白噪声或者异常数据都起到很好的仰制作用并提高了重力测量精度。

绝对重力仪；抗差估计；加权递推；最小二乘

绝对重力测量在计量学，地球物理和大地测量有着重要的作用，其精度的提高对地震活动、地质勘测都有不可估量的作用[1]。一般测量绝对重力加速度g的方法是基于落体每下落半个波长距离产生一个干涉条纹的规律，获得物体下落的位移信息，对比铷原子钟给出的时间，得到时间距离对，拟合得到重力加速度值[2]。

目前普遍采用最小二乘和非线性拟合时间距离对[3]，与非线性相比，最小二乘原理简单，易于理解和掌握，且在编程中易于实现，但是最小二乘算法适用于近似服从正态分布的测量数据，参数拟合结果对非正态分布的观测粗差极为敏感[4]，没有考虑观测误差的不同对算法解的影响[5]。用双采样过零检测，时刻数据将越来越密集，时间间隔越来越短，测量的误差影响将会放大，因此需要对最小二乘进行改进。文献[4]提出了抗差估计，但是涉及到求解复杂多维的逆矩阵。基于上述问题，为了进一步提供算法的可行性和精度要求，提出加权递推形式，避免了复杂逆矩阵的求解。文中基于文献[4]中的迭代算法，采用文献[5]的加权阵，对拟合的时间距离对进行加权处理，拟合误差方差越大的赋予的权重越小，误差方差小的赋予的权重大，这也充分体现了加权递推最小二乘算法的优越性[5]。

1 加权递推最小二乘算法

1.1 加权递推最小二乘

绝对重力仪研究的对象是做自由落体运动的落体，针对落体运动，绝对重力仪数据符合二次多项式曲线（包含重力梯度和光速有限性请参考文献[6]中的算法，在拟合公式中增加正弦项和余弦项以仰制激光器调制频率对输出的影响参考文献[4]，将其整理后依旧可以带入最小二乘拟合参数，并不影响算法的通用性），通常假设时间距离对数据为(ti,si)(i=1,2,…m)，最小二乘的基本思想即选择一个n次多项式s(t)用一种使均方误差极小的方式来逼近(ti,si)，其中n远小于m。

写成向量模型，模型偏差为：

其中ε为偏差，。

有m个等式，将其合成矩阵向量形式：

选取适当的对角加权矩阵W，得到加权非递推最小二乘解：

初始化：P(0)=aIm，a为很大的实数，Im为m×m单位矩阵，（0）=0 。

图1给出了对应的加权递推算法计算框图，绝对重力仪数据处理系统每读取一个时间距离对(ti,si)，根据式（2）算出偏差ε(ti)，取权值带入加权递推最小二乘中，通过迭代，拟合出高精度的参数估计。对式（4）取加权矩阵W(m)=Im，则上述递推公式则为递推最小二乘。

图1 加权递推算法计算框图

1.2 加权递推算法改进

在上述递推公式中，每次读取到新的时间距离对数据，立刻参与计算得到参数，更新的参数取代旧参数参与下一步的拟合，即可以立刻判断参数的优异程度，但是数据精度与初始化参数有关，并且需要大量拟合数据来消去初始化参数带来的影响。实际运用中，先存取所有的时间距离对，再取出合适数据来参与计算，很明显这种分段式处理数据，其精度将得到提高。文中将结合这种方式，对加权递推方式做出改进。

将公式（5）做变形处理

上式中s0,v0,g0分别为初始位移、初始速度、当地重力加速度，在实际中只需要求出重力加速度参数g值，因此对式（8）做变形写成矩阵形式：

式中δi=w(i)为权值 ，当δi=1 时，上式（9）（10）为改进最小二乘递推。明显看出改进之后的算法避免了求逆矩阵，由于时间距离对的选取可以是非连续的，在数据选取和处理上是非常容易实现的。

2 仿真实验

为验证算法在绝对重力仪数据处理中的可行性，文中基于MATLAB软件分别针对最小二乘，加权最小二乘递推，改进最小二乘递推和改进加权最小二乘递推算法这4种情况进行了模拟仿真。对文献[7]中的扫频正弦信号模拟干涉条纹进行双采样提取时间距离对，往信号中添加18%的幅值噪声。

图2给出了这4种算法的g值仿真结果，其中直线代表真实值g=9.8 m/s2，曲线代表估计值，用上述4种算法及其抗差最小二乘得到的数据处理结果见表1。

表1 4种算法及抗差加权仿真结果

如图2和表1可以看出，加权后的算法精度和准度比未加权的高出两个量级，同时改进后的算法比未改进的算法更快速的收敛于真实值，即改进算法只需要相对少量数据便可得到高精度的估计值。在干扰严重时，无论是抗差还是加权递推都有效地提高了精度，当只加入幅值白噪声时，抗差效果比最小二乘提高不大，但是加权递推依旧提高了精度，这充分的说明了本文提出的算法能够满足准确性。

图2 4种算法参数估计收敛图

3 结束语

文中借鉴文献[4]中抗差的思想和迭代流程，就算法的改进提出了加权递推最小二乘算法，采用平差准则作为加权阵，经过改进算得到的估计g值不仅具有抗差参数估计的优越性，还具有运算简单方便，易于实现、移植的特征。结合计算机运算特点，改进后的算法通过仿真实验验证，具有快速收敛特性，不依赖初始值设定的特点，经过多次实验，改进算法非常稳定，具有工程应用的可行性。该算法还可以在参数细节上进行修改，以达到实际绝对重力仪数据处理中最佳的估计值。

[1]Svitlov S.Frequency-domain analysis of absolute gravimeters[J].Metrologia，2012，49（6）:706-726.

[2]吴琼.高精度绝对重力仪关键技术研究[D].北京：中国地震局地球物理研究所，2011.

[3]Svitlov S，Maslyk P，Rothleithner C，et al.Comparison of Three Digital Fringe Signal Processing Methods in a Ballistic Free-Fall Absolute Gravimeter[J].Metrologia，2010，47（6）:677-689.

[4]胡明，张为民.抗差估计在自由落体式绝对重力仪中的应用[J].大地测量与地球动力学，2016，9（9）:822-825.

[5]刘谢进.递推加权最小二乘算法的研究[J].系统仿真学报，2007，21（14）:4248-4250.

[6]吴书清，吉望西，刘达伦，等.绝对重力仪数值计算及不确定度评定方法[J].计量学报，2009，30（3）：212—215.

[7]杨厚丽，邹彤，郭唐永，等.复外差法处理绝对重力仪数字干涉信号[J].大地测量与地球动力学，2015，4（2）:346-348.

[8]周姗姗，柴金广，李丹.实时递推的最小二乘预测跟踪算法[J].电子设计工程，2011，19（10）：53-55.

[9]强明辉，张京娥.基于MATLAB的递推最小二乘法识别与仿真[J].自动化与仪器仪表，2008（6）:4-5，39.

[10]贾沛璋.误差分析与数据处理[M].北京：国防工业出版社，1992.

[11]秦延，陈宗海，李衍杰.递推最小二乘算法的补充证明[J].系统仿真学报，2004，16（10）:2159-2160.

[12]Orlob M，Braun A.Impact Estimation and Filtering ofDisturbancesin FG5 AbsoluteGravimeter Observations[J].International Journal of Geosciences，2013，4（2）:302-308.

[13]邓自立，王欣，高媛.建模与估计[M].北京：科学出版社，2007.

[14]邵建新.最小二乘法线性拟合中参数的确定问题[J].大学物理，2003，22（1）.

[15]卓金武.MATLAB在数学建模中的应用[M].北京：北京航空航天大学出版社，2014.

Application of weighted recursive least square to data processing of absolute gravimeter

LUO Qi，PANG Cong，LI Cha⁃wei
（Institute of Seismology，CEA，Wuhan430071，China）

In view of the high precision demand of absolute gravimeter，the weighted recursive algorithm and its improved algorithm are proposed based on the robust parameter estimation，weighted recursive algorithm and improved robust weighting factor and inverse matrix for complex multidimensional.Simulation results show that the proposed algorithm has a good effect on the white noise or abnormal data and improves the accuracy of gravity measurement.

absolute gravimeter；robust estimation；weighted recursive；least square

P315

A

1674-6236（2017）22-0065-04

2016-09-21稿件编号：201609188

罗棋（1990—），男，湖北天门人，硕士。研究方向：绝对重力仪落体控制。