丁學(xué)振，李卓軒，李予國(guó),2,3*，劉浩

1 中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島 266100 2 中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100 3 青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島 266237

0 引言

磁法探測(cè)在地質(zhì)調(diào)查、資源探測(cè)、軍事、環(huán)境等領(lǐng)域有重要的應(yīng)用價(jià)值（管志寧，1997；林君等，2017）.磁法探測(cè)的發(fā)展可分為總磁場(chǎng)測(cè)量、總場(chǎng)梯度測(cè)量、矢量場(chǎng)和磁梯度張量場(chǎng)測(cè)量等3個(gè)階段（張昌達(dá),2006）.相對(duì)于磁場(chǎng)數(shù)據(jù)，磁梯度張量數(shù)據(jù)具有異常信號(hào)更敏感、空間分辨率高、降低地磁噪聲、減少對(duì)基站的需求等優(yōu)勢(shì)（Pedersen and Rasmussen,1990;Schmidt and Clark,2006;Clark,2012）.

磁梯度張量是磁場(chǎng)矢量在3個(gè)正交方向上的空間變化率（Schmidt et al.,2004）.在實(shí)際測(cè)量中，可以采用差分近似法、旋轉(zhuǎn)調(diào)制法（Tilbrook,2009）和直線測(cè)量法（Sunderland et al.,2009）測(cè)得磁梯度張量數(shù)據(jù).差分近似法易于實(shí)現(xiàn)，并且能夠獲得全張量信息，因此得到了廣泛應(yīng)用（苗紅松,2017）.目前，基于差分近似原理的磁梯度張量數(shù)據(jù)采集系統(tǒng)大致可分為基于超導(dǎo)量子干涉儀的磁梯度張量系統(tǒng)（Schmidt et al.,2004;Gamey,2008;Keenan et al.,2010;汪瀛等,2015）和基于磁通門傳感器的磁梯度張量系統(tǒng)（Koch et al.,1996;Sulzberger et al.,2006;Sui et al.,2014;Yin et al.,2014）.磁梯度張量超導(dǎo)量子干涉系統(tǒng)具有很高的靈敏度，但量程范圍較小、成本高，并且對(duì)系統(tǒng)安裝工藝和測(cè)量環(huán)境要求也較高.而基于磁通門傳感器搭建的磁梯度張量系統(tǒng)成本較低，對(duì)安裝工藝要求相對(duì)較低，利于大批量生產(chǎn)制造（李青竹等,2017）.由于磁通門傳感器制作工藝和磁梯度張量系統(tǒng)安裝精度的限制，張量測(cè)量精度會(huì)受到單個(gè)傳感器誤差（零偏、三軸非正交誤差、三軸靈敏度不一致）和多傳感器非對(duì)準(zhǔn)誤差的影響，磁梯度張量系統(tǒng)輸出誤差有時(shí)可達(dá)上千nT/m （李青竹等,2018），因此必須對(duì)其進(jìn)行校正.

目前，大部分磁梯度張量系統(tǒng)校正方法為兩步校正方法，該方法包括兩個(gè)步驟（Pang et al.,2013a,b;Yin et al.,2014,2015a;遲鋮等,2017）.第一步是對(duì)單傳感器的誤差進(jìn)行校正.該校正可分為矢量校正方法（Pang et al.,2014）和標(biāo)量校正方法（Pang et al.,2013b;Yin et al.,2014,2015a;遲鋮等,2017）.矢量校正方法以高精度磁場(chǎng)矢量為標(biāo)準(zhǔn)對(duì)磁通門傳感器進(jìn)行校正，該方法需要測(cè)量準(zhǔn)確的磁場(chǎng)矢量，對(duì)試驗(yàn)條件要求較高.標(biāo)量校正方法以高精度磁場(chǎng)總場(chǎng)為標(biāo)準(zhǔn)對(duì)磁通門傳感器進(jìn)行校正，該方法因較為方便且成本低而得到廣泛應(yīng)用.第二步是對(duì)磁梯度張量系統(tǒng)的非對(duì)準(zhǔn)誤差進(jìn)行校正.以磁梯度張量系統(tǒng)中的一個(gè)傳感器或搭載平臺(tái)坐標(biāo)系（Yin et al.,2014;Li et al.,2018）為基準(zhǔn)，將多個(gè)傳感器校正至同一坐標(biāo)系中.兩步校正方法的缺點(diǎn)在于第二步的校正結(jié)果依賴于第一步校正結(jié)果的精度，第一步校正結(jié)果的誤差直接影響著第二步校正效果.為了解決兩步校正方法存在的問題，Yin等（2015b）和Li等（2018）提出了基于最小二乘擬合方法的一步校正方法.一步校正方法將磁梯度張量系統(tǒng)的理想輸出校正至搭載平臺(tái)坐標(biāo)系中，需要令張量系統(tǒng)繞三個(gè)正交的坐標(biāo)軸旋轉(zhuǎn)采集數(shù)據(jù)，該方法依賴于高精度的三軸無磁旋轉(zhuǎn)臺(tái)且試驗(yàn)難度較大.上述兩步校正方法和一步校正方法大多基于最小二乘擬合方法估計(jì)磁梯度張量系統(tǒng)的誤差參數(shù)，而最小二乘擬合反演方法對(duì)初始參數(shù)比較敏感，為獲得較好的校正結(jié)果，需對(duì)初始參數(shù)進(jìn)行仔細(xì)調(diào)整（Pang et al.,2013c）.

為了進(jìn)一步解決一步和兩步校正方法存在的問題，本文提出了一種基于改進(jìn)差分進(jìn)化算法的磁梯度張量系統(tǒng)校正方法.差分進(jìn)化算法具有全局尋優(yōu)、不受初始模型選取限制等優(yōu)點(diǎn)，可以有效解決基于最小二乘擬合的校正方法依賴于初始參數(shù)的問題.我們將張量分量和張量不變量引入反演目標(biāo)函數(shù)中，只需一步反演即可實(shí)現(xiàn)單傳感器誤差和張量系統(tǒng)非對(duì)準(zhǔn)誤差的校正，有效避免了兩步校正方法中第一步校正精度對(duì)第二步校正結(jié)果的影響問題，提高了校正精度.

1 磁梯度張量系統(tǒng)及誤差校正

1.1 張量數(shù)據(jù)采集系統(tǒng)

磁梯度張量系統(tǒng)包括十字形結(jié)構(gòu)、正方形結(jié)構(gòu)、三角形結(jié)構(gòu)、直角四面體和正四面體結(jié)構(gòu)等5種基本結(jié)構(gòu)（劉麗敏,2012）.磁梯度張量系統(tǒng)采用差分近似的原理測(cè)量張量數(shù)據(jù)，不同結(jié)構(gòu)的張量系統(tǒng)會(huì)帶來測(cè)量誤差，其中，十字形結(jié)構(gòu)的磁梯度張量系統(tǒng)測(cè)量誤差最小.圖1為十字形磁梯度張量系統(tǒng)示意圖，它由4個(gè)三軸磁通門傳感器組成，其中1號(hào)和3號(hào)傳感器沿x軸分布，2號(hào)和4號(hào)傳感器沿y軸分布，基線距離為d.根據(jù)差分近似原理，磁梯度張量系統(tǒng)中心位置o處的張量矩陣為（Li et al.,2018;Yin et al.,2015b）

（1）

式中，Bij表示第i（i=1,2,3,4）個(gè)傳感器測(cè)得的j（j=x,y,z）方向上的磁場(chǎng)分量.在磁法測(cè)量中，通常將磁性物體產(chǎn)生的異常場(chǎng)看作為無源的靜磁場(chǎng)，異常場(chǎng)的散度和旋度都為零，故張量矩陣G中只有Bxx,Bxy,Bxz,Byy,Byz5個(gè)獨(dú)立分量（尹剛等,2016）.

由張量矩陣G的特征值可以計(jì)算得到磁梯度張量不變量，三個(gè)張量不變量的表達(dá)式如下（Clark,2012;Mu et al.,2019）：

（2）

式中，λ1,λ2,λ3為張量矩陣G的3個(gè)特征值，且滿足λ1≥λ2≥λ3,|λ1|≥|λ2|,|λ3|≥|λ2|.

圖1 十字形磁梯度張量系統(tǒng)示意圖Fig.1 Schematic diagram of the cross-shaped magnetic gradient tensor system

1.2 誤差校正

在理想情況下，三軸磁通門傳感器的三個(gè)軸保持兩兩完全正交，但是由于制作工藝的限制，傳感器普遍存在非正交誤差.建立傳感器非正交模型如圖2所示，O-XYZ為標(biāo)準(zhǔn)正交坐標(biāo)系，O-X1Y1Z1為傳感器實(shí)際坐標(biāo)系，O-X2Y2Z2為理想正交坐標(biāo)系.假設(shè)OZ1和OZ2共軸，平面Y1OZ1和平面Y2OZ2共面，則OY1和OY2之間的夾角為ψ，OX1與平面X2OY2的夾角為φ，OX1在平面XOY的投影OX′1與OX2的夾角為θ.一旦非正交角φ,θ,ψ確定后，傳感器理想正交坐標(biāo)系O-X2Y2Z2即可唯一確定.

圖2 單傳感器非正交誤差示意圖Fig.2 Schematic diagram of the non-orthogonal error for single sensor

在理想情況下，磁通門傳感器在無磁環(huán)境下輸出為零，但是由于磁芯存在剩磁或者電路部分存在零位偏移導(dǎo)致傳感器三軸輸出不為零.一般情況下，傳感器的零偏可達(dá)幾十nT.另外，由于傳感器磁敏元件制作工藝的限制，可能會(huì)導(dǎo)致三軸靈敏度不一致，從而使得三軸在相同磁場(chǎng)環(huán)境下的測(cè)量值不同.傳感器三軸實(shí)際測(cè)量值與真實(shí)值的比值稱為靈敏度標(biāo)度因子（李青竹等,2017）.假定傳感器三軸零偏為b=（bx,by,bz），三軸靈敏度標(biāo)度因子為kx,ky,kz.綜合考慮零偏、靈敏度不一致誤差和非正交誤差，建立的單傳感器誤差模型為（Li et al.,2018）：

Bc=KCB+b，

（3）

在實(shí)際情況下，傳感器的實(shí)際輸出Bc是已知的，傳感器誤差的校正就是由其實(shí)際輸出求取理想輸出B，即

B=（KC）-1（Bc-b）.

（4）

由于傳感器安裝精度的限制，構(gòu)成磁梯度張量系統(tǒng)4個(gè)矢量傳感器的測(cè)量軸指向可能不一致，張量系統(tǒng)非對(duì)準(zhǔn)誤差可通過旋轉(zhuǎn)方式進(jìn)行校正（Pang et al.,2013b;Li et al.,2018）.定義傳感器繞x軸旋轉(zhuǎn)為橫傾角α，繞y軸旋轉(zhuǎn)為俯仰角β，繞z軸旋轉(zhuǎn)為方位角γ.如圖3所示，O-X1Y1Z1表示1號(hào)傳感器的正交坐標(biāo)系，O-XiYiZi表示i（i=2,3,4）號(hào)傳感器的正交坐標(biāo)系.假定i號(hào)傳感器先繞y軸旋轉(zhuǎn)，再繞x軸旋轉(zhuǎn)，最后繞z軸旋轉(zhuǎn)，可以令i號(hào)傳感器坐標(biāo)系與1號(hào)傳感器坐標(biāo)系重合，即可實(shí)現(xiàn)非對(duì)準(zhǔn)誤差校正.

圖3 傳感器非對(duì)準(zhǔn)誤差校正示意圖Fig.3 Schematic diagram of misalignment error calibration

由式（4）可知，單傳感器誤差校正后的輸出為B，經(jīng)非對(duì)準(zhǔn)誤差校正后得到的輸出Br可表示為（遲鋮等,2017）

Br=TγTαTβB，

（5）

綜合考慮單傳感器誤差和多傳感器非對(duì)準(zhǔn)誤差，可得

Br=TγTαTβ（KC）-1（Bc-b）.

（6）

由上式可知，如果得到單傳感器誤差參數(shù)和多傳感器非對(duì)準(zhǔn)誤差參數(shù)即可實(shí)現(xiàn)張量系統(tǒng)誤差校正.

2 差分進(jìn)化算法

差分進(jìn)化算法是一種模擬自然界生物遺傳和進(jìn)化機(jī)理發(fā)展起來的基于種群的全局優(yōu)化算法，它是由Storn和Price于1995年提出的（Storn and Price,1995）.差分進(jìn)化算法具有原理簡(jiǎn)單，控制參數(shù)較少，魯棒性較好的優(yōu)點(diǎn)（Storn and Price,1995,1997），已得到較為廣泛的應(yīng)用（Das and Suganthan,2011）.差分進(jìn)化算法主要包括參數(shù)初始化、變異、交叉和選擇4個(gè)步驟.參數(shù)初始化就是在預(yù)先設(shè)定的搜索空間中隨機(jī)產(chǎn)生NP個(gè)維度為P的參數(shù)向量，在該步驟中隨機(jī)產(chǎn)生的參數(shù)向量應(yīng)盡可能覆蓋整個(gè)搜索空間.變異操作的目的是產(chǎn)生新的個(gè)體以增加種群多樣性，當(dāng)前種群向量稱為目標(biāo)向量，通過使用某種變異策略對(duì)目標(biāo)向量進(jìn)行變異操作產(chǎn)生新的變異向量.在交叉操作中，令目標(biāo)向量和變異向量隨機(jī)交換參數(shù)，產(chǎn)生新的試驗(yàn)向量.選擇的目的是保持種群數(shù)量不變，使用“貪婪性”選擇策略將優(yōu)良個(gè)體保留至下一代種群中.下面，結(jié)合磁梯度張量系統(tǒng)誤差校正的目的介紹差分進(jìn)化算法的實(shí)現(xiàn)步驟.

（1）參數(shù)初始化

在參數(shù)初始化這一步驟中，隨機(jī)產(chǎn)生NP個(gè)維度為P的參數(shù)向量.針對(duì)磁梯度張量系統(tǒng)誤差校正這一應(yīng)用，每個(gè)傳感器有12個(gè)誤差參數(shù)，于是參數(shù)向量的維度為48.假設(shè)當(dāng)前進(jìn)化代數(shù)為g，第m個(gè)體可以表示為

m=1,2,…,NPi=1,2,3,4

（7）

式中，φ,θ,ψ表示傳感器的非正交角（圖2），kx,ky,kz表示三軸靈敏度標(biāo)度因子，bx,by,bz表示零偏誤差，α,β,γ表示非對(duì)準(zhǔn)角（圖3），i表示傳感器序號(hào).

在參數(shù)初始化過程中，應(yīng)盡量保證隨機(jī)產(chǎn)生的參數(shù)向量能夠覆蓋整個(gè)搜索范圍.假設(shè)參數(shù)向量的取值上下界為Xmax=（x1,max,x2,max,…,xP,max）和Xmin=（x1,min,x2,min,…,xP,min），則初始化過程可表示為

（2）變異操作

（9）

式中，F(xiàn)為處于[0,2]之間變異因子（Storn and Price,1997），r1,r2,r3為[1,NP]之間隨機(jī)選擇的不同于m的整數(shù).變異策略可以表示為DE/X/Y/Z的形式，其中DE表示差分進(jìn)化算法，X表示指定的變異向量，Y表示變異過程中使用的差分向量個(gè)數(shù)，Z表示交叉方式（Das and Suganthan,2011）.

（3）交叉操作

（10）

（4）選擇操作

為了使下一代中種群數(shù)量保持不變，采用“貪婪性”策略將試驗(yàn)向量與目標(biāo)向量進(jìn)行對(duì)比，優(yōu)良個(gè)體將會(huì)保留至下一代種群中，選擇操作可以表示為

（11）

式中，f（x）為目標(biāo)函數(shù)，其值越小表明個(gè)體表現(xiàn)越好.

在傳統(tǒng)的差分進(jìn)化算法中，變異因子F和交叉概率因子CR均為固定值，且只采用DE/rand/1/bin變異策略.為了提高算法的搜索效率和精度，避免陷入局部最優(yōu)，本文使用DE/rand/1/bin、DE/target-to-best/2/bin變異策略（陳亮,2012）、時(shí)變變異因子（顏學(xué)峰等，2006）和時(shí)變交叉概率因子（王天意,2015）對(duì)差分進(jìn)化算法進(jìn)行了改進(jìn).

改進(jìn)差分進(jìn)化算法迭代過程中變異策略選取方法為

（12）

式中，rand（0,1）為處于[0,1]之間的隨機(jī)數(shù)，φ為隨迭代次數(shù)變化的閾值，可以表示為

（13）

其中，φmax=1,φmin=0.1分別為φ取值上下界（陳亮,2012），itermax和iter分別為最大迭代次數(shù)和當(dāng)前迭代次數(shù).

時(shí)變變異因子可以表示為（顏學(xué)峰等，2006）

（14）

式中，F(xiàn)0為初始變異因子.

時(shí)變交叉概率因子可以表示為（王天意,2015）

（15）

式中，CR0為初始交叉概率因子.

如果磁梯度張量系統(tǒng)不存在誤差，在勻強(qiáng)磁場(chǎng)環(huán)境下旋轉(zhuǎn)其姿態(tài)測(cè)量數(shù)據(jù)，由磁通門傳感器三軸讀數(shù)計(jì)算得到的總磁場(chǎng)為一定值，且不同傳感器相同方向的測(cè)量軸讀數(shù)一致.磁梯度張量系統(tǒng)誤差校正標(biāo)準(zhǔn)為

（16）

式中，Bxi,Byi,Bzi（i=1,2,3,4）分別為由每個(gè)個(gè)體向量計(jì)算得到的第i個(gè)傳感器的三分量數(shù)據(jù)，Bti為由三分量數(shù)據(jù)計(jì)算得到的總磁場(chǎng)，Bc為參考總磁場(chǎng).

本文針對(duì)磁梯度張量系統(tǒng)校正目的，設(shè)計(jì)如下目標(biāo)函數(shù)：

（17）

由式（17）可知，目標(biāo)函數(shù)包括兩項(xiàng)：第一項(xiàng)為由4個(gè)磁通門傳感器3分量數(shù)據(jù)構(gòu)成的總磁場(chǎng)與質(zhì)子磁力儀測(cè)得的參考總磁場(chǎng)之間的偏差，它主要用于校正單傳感器誤差（零偏、三軸靈敏度不一致和三軸非正交誤差）.第二項(xiàng)為傳感器三分量數(shù)據(jù)的偏差、張量分量和張量不變量，主要用于校正磁梯度張量系統(tǒng)的三軸非對(duì)準(zhǔn)誤差.

3 仿真試驗(yàn)

3.1 仿真試驗(yàn)1

為驗(yàn)證本文方法的有效性，設(shè)計(jì)了仿真模擬試驗(yàn).如表1所示，磁梯度張量系統(tǒng)零偏誤差在-200～200 nT之間隨機(jī)產(chǎn)生，非正交誤差和非對(duì)準(zhǔn)誤差在-2°～2°之間隨機(jī)產(chǎn)生，靈敏度參數(shù)在0.9～1.1之間隨機(jī)產(chǎn)生.假定磁梯度張量系統(tǒng)的基線距離為0.4 m，試驗(yàn)地點(diǎn)處總磁場(chǎng)為52000 nT，磁偏角和磁傾角分別為-7°和54°，使用磁梯度張量系統(tǒng)在三維空間內(nèi)隨機(jī)旋轉(zhuǎn)采集200組數(shù)據(jù).

在實(shí)際應(yīng)用中，磁梯度張量系統(tǒng)采集數(shù)據(jù)可能會(huì)存在一定的測(cè)量誤差，且不可避免地受到非均勻背景場(chǎng)的影響.因此，在仿真試驗(yàn)中加入隨機(jī)噪聲，在4個(gè)磁通門傳感器的3個(gè)測(cè)量軸中分別加入了最大幅值為0 nT、2 nT、5 nT和10 nT的隨機(jī)噪聲.加入不同大小噪聲時(shí)，反演得到的磁梯度張量系統(tǒng)誤差參數(shù)如表1所示，校正后的總磁場(chǎng)和張量分量如圖4所示.當(dāng)仿真模擬數(shù)據(jù)不含隨機(jī)噪聲時(shí)，反演得到的誤差參數(shù)和預(yù)設(shè)參數(shù)完全一致，校正后的總磁場(chǎng)和張量分量分別為52000 nT和0 nT/m.當(dāng)仿真模擬數(shù)據(jù)加入2 nT、5 nT和10 nT的隨機(jī)噪聲時(shí)，反演得到的零偏誤差參數(shù)與預(yù)設(shè)參數(shù)的最大偏差分別為0.367 nT、0.909 nT和1.825 nT，非正交誤差與預(yù)設(shè)參數(shù)的最大偏差為0.001°、0.003°和0.006°，對(duì)靈敏度不一致誤差和傳感器非對(duì)準(zhǔn)誤差6個(gè)誤差參數(shù)的反演結(jié)果沒有影響.由圖4可知，當(dāng)加入最大幅值為10 nT的隨機(jī)噪聲時(shí)，校正后的張量分量仍小于0.2 nT/m.

表1 預(yù)設(shè)誤差參數(shù)與反演得到的誤差參數(shù)Table 1 Preset and estimated error parameters

圖4 加入不同大小隨機(jī)噪聲得到的總磁場(chǎng)和張量分量校正結(jié)果Fig.4 Calibration results of total magnetic intensity and tensor component

圖5 不同種群大小得到的校正數(shù)據(jù)RMS隨迭代次數(shù)的變化情況Fig.5 RMS error of calibration results with different population size

差分進(jìn)化算法的反演精度與設(shè)定的種群數(shù)量NP、變異因子F和交叉概率因子CR有關(guān)（Balkaya et al.,2017）.本文分析了這些參數(shù)對(duì)反演結(jié)果的影響.將種群數(shù)量分別設(shè)定為10、15、20、50、100和200，分別反演200組不含噪聲的磁場(chǎng)數(shù)據(jù).令最大迭代次數(shù)為10000，利用每次迭代中最佳個(gè)體校正磁場(chǎng)數(shù)據(jù)，并計(jì)算校正數(shù)據(jù)的RMS.RMS計(jì)算公式如式（18）所示.圖5為不同種群數(shù)量時(shí)校正數(shù)據(jù)RMS隨迭代次數(shù)的變化情況，可以看出隨著種群數(shù)量增大，差分進(jìn)化算法陷入局部最優(yōu)解的可能性降低.當(dāng)種群數(shù)量為50時(shí)即可實(shí)現(xiàn)精度較高的反演.

RMS=

（18）

式中，N為磁場(chǎng)數(shù)據(jù)量.

假定變異因子和交叉概率因子的取值范圍分別為[0.2,0.9]和[0.2,1]，我們分別采用傳統(tǒng)差分進(jìn)化算法和改進(jìn)差分進(jìn)化算法估計(jì)張量系統(tǒng)誤差參數(shù)，并計(jì)算校正后總磁場(chǎng)與張量分量RMS的和，結(jié)果如圖6所示.假定校正后總磁場(chǎng)與張量分量RMS之和小于0.01時(shí)校正精度是可以接受的.由圖6可知，與傳統(tǒng)差分進(jìn)化算法相比，改進(jìn)差分進(jìn)化算法對(duì)變異因子和交叉概率因子依賴性較小.本文算例中F0取值為0.3，CR0取值為0.8.

參與反演的數(shù)據(jù)量可能也會(huì)影響誤差參數(shù)的反演精度，為分析數(shù)據(jù)量對(duì)反演結(jié)果的影響，令磁梯度張量系統(tǒng)在三維空間中隨機(jī)旋轉(zhuǎn)采集10組、20組、50組、100組、200組和500組磁場(chǎng)數(shù)據(jù)，并將其輸入改進(jìn)差分進(jìn)化算法反演誤差參數(shù).令種群數(shù)量為200，最大迭代次數(shù)為10000，利用每次迭代中最佳個(gè)體校正磁場(chǎng)數(shù)據(jù)，并計(jì)算校正數(shù)據(jù)的RMS.由圖7可知，當(dāng)數(shù)據(jù)量為10組時(shí)，無法得到良好的校正效果，而當(dāng)數(shù)據(jù)量為50組時(shí)即可實(shí)現(xiàn)誤差參數(shù)的精確反演.

將本文方法所得結(jié)果與遺傳算法、最小二乘擬合方法的校正結(jié)果進(jìn)行了對(duì)比.在反演中使用了200組磁場(chǎng)數(shù)據(jù)，最大迭代次數(shù)為10000，種群數(shù)量分別設(shè)置為50和200，利用每次迭代中最佳個(gè)體校正磁場(chǎng)數(shù)據(jù)，并計(jì)算校正數(shù)據(jù)的RMS.改進(jìn)差分進(jìn)化算法和遺傳算法反演過程中RMS隨迭代次數(shù)變化情況如圖8所示.由圖8可知，改進(jìn)差分進(jìn)化算法的搜索效率與反演精度明顯高于遺傳算法.當(dāng)種群數(shù)量為50時(shí)，改進(jìn)差分進(jìn)化算法迭代次數(shù)約為2000時(shí)反演得到的RMS收斂于10-4量級(jí)，而遺傳算法迭代10000次得到的RMS仍為17.68.

校正前后磁場(chǎng)數(shù)據(jù)的均方根誤差如表2所示.由表2可知，當(dāng)不存在隨機(jī)噪聲時(shí)，本文方法和最小二乘方法的校正結(jié)果相同，均能準(zhǔn)確反演出磁梯度張量系統(tǒng)的誤差參數(shù).添加隨機(jī)噪聲后，本文方法反演結(jié)果的RMS明顯小于最小二乘方法，本文方法具有較強(qiáng)的抗噪能力.與遺傳算法相比，本文方法反演結(jié)果的RMS小于遺傳算法反演結(jié)果的RMS，這意味著本文方法的反演精度較高.

表2 磁場(chǎng)數(shù)據(jù)校正前后的均方根誤差Table 2 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration

圖6 兩種差分進(jìn)化算法校正后磁場(chǎng)數(shù)據(jù)的均方根誤差（a） 傳統(tǒng)差分進(jìn)化算法；（b） 改進(jìn)差分進(jìn)化算法.Fig.6 RMS error of calibration results with differential evolution algorithm（a） Traditional differential evolution algorithm；（b） Improved differential evolution algorithm.

圖7 不同數(shù)據(jù)量得到的校正數(shù)據(jù)RMS隨迭代次數(shù)的變化情況Fig.7 RMS error of calibration results with different magnetic data

圖8 改進(jìn)差分進(jìn)化算法和遺傳算法反演結(jié)果對(duì)比Fig.8 Comparison of RMS error of calibration results with improved differential evolution algorithm and genetic algorithm

3.2 仿真試驗(yàn)2

為進(jìn)一步測(cè)試本文方法的實(shí)用性和可靠性，依據(jù)某國(guó)產(chǎn)磁通門傳感器測(cè)定的誤差大致取值范圍，設(shè)計(jì)了仿真模擬試驗(yàn).如表3所示，磁梯度張量系統(tǒng)的零偏誤差在-100～100 nT之間隨機(jī)產(chǎn)生，非正交誤差在-0.1～0.1度之間隨機(jī)產(chǎn)生，靈敏度參數(shù)在0.9999～1.0001之間隨機(jī)產(chǎn)生，非對(duì)準(zhǔn)誤差在-1～1度之間隨機(jī)產(chǎn)生.同樣，假定磁梯度張量系統(tǒng)的基線距離為0.4 m，試驗(yàn)地點(diǎn)處總磁場(chǎng)為52000 nT，磁偏角和磁傾角分別為-7°和54°，使用磁梯度張量系統(tǒng)在三維空間內(nèi)隨機(jī)選擇采集200組數(shù)據(jù)，在磁場(chǎng)數(shù)據(jù)中加入了最大幅值為0 nT、2 nT、5 nT和10 nT的隨機(jī)噪聲.

由表4可知，當(dāng)仿真模擬數(shù)據(jù)不含隨機(jī)噪聲時(shí)，磁梯度張量系統(tǒng)的誤差參數(shù)可以得到精確的反演.當(dāng)仿真模擬數(shù)據(jù)加入2 nT、5 nT和10 nT的隨機(jī)噪聲時(shí)，反演得到的零偏誤差參數(shù)與預(yù)設(shè)參數(shù)的最大偏差分別為0.339 nT、0.844 nT和1.715 nT，非正交誤差與預(yù)設(shè)參數(shù)的最大偏差為0.001°、0.003°和0.006°，靈敏度誤差與預(yù)設(shè)參數(shù)的最大偏差為2.2×10-5、5.6×10-5和1.13×10-4.由表5可知，經(jīng)過校正后，磁梯度張量系統(tǒng)的測(cè)量誤差得到了較大程度地壓制，當(dāng)仿真模擬數(shù)據(jù)加入10 nT的隨機(jī)噪聲時(shí)，校正后總磁場(chǎng)、張量分量Bxx,Bxy,Bxz,Byy,Byz的均方根誤差分別為7.327、0.379、0.211、0.151、0.044和0.075.

表3 預(yù)設(shè)誤差參數(shù)與反演得到的誤差參數(shù)Table 3 Preset and estimated error parameters

表4 磁場(chǎng)數(shù)據(jù)校正前后的均方根誤差Table 4 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration

4 實(shí)測(cè)試驗(yàn)

為驗(yàn)證本文方法的實(shí)用性，在地磁干擾比較小的夜間開展了試驗(yàn)工作.試驗(yàn)中將十字形磁梯度張量系統(tǒng)安裝于三軸無磁旋轉(zhuǎn)臺(tái)上（圖9a）隨機(jī)旋轉(zhuǎn)采集數(shù)據(jù)，采用質(zhì)子磁力儀（圖9c）測(cè)量試驗(yàn)場(chǎng)地的地磁場(chǎng).磁梯度張量系統(tǒng)基線距離為0.16 m，質(zhì)子磁力儀分辨率為0.1 nT、測(cè)量精度為1 nT，三軸無磁旋轉(zhuǎn)臺(tái)旋轉(zhuǎn)角度可讀分辨率為1°.為降低試驗(yàn)中地磁場(chǎng)的天然變化對(duì)試驗(yàn)結(jié)果的影響，本次試驗(yàn)中只采集了60組數(shù)據(jù)，以縮短試驗(yàn)所需時(shí)間.

圖9 試驗(yàn)所用儀器（a） 張量系統(tǒng)與無磁三軸旋轉(zhuǎn)臺(tái)；（b） 張量系統(tǒng)記錄儀；（c） 質(zhì)子磁力儀.Fig.9 Test equipment（a） Magnetic gradient tensor system and three-axis turntable；（b） Data acquisition system；（c） Proton magnetometer.

利用本文提出的校正方法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行了校正，校正前和校正后的總磁場(chǎng)和張量分量如圖10和圖11所示.磁梯度張量系統(tǒng)誤差未校正前測(cè)量值波動(dòng)較大，地磁場(chǎng)總磁場(chǎng)的波動(dòng)值可達(dá)數(shù)百 nT，磁梯度張量分量的最大偏差可達(dá)1000 nT/m，這說明磁梯度張量系統(tǒng)的誤差會(huì)對(duì)測(cè)量結(jié)果產(chǎn)生較為嚴(yán)重的影響，必須對(duì)其進(jìn)行校正.從圖10和圖11中可以看出，經(jīng)過本文方法校正后，總磁場(chǎng)均收斂至試驗(yàn)地點(diǎn)處地磁場(chǎng)（51600 nT）附近，張量分量均收斂至0 nT/m附近.

圖10 總磁場(chǎng)校正前后對(duì)比Fig.10 Comparison of total magnetic intensity before and after calibration

圖11 張量分量校正前后對(duì)比Fig.11 Comparison of tensor components before and after calibration

校正前后總磁場(chǎng)與張量分量的均方根誤差如表5所示.由表5可以看出，經(jīng)過校正后，磁場(chǎng)數(shù)據(jù)的均方根誤差明顯減小，且本文方法的校正結(jié)果明顯優(yōu)于基于遺傳算法和最小二乘方法的誤差校正方法.以總磁場(chǎng)為例，校正前總磁場(chǎng)均方根誤差為179.566 nT，基于本文方法、遺傳算法和最小二乘方法的總磁場(chǎng)均方根誤差分別為3.398 nT、18.197 nT和14.210 nT.

表5 磁場(chǎng)數(shù)據(jù)校正前后的均方根誤差Table 5 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration

5 結(jié)論

磁梯度張量系統(tǒng)的測(cè)量精度受到磁場(chǎng)傳感器零偏、靈敏度不一致、非正交和傳感器非對(duì)準(zhǔn)誤差的影響.本文提出了一種基于改進(jìn)差分進(jìn)化算法的磁梯度張量系統(tǒng)誤差校正方法，只需一步反演即可實(shí)現(xiàn)磁梯度張量系統(tǒng)48個(gè)誤差參數(shù)的估計(jì).用兩個(gè)仿真模擬算例驗(yàn)證了本文方法的有效性，仿真結(jié)果表明本文方法與遺傳算法相比具有更高的求解精度和搜索效率，本文方法較最小二乘擬合方法具有更高的抗噪能力.改進(jìn)的差分進(jìn)化算法能有效避免傳統(tǒng)差分進(jìn)化算法搜索效率低、易陷入局部最優(yōu)解、對(duì)變異因子和交叉概率因子取值依賴大等問題.實(shí)測(cè)數(shù)據(jù)總磁場(chǎng)、張量分量Bxx,Bxy,Bxz,Byy,Byz校正后RMS改善比分別為52.845、54.530、25.415、40.057、38.330和25.586，本文方法能夠有效降低磁梯度張量系統(tǒng)誤差引起的測(cè)量誤差，提高測(cè)量精度.本論文的不足在于實(shí)際試驗(yàn)場(chǎng)地磁場(chǎng)環(huán)境不夠理想，導(dǎo)致實(shí)測(cè)數(shù)據(jù)總磁場(chǎng)校正后的均方根誤差仍為3.398 nT.在下一步的工作中需尋找環(huán)境干擾更小的場(chǎng)地進(jìn)行試驗(yàn)，以更好地估計(jì)磁梯度張量系統(tǒng)的誤差參數(shù).

致謝三位審稿人提出了建設(shè)性意見，對(duì)本文的質(zhì)量提升幫助很大，在此表示誠(chéng)摯的感謝.