基于改进的约束总体最小二乘的地磁测量误差补偿方法与流程

本发明涉及地磁测量误差补偿领域，特别是涉及基于改进的约束总体最小二乘的地磁测量误差补偿方法。

背景技术：

现有地磁测量误差补偿方法中，多采用约束总体最小二乘法(CTLS)等模型对测量误差进行补偿。CTLS算法是最小二乘法的扩展，适用于磁场测量中误差矩阵和误差测量值相关的情况，其前提条件为标定后的测量值分布在以当地磁场强度为半径的圆球上，利用CTLS算法获取误差转移矩阵的最优估计值，实现对地磁测量误差的补偿。

但是，当随机噪声存在时，该球体会变为表面不规则的球体，无法满足CTLS方法的前提条件，使得最优估计值失去参考价值。去除随机噪声的常见方法有卡尔曼滤波、小波变换等，但均无法满足磁场测量非线性、非周期性的需求。

技术实现要素：

发明目的：本发明的目的是提供一种能够解决现有技术中存在的缺陷的基于改进的约束总体最小二乘的地磁测量误差补偿方法。

技术方案：为达到此目的，本发明采用以下技术方案：

本发明所述的基于改进的约束总体最小二乘的地磁测量误差补偿方法，包括如下步骤：

S1：利用三轴磁力仪试验平台釆集多组地磁数据序列{x₀(t)}、{y₀(t)}和{z₀(t)}，其中{x₀(t)}为x轴方向上的磁场强度，{y₀(t)}为y轴方向上的磁场强度，{z₀(t)}为z轴方向上的磁场强度；

S2：对{x₀(t)}、{y₀(t)}和{z₀(t)}进行重构，得到重构信号{x₁(t)}、{y₁(t)}和{z₁(t)}；

S3：针对步骤S2得到的重构信号{x₁(t)}、{y₁(t)}和{z₁(t)}，消除各种误差后获得地磁测量数据的最优估计值。

进一步，所述步骤S2包括以下步骤：

S2.1：对{x₀(t)}进行模态分解，包括以下步骤：

a：令xk＝1，xp＝1；

b：对{x₀(t)}进行模态分解，寻找信号{x₀(t)}中所有极大值点，利用三次样条插 值拟合得到极大值包络线m_xk+(t)，寻找信号{x₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_xk-(t)，计算均值包络线为

c：将信号{x₀(t)}减去m_xk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第xp个本征信号函数，记为IMFX_xp(t)，然后进行步骤d；如果不满足，则令xk＝xk+1，xp＝1，然后返回步骤b；

d：将信号{x₀(t)}减去IMFX_xp(t)，得到一个去掉高频信号的新信号r_xp(t)；

e：判断新信号r_xp(t)是否还能再进行分解：如果能，则将r_xp(t)赋值给{x₀(t)}，令xp＝xp+1，xk＝1，然后返回步骤b；如果不能，则记录分解的总次数为nx，最终得到的不能分解的新信号为r_x(t)，信号{x₀(t)}如式(1)所示：

S2.2：确定信号{x₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(2)计算IMFX₁(t),…,IMFX_nx(t)中相邻两个本征信号函数之间的均方差CMSE，如xq＝u时，CMSE发生突变，u即为噪声与信号起主导作用的模态分界点；

S2.3：利用小波软阈值法去噪；

对第1个到第u个IMFX信号利用小波软阈值法去噪：利用式(3)计算第i个IMFX信号的阈值xt_i，1≤i≤u，对第i个IMFX信号进行重构，得到IMFX_i'，IMFX_i'中第ii个点为IMFX_i'(ii)，如式(4)所示；

式(3)中的media(·)为中值；

S2.4：对信号{x₀(t)}进行重构，重构信号为

S2.5：对{y₀(t)}进行模态分解，包括以下步骤：

S2.51：令yk＝1，yp＝1；

S2.52：对{y₀(t)}进行模态分解，寻找信号{y₀(t)}中所有极大值点，利用三次样条插值拟合得到极大值包络线m_yk+(t)，寻找信号{y₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_yk-(t)，计算均值包络线为

S2.53：将信号{y₀(t)}减去m_yk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第yp个本征信号函数，记为IMFY_yp(t)，然后进行步骤S2.54；如果不满足，则令yk＝yk+1，yp＝1，然后返回步骤S2.52；

S2.54：将信号{y₀(t)}减去IMFY_yp(t)，得到一个去掉高频信号的新信号r_yp(t)；

S2.55：判断新信号r_yp(t)是否还能再进行分解：如果能，则将r_yp(t)赋值给{y₀(t)}，令yp＝yp+1，yk＝1，然后返回步骤S2.52；如果不能，则记录分解的总次数为ny，最终得到的不能分解的新信号为r_y(t)，信号{y₀(t)}如式(5)所示：

S2.6：确定{y₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(6)计算IMFY₁(t),…,IMFY_ny(t)中相邻两个本征信号函数之间的均方差CMSE，如yq＝v时，CMSE发生突变，v即为噪声与信号起主导作用的模态分界点。

S2.7：利用小波软阈值法去噪；

对第1个到第v个IMFY信号利用小波软阈值法去噪：利用式(7)计算第j个IMFY信号的阈值yt_j，1≤j≤v，对第j个IMFY信号进行重构，得到IMFY_j'，IMFY_j'中第jj个点为IMFY_j'(jj)，如式(8)所示；

式(7)中的media(·)为中值；

S2.8：对信号{y₀(t)}进行重构，重构信号为

S2.9：对{z₀(t)}进行模态分解，包括以下步骤：

S2.91：令zk＝1，zp＝1；

S2.92：对{z₀(t)}进行模态分解，寻找信号{z₀(t)}中所有极大值点，利用三次样条插值拟合得到极大值包络线m_zk+(t)，寻找信号{z₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_zk-(t)，计算均值包络线为

S2.93：将信号{z₀(t)}减去m_zk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第zp个本征信号函数，记为IMFZ_zp(t)，然后进行步骤S2.94；如果不满足，则令zk＝zk+1，zp＝1，然后返回步骤S2.92；

S2.94：将信号{z₀(t)}减去IMFZ_zp(t)，得到一个去掉高频信号的新信号r_zp(t)；

S2.95：判断新信号r_zp(t)是否还能再进行分解：如果能，则将r_zp(t)赋值给{z₀(t)}，令zp＝zp+1，zk＝1，然后返回步骤S2.92；如果不能，则记录分解的总次数为nz，最终得到的不能分解的新信号为r_z(t)，信号{z₀(t)}如式(9)所示：

S2.10：确定信号{z₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(10)计算IMFZ₁(t),…,IMFZ_nz(t)中相邻两个本征信号函数之间的均方差CMSE，如zq＝w时，CMSE发生突变，w即为噪声与信号起主导作用的模态分界点；

S2.11：利用小波软阈值法去噪；

对第1个到第w个IMFZ信号利用小波软阈值法去噪：利用式(11)计算第k个IMFZ信号的阈值zt_j，1≤k≤w，对第k个IMFZ信号进行重构，得到IMFZ'_k，IMFZ'_k中第kk个点为IMFZ'_k(kk)，如式(12)所示；

式(11)中的media(·)为中值；

S2.12：对信号{z₀(t)}进行重构，重构信号为

进一步，所述步骤S3中，采用约束总体最小二乘法对重构信号{x₁(t)}、{y₁(t)}和{z₁(t)}进行误差补偿，消除各种误差项。

有益效果：本发明采用改进的约束总体最小二乘的方法，能够更加精确地补偿三轴磁力仪的测量误差，有效改善了约束总体最小二乘法在随机噪声存在时无法正常工作的缺点，提高了磁场测量的精度；采用改进的经验模态分解法，极大地保留了各个本征信号的有效部分，有效减少了滤波带来的信号损失，并使得预处理后的数据满足约束总体最小二乘法的使用条件。

附图说明

图1为本发明具体实施方式的无随机噪声时地磁场测量值满足的约束分布；

图2为本发明具体实施方式的存在随机噪声时地磁场测量值满足的约束分布；

图3为本发明具体实施方式的方法流程图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步的介绍。

如图1所示，无随机噪声存在时，地磁场的测量值分布在以当地地磁场强度为半径的球体表面上，满足约束总体最小二乘法的使用前提。

如图2所示，存在随机噪声时，地磁场的测量值分布在以当地地磁场强度为半径的不光滑球体表面上，无法满足约束总体最小二乘法的使用前提。

如图3所示，本具体实施方式公开了一种基于改进的约束总体最小二乘的地磁 测量误差补偿方法，包括如下步骤：

S1：利用三轴磁力仪试验平台釆集多组地磁数据序列{x₀(t)}、{y₀(t)}和{z₀(t)}，其中{x₀(t)}为x轴方向上的磁场强度，{y₀(t)}为y轴方向上的磁场强度，{z₀(t)}为z轴方向上的磁场强度；

S2：对{x₀(t)}、{y₀(t)}和{z₀(t)}进行重构，得到重构信号{x₁(t)}、{y₁(t)}和{z₁(t)}；

S3：针对步骤S2得到的重构信号{x₁(t)}、{y₁(t)}和{z₁(t)}，消除各种误差后获得地磁测量数据的最优估计值。

其中，步骤S2包括以下步骤：

S2.1：对{x₀(t)}进行模态分解，包括以下步骤：

a：令xk＝1，xp＝1；

b：对{x₀(t)}进行模态分解，寻找信号{x₀(t)}中所有极大值点，利用三次样条插值拟合得到极大值包络线m_xk+(t)，寻找信号{x₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_xk-(t)，计算均值包络线为

c：将信号{x₀(t)}减去m_xk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第xp个本征信号函数，记为IMFX_xp(t)，然后进行步骤d；如果不满足，则令xk＝xk+1，xp＝1，然后返回步骤b；

d：将信号{x₀(t)}减去IMFX_xp(t)，得到一个去掉高频信号的新信号r_xp(t)；

e：判断新信号r_xp(t)是否还能再进行分解：如果能，则将r_xp(t)赋值给{x₀(t)}，令xp＝xp+1，xk＝1，然后返回步骤b；如果不能，则记录分解的总次数为nx，最终得到的不能分解的新信号为r_x(t)，信号{x₀(t)}如式(1)所示：

S2.2：确定信号{x₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(2)计算IMFX₁(t),…,IMFX_nx(t)中相邻两个本征信号函数之间的均方差CMSE，如xq＝u时，CMSE发生突变，u即为噪声与信号起主导作用的模态分界点；

S2.3：利用小波软阈值法去噪；

对第1个到第u个IMFX信号利用小波软阈值法去噪：利用式(3)计算第i个IMFX信号的阈值xt_i，1≤i≤u，对第i个IMFX信号进行重构，得到IMFX_i'，IMFX_i'中第ii个点为IMFX_i'(ii)，如式(4)所示；

式(3)中的media(·)为中值；

S2.4：对信号{x₀(t)}进行重构，重构信号为

S2.5：对{y₀(t)}进行模态分解，包括以下步骤：

S2.51：令yk＝1，yp＝1；

S2.52：对{y₀(t)}进行模态分解，寻找信号{y₀(t)}中所有极大值点，利用三次样条插值拟合得到极大值包络线m_yk+(t)，寻找信号{y₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_yk-(t)，计算均值包络线为

S2.53：将信号{y₀(t)}减去m_yk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第yp个本征信号函数，记为IMFY_yp(t)，然后进行步骤S2.54；如果不满足，则令yk＝yk+1，yp＝1，然后返回步骤S2.52；

S2.54：将信号{y₀(t)}减去IMFY_yp(t)，得到一个去掉高频信号的新信号r_yp(t)；

S2.55：判断新信号r_yp(t)是否还能再进行分解：如果能，则将r_yp(t)赋值给{y₀(t)}，令yp＝yp+1，yk＝1，然后返回步骤S2.52；如果不能，则记录分解的总次数为ny，最终得到的不能分解的新信号为r_y(t)，信号{y₀(t)}如式(5)所示：

S2.6：确定{y₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(6)计算IMFY₁(t),…,IMFY_ny(t)中相邻两个本征信号函数之间的均方差CMSE，如yq＝v时，CMSE发生突变，v即为噪声与信号起主导作用的模态分界点。

S2.7：利用小波软阈值法去噪；

对第1个到第v个IMFY信号利用小波软阈值法去噪：利用式(7)计算第j个IMFY信号的阈值yt_j，1≤j≤v，对第j个IMFY信号进行重构，得到IMFY_j'，IMFY_j'中第jj个点为IMFY_j'(jj)，如式(8)所示；

式(7)中的media(·)为中值；

S2.8：对信号{y₀(t)}进行重构，重构信号为

S2.9：对{z₀(t)}进行模态分解，包括以下步骤：

S2.91：令zk＝1，zp＝1；

S2.92：对{z₀(t)}进行模态分解，寻找信号{z₀(t)}中所有极大值点，利用三次样条插值拟合得到极大值包络线m_zk+(t)，寻找信号{z₀(t)}中所有极小值点，利用三次样条插值拟合得到极小值包络线m_zk-(t)，计算均值包络线为

S2.93：将信号{z₀(t)}减去m_zk(t)得到新信号判断新信号是否满足本征信号函数的条件：如果满足，则令为第zp个本征信号函数，记为IMFZ_zp(t)，然后进行步骤S2.94；如果不满足，则令zk＝zk+1，zp＝1，然后返回步骤S2.92；

S2.94：将信号{z₀(t)}减去IMFZ_zp(t)，得到一个去掉高频信号的新信号r_zp(t)；

S2.95：判断新信号r_zp(t)是否还能再进行分解：如果能，则将r_zp(t)赋值给{z₀(t)}，令zp＝zp+1，zk＝1，然后返回步骤S2.92；如果不能，则记录分解的总次数为nz，最终得到的不能分解的新信号为r_z(t)，信号{z₀(t)}如式(9)所示：

S2.10：确定信号{z₀(t)}中噪声与信号起主导作用的模态分界点；

按照式(10)计算IMFZ₁(t),…,IMFZ_nz(t)中相邻两个本征信号函数之间的均方差CMSE，如zq＝w时，CMSE发生突变，w即为噪声与信号起主导作用的模态分界点；

S2.11：利用小波软阈值法去噪；

对第1个到第w个IMFZ信号利用小波软阈值法去噪：利用式(11)计算第k个IMFZ信号的阈值zt_j，1≤k≤w，对第k个IMFZ信号进行重构，得到IMFZ'_k，IMFZ'_k中第kk个点为IMFZ'_k(kk)，如式(12)所示；

式(11)中的media(·)为中值；

S2.12：对信号{z₀(t)}进行重构，重构信号为

。