本發(fā)明涉及地磁測(cè)量誤差補(bǔ)償領(lǐng)域，特別是涉及基于改進(jìn)的約束總體最小二乘的地磁測(cè)量誤差補(bǔ)償方法。

背景技術(shù)：

現(xiàn)有地磁測(cè)量誤差補(bǔ)償方法中，多采用約束總體最小二乘法(CTLS)等模型對(duì)測(cè)量誤差進(jìn)行補(bǔ)償。CTLS算法是最小二乘法的擴(kuò)展，適用于磁場(chǎng)測(cè)量中誤差矩陣和誤差測(cè)量值相關(guān)的情況，其前提條件為標(biāo)定后的測(cè)量值分布在以當(dāng)?shù)卮艌?chǎng)強(qiáng)度為半徑的圓球上，利用CTLS算法獲取誤差轉(zhuǎn)移矩陣的最優(yōu)估計(jì)值，實(shí)現(xiàn)對(duì)地磁測(cè)量誤差的補(bǔ)償。

但是，當(dāng)隨機(jī)噪聲存在時(shí)，該球體會(huì)變?yōu)楸砻娌灰?guī)則的球體，無法滿足CTLS方法的前提條件，使得最優(yōu)估計(jì)值失去參考價(jià)值。去除隨機(jī)噪聲的常見方法有卡爾曼濾波、小波變換等，但均無法滿足磁場(chǎng)測(cè)量非線性、非周期性的需求。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：本發(fā)明的目的是提供一種能夠解決現(xiàn)有技術(shù)中存在的缺陷的基于改進(jìn)的約束總體最小二乘的地磁測(cè)量誤差補(bǔ)償方法。

技術(shù)方案：為達(dá)到此目的，本發(fā)明采用以下技術(shù)方案：

本發(fā)明所述的基于改進(jìn)的約束總體最小二乘的地磁測(cè)量誤差補(bǔ)償方法，包括如下步驟：

S1：利用三軸磁力儀試驗(yàn)平臺(tái)釆集多組地磁數(shù)據(jù)序列{x₀(t)}、{y₀(t)}和{z₀(t)}，其中{x₀(t)}為x軸方向上的磁場(chǎng)強(qiáng)度，{y₀(t)}為y軸方向上的磁場(chǎng)強(qiáng)度，{z₀(t)}為z軸方向上的磁場(chǎng)強(qiáng)度；

S2：對(duì){x₀(t)}、{y₀(t)}和{z₀(t)}進(jìn)行重構(gòu)，得到重構(gòu)信號(hào){x₁(t)}、{y₁(t)}和{z₁(t)}；

S3：針對(duì)步驟S2得到的重構(gòu)信號(hào){x₁(t)}、{y₁(t)}和{z₁(t)}，消除各種誤差后獲得地磁測(cè)量數(shù)據(jù)的最優(yōu)估計(jì)值。

進(jìn)一步，所述步驟S2包括以下步驟：

S2.1：對(duì){x₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

a：令xk＝1，xp＝1；

b：對(duì){x₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){x₀(t)}中所有極大值點(diǎn)，利用三次樣條插 值擬合得到極大值包絡(luò)線m_xk+(t)，尋找信號(hào){x₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_xk-(t)，計(jì)算均值包絡(luò)線為

c：將信號(hào){x₀(t)}減去m_xk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第xp個(gè)本征信號(hào)函數(shù)，記為IMFX_xp(t)，然后進(jìn)行步驟d；如果不滿足，則令xk＝xk+1，xp＝1，然后返回步驟b；

d：將信號(hào){x₀(t)}減去IMFX_xp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_xp(t)；

e：判斷新信號(hào)r_xp(t)是否還能再進(jìn)行分解：如果能，則將r_xp(t)賦值給{x₀(t)}，令xp＝xp+1，xk＝1，然后返回步驟b；如果不能，則記錄分解的總次數(shù)為nx，最終得到的不能分解的新信號(hào)為r_x(t)，信號(hào){x₀(t)}如式(1)所示：

S2.2：確定信號(hào){x₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(2)計(jì)算IMFX₁(t),…,IMFX_nx(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如xq＝u時(shí)，CMSE發(fā)生突變，u即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

S2.3：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第u個(gè)IMFX信號(hào)利用小波軟閾值法去噪：利用式(3)計(jì)算第i個(gè)IMFX信號(hào)的閾值xt_i，1≤i≤u，對(duì)第i個(gè)IMFX信號(hào)進(jìn)行重構(gòu)，得到IMFX_i'，IMFX_i'中第ii個(gè)點(diǎn)為IMFX_i'(ii)，如式(4)所示；

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

S2.4：對(duì)信號(hào){x₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

S2.5：對(duì){y₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

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

S2.52：對(duì){y₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){y₀(t)}中所有極大值點(diǎn)，利用三次樣條插值擬合得到極大值包絡(luò)線m_yk+(t)，尋找信號(hào){y₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_yk-(t)，計(jì)算均值包絡(luò)線為

S2.53：將信號(hào){y₀(t)}減去m_yk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第yp個(gè)本征信號(hào)函數(shù)，記為IMFY_yp(t)，然后進(jìn)行步驟S2.54；如果不滿足，則令yk＝y(tǒng)k+1，yp＝1，然后返回步驟S2.52；

S2.54：將信號(hào){y₀(t)}減去IMFY_yp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_yp(t)；

S2.55：判斷新信號(hào)r_yp(t)是否還能再進(jìn)行分解：如果能，則將r_yp(t)賦值給{y₀(t)}，令yp＝y(tǒng)p+1，yk＝1，然后返回步驟S2.52；如果不能，則記錄分解的總次數(shù)為ny，最終得到的不能分解的新信號(hào)為r_y(t)，信號(hào){y₀(t)}如式(5)所示：

S2.6：確定{y₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(6)計(jì)算IMFY₁(t),…,IMFY_ny(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如yq＝v時(shí)，CMSE發(fā)生突變，v即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)。

S2.7：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第v個(gè)IMFY信號(hào)利用小波軟閾值法去噪：利用式(7)計(jì)算第j個(gè)IMFY信號(hào)的閾值yt_j，1≤j≤v，對(duì)第j個(gè)IMFY信號(hào)進(jìn)行重構(gòu)，得到IMFY_j'，IMFY_j'中第jj個(gè)點(diǎn)為IMFY_j'(jj)，如式(8)所示；

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

S2.8：對(duì)信號(hào){y₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

S2.9：對(duì){z₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

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

S2.92：對(duì){z₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){z₀(t)}中所有極大值點(diǎn)，利用三次樣條插值擬合得到極大值包絡(luò)線m_zk+(t)，尋找信號(hào){z₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_zk-(t)，計(jì)算均值包絡(luò)線為

S2.93：將信號(hào){z₀(t)}減去m_zk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第zp個(gè)本征信號(hào)函數(shù)，記為IMFZ_zp(t)，然后進(jìn)行步驟S2.94；如果不滿足，則令zk＝zk+1，zp＝1，然后返回步驟S2.92；

S2.94：將信號(hào){z₀(t)}減去IMFZ_zp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_zp(t)；

S2.95：判斷新信號(hào)r_zp(t)是否還能再進(jìn)行分解：如果能，則將r_zp(t)賦值給{z₀(t)}，令zp＝zp+1，zk＝1，然后返回步驟S2.92；如果不能，則記錄分解的總次數(shù)為nz，最終得到的不能分解的新信號(hào)為r_z(t)，信號(hào){z₀(t)}如式(9)所示：

S2.10：確定信號(hào){z₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(10)計(jì)算IMFZ₁(t),…,IMFZ_nz(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如zq＝w時(shí)，CMSE發(fā)生突變，w即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

S2.11：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第w個(gè)IMFZ信號(hào)利用小波軟閾值法去噪：利用式(11)計(jì)算第k個(gè)IMFZ信號(hào)的閾值z(mì)t_j，1≤k≤w，對(duì)第k個(gè)IMFZ信號(hào)進(jìn)行重構(gòu)，得到IMFZ'_k，IMFZ'_k中第kk個(gè)點(diǎn)為IMFZ'_k(kk)，如式(12)所示；

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

S2.12：對(duì)信號(hào){z₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

進(jìn)一步，所述步驟S3中，采用約束總體最小二乘法對(duì)重構(gòu)信號(hào){x₁(t)}、{y₁(t)}和{z₁(t)}進(jìn)行誤差補(bǔ)償，消除各種誤差項(xiàng)。

有益效果：本發(fā)明采用改進(jìn)的約束總體最小二乘的方法，能夠更加精確地補(bǔ)償三軸磁力儀的測(cè)量誤差，有效改善了約束總體最小二乘法在隨機(jī)噪聲存在時(shí)無法正常工作的缺點(diǎn)，提高了磁場(chǎng)測(cè)量的精度；采用改進(jìn)的經(jīng)驗(yàn)?zāi)B(tài)分解法，極大地保留了各個(gè)本征信號(hào)的有效部分，有效減少了濾波帶來的信號(hào)損失，并使得預(yù)處理后的數(shù)據(jù)滿足約束總體最小二乘法的使用條件。

附圖說明

圖1為本發(fā)明具體實(shí)施方式的無隨機(jī)噪聲時(shí)地磁場(chǎng)測(cè)量值滿足的約束分布；

圖2為本發(fā)明具體實(shí)施方式的存在隨機(jī)噪聲時(shí)地磁場(chǎng)測(cè)量值滿足的約束分布；

圖3為本發(fā)明具體實(shí)施方式的方法流程圖。

具體實(shí)施方式

下面結(jié)合具體實(shí)施方式對(duì)本發(fā)明的技術(shù)方案作進(jìn)一步的介紹。

如圖1所示，無隨機(jī)噪聲存在時(shí)，地磁場(chǎng)的測(cè)量值分布在以當(dāng)?shù)氐卮艌?chǎng)強(qiáng)度為半徑的球體表面上，滿足約束總體最小二乘法的使用前提。

如圖2所示，存在隨機(jī)噪聲時(shí)，地磁場(chǎng)的測(cè)量值分布在以當(dāng)?shù)氐卮艌?chǎng)強(qiáng)度為半徑的不光滑球體表面上，無法滿足約束總體最小二乘法的使用前提。

如圖3所示，本具體實(shí)施方式公開了一種基于改進(jìn)的約束總體最小二乘的地磁 測(cè)量誤差補(bǔ)償方法，包括如下步驟：

S1：利用三軸磁力儀試驗(yàn)平臺(tái)釆集多組地磁數(shù)據(jù)序列{x₀(t)}、{y₀(t)}和{z₀(t)}，其中{x₀(t)}為x軸方向上的磁場(chǎng)強(qiáng)度，{y₀(t)}為y軸方向上的磁場(chǎng)強(qiáng)度，{z₀(t)}為z軸方向上的磁場(chǎng)強(qiáng)度；

S2：對(duì){x₀(t)}、{y₀(t)}和{z₀(t)}進(jìn)行重構(gòu)，得到重構(gòu)信號(hào){x₁(t)}、{y₁(t)}和{z₁(t)}；

S3：針對(duì)步驟S2得到的重構(gòu)信號(hào){x₁(t)}、{y₁(t)}和{z₁(t)}，消除各種誤差后獲得地磁測(cè)量數(shù)據(jù)的最優(yōu)估計(jì)值。

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

S2.1：對(duì){x₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

a：令xk＝1，xp＝1；

b：對(duì){x₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){x₀(t)}中所有極大值點(diǎn)，利用三次樣條插值擬合得到極大值包絡(luò)線m_xk+(t)，尋找信號(hào){x₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_xk-(t)，計(jì)算均值包絡(luò)線為

c：將信號(hào){x₀(t)}減去m_xk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第xp個(gè)本征信號(hào)函數(shù)，記為IMFX_xp(t)，然后進(jìn)行步驟d；如果不滿足，則令xk＝xk+1，xp＝1，然后返回步驟b；

d：將信號(hào){x₀(t)}減去IMFX_xp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_xp(t)；

e：判斷新信號(hào)r_xp(t)是否還能再進(jìn)行分解：如果能，則將r_xp(t)賦值給{x₀(t)}，令xp＝xp+1，xk＝1，然后返回步驟b；如果不能，則記錄分解的總次數(shù)為nx，最終得到的不能分解的新信號(hào)為r_x(t)，信號(hào){x₀(t)}如式(1)所示：

S2.2：確定信號(hào){x₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(2)計(jì)算IMFX₁(t),…,IMFX_nx(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如xq＝u時(shí)，CMSE發(fā)生突變，u即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

S2.3：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第u個(gè)IMFX信號(hào)利用小波軟閾值法去噪：利用式(3)計(jì)算第i個(gè)IMFX信號(hào)的閾值xt_i，1≤i≤u，對(duì)第i個(gè)IMFX信號(hào)進(jìn)行重構(gòu)，得到IMFX_i'，IMFX_i'中第ii個(gè)點(diǎn)為IMFX_i'(ii)，如式(4)所示；

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

S2.4：對(duì)信號(hào){x₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

S2.5：對(duì){y₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

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

S2.52：對(duì){y₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){y₀(t)}中所有極大值點(diǎn)，利用三次樣條插值擬合得到極大值包絡(luò)線m_yk+(t)，尋找信號(hào){y₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_yk-(t)，計(jì)算均值包絡(luò)線為

S2.53：將信號(hào){y₀(t)}減去m_yk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第yp個(gè)本征信號(hào)函數(shù)，記為IMFY_yp(t)，然后進(jìn)行步驟S2.54；如果不滿足，則令yk＝y(tǒng)k+1，yp＝1，然后返回步驟S2.52；

S2.54：將信號(hào){y₀(t)}減去IMFY_yp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_yp(t)；

S2.55：判斷新信號(hào)r_yp(t)是否還能再進(jìn)行分解：如果能，則將r_yp(t)賦值給{y₀(t)}，令yp＝y(tǒng)p+1，yk＝1，然后返回步驟S2.52；如果不能，則記錄分解的總次數(shù)為ny，最終得到的不能分解的新信號(hào)為r_y(t)，信號(hào){y₀(t)}如式(5)所示：

S2.6：確定{y₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(6)計(jì)算IMFY₁(t),…,IMFY_ny(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如yq＝v時(shí)，CMSE發(fā)生突變，v即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)。

S2.7：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第v個(gè)IMFY信號(hào)利用小波軟閾值法去噪：利用式(7)計(jì)算第j個(gè)IMFY信號(hào)的閾值yt_j，1≤j≤v，對(duì)第j個(gè)IMFY信號(hào)進(jìn)行重構(gòu)，得到IMFY_j'，IMFY_j'中第jj個(gè)點(diǎn)為IMFY_j'(jj)，如式(8)所示；

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

S2.8：對(duì)信號(hào){y₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

S2.9：對(duì){z₀(t)}進(jìn)行模態(tài)分解，包括以下步驟：

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

S2.92：對(duì){z₀(t)}進(jìn)行模態(tài)分解，尋找信號(hào){z₀(t)}中所有極大值點(diǎn)，利用三次樣條插值擬合得到極大值包絡(luò)線m_zk+(t)，尋找信號(hào){z₀(t)}中所有極小值點(diǎn)，利用三次樣條插值擬合得到極小值包絡(luò)線m_zk-(t)，計(jì)算均值包絡(luò)線為

S2.93：將信號(hào){z₀(t)}減去m_zk(t)得到新信號(hào)判斷新信號(hào)是否滿足本征信號(hào)函數(shù)的條件：如果滿足，則令為第zp個(gè)本征信號(hào)函數(shù)，記為IMFZ_zp(t)，然后進(jìn)行步驟S2.94；如果不滿足，則令zk＝zk+1，zp＝1，然后返回步驟S2.92；

S2.94：將信號(hào){z₀(t)}減去IMFZ_zp(t)，得到一個(gè)去掉高頻信號(hào)的新信號(hào)r_zp(t)；

S2.95：判斷新信號(hào)r_zp(t)是否還能再進(jìn)行分解：如果能，則將r_zp(t)賦值給{z₀(t)}，令zp＝zp+1，zk＝1，然后返回步驟S2.92；如果不能，則記錄分解的總次數(shù)為nz，最終得到的不能分解的新信號(hào)為r_z(t)，信號(hào){z₀(t)}如式(9)所示：

S2.10：確定信號(hào){z₀(t)}中噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

按照式(10)計(jì)算IMFZ₁(t),…,IMFZ_nz(t)中相鄰兩個(gè)本征信號(hào)函數(shù)之間的均方差CMSE，如zq＝w時(shí)，CMSE發(fā)生突變，w即為噪聲與信號(hào)起主導(dǎo)作用的模態(tài)分界點(diǎn)；

S2.11：利用小波軟閾值法去噪；

對(duì)第1個(gè)到第w個(gè)IMFZ信號(hào)利用小波軟閾值法去噪：利用式(11)計(jì)算第k個(gè)IMFZ信號(hào)的閾值z(mì)t_j，1≤k≤w，對(duì)第k個(gè)IMFZ信號(hào)進(jìn)行重構(gòu)，得到IMFZ'_k，IMFZ'_k中第kk個(gè)點(diǎn)為IMFZ'_k(kk)，如式(12)所示；

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

S2.12：對(duì)信號(hào){z₀(t)}進(jìn)行重構(gòu)，重構(gòu)信號(hào)為

。