李青竹,石志勇,李志寧,范紅波
(陸軍工程大學 車輛與電氣工程系,河北 石家莊 050003)
磁異常探測技術可應用于地下小尺度磁目標的定位與識別[1-2]。與磁場矢量和總磁場強度(Total Magnetic Intensity,TMI)相比,磁梯度張量(Magnetic Gradient Tensor,MGT)可以提供更豐富的目標體方位信息,抗干擾能力強,能夠更好地適應復雜的測量環境[3]。MGT探測具有廣泛的應用前景,如航空磁測、礦產勘探、未爆彈藥搜索和排雷等[4-5]。
直接測量磁場梯度本質上很困難,但可利用矢量傳感器在基線兩端的讀數差異近似估計磁標勢的二次偏微分[6]。目前,常以磁傳感器陣列形式構建磁梯度張量系統(Magnetic Gradient Tensor System,MGTS),并實現張量差分測量。MGTS主要分為兩類:(1)基于超導效應[7],這類系統由具有極高靈敏度和較小量程的超導量子干涉裝置(Superconducting Quantum Interference Devices,SQUID)組成,其制造成本高,對測量環境要求嚴格,適用于生物磁檢測、金屬無損檢測、航空磁測等靈敏度要求高但磁異常較弱的工況;(2)基于磁通門法[8-9],這類系統由多個磁通門傳感器組成,利于批量生產和制造,成本較低,安裝要求較簡單。目前,最先進的磁通門探頭的靈敏度噪聲可以達到6 p Trms/ Hz@1 Hz[10],量程是SQUID的數千倍[11-12]。
近年來,國內外研究團隊搭建了各類磁通門法MGTS,包括直角四面體、正四面體、正方形、十字形、三角形等結構[13-16],但針對特定結構MGTS探測極限的研究卻鮮有報道。實測經驗表明,差分方法測量磁梯度張量時系統的探測極限常受到結構誤差、噪聲、磁源強弱和探測方位等因素的影響[17]。為了定性且定量研究MGTS的理論探測極限,本文利用磁偶極子正演方程、張量矩陣特征方程和張量不變量聯合推導出在系統基線距離、傳感器測量精度、系統觀測方位和磁源強度等相關參數約束下MGTS空間理論探測范圍的計算公式。針對單一磁源,張量衍生不變關系定位方法[18]利用MGTS實現單點精確定位并準確估計出磁源的磁矩強度,使實測MGTS的理論探測極限成為可能。
MGT為3個正交方向的磁場矢量空間變化率[3],共9個元素,可表示為:

其中:B為磁場強度矢量,G表示MGT矩陣,φm是磁標勢,Bm(m=x,y,z)表示B的正交分量,Bij(i,j=x,y,z)表示MGT分量。在沒有電流的靜態磁環境中,麥克斯韋方程約束下磁場的散度和旋度為0,有?·B=0,?×B=0,因此G是對稱且無跡的。故G的9個元素可以用5個獨立的分量來表示,即Bxx,Bxy,Bxz,Byy和Byz。
MGT是磁標量勢φm的二階偏微分,也是磁場矢量B的偏微分,難以直接測量。事實上,常利用跨測量基線的分量讀數差異來代替磁矢量偏微分,以近似MGT分量,如:

其中:ΔBi是兩個相鄰磁傳感器測量的i分量讀數差異,Δdj是兩個磁傳感器在j方向上的距離,定義為基線距離。
MGT不變量是當觀測點固定不隨測量系統不同方向而改變的量或對應關系。一些基本的MGT不變量包括矩陣跡、特征值、特征方程系數和Frobenius范數等[19]。MGT矩陣是實對稱矩陣,可對角化。令λ1,λ2和λ3表示G的特征值,滿足特征方程λ3-I0λ2+I1λ-I2=0,其中:

其中:I0是G的跡,I1是關于G對稱且無跡的量,I2是矩陣G的行列式,CT是張量收縮。給定一個磁偶極子,磁矩為m=(mx,my,mz)[20],偶極子到觀測點的位置矢量為r=(x,y,z)。設M=||m||,r=||r||,r為觀測距離,則磁場矢量和5個獨立MGT分量的正演方程為[19-20]:

其中:μ為介質的磁導率,空氣中μ≈μ0=4π×10-7N·A-2,μ0為真空磁導率。
空間衍生不變關系如下[18]:

其中:θ是r和m之間的角度,0≤θ≤180°,u是歸一 化 源 強 度(Normalized Source Strength,NSS)[19]。
由式(3)、式(4)、式(6)分別得到MGT矩陣的特征值λ1,λ2和λ3,即有:

則NSS可以表示為:

其 中:|λ1|≥|λ3|,|λ2|≥|λ3|和λ2≥λ3≥λ1。式(8)可為MGTS的理論探測極限估計提供條件約束。
相關研究表明,相較于正四面體、直角四面體、三角形和正方形等結構,平面十字形MGTS不僅結構簡單、傳感器安裝方便、對準精度高,其結構誤差較小且同基線距離下的最近觀測范圍最大[17]。圖1為平面十字形MGTS的結構設計,由4個磁通門傳感器和十字架構成。設d為基線距離,表示傳感器a在b軸方向的讀數。結合式(2),平面十字形MGTS的張量測量矩陣Gm為:


圖1 平面十字形MGTS的結構Fig.1 Structure of planar-cross MGTS
式中:Byx和Bxy的測量方式存在差異(事實上二者實測值幾乎相同),故Gm共有6個獨立分量。
令n為一單位空間向量,?(B)/?(n)為B對n的方向導數。用MGT矩陣G的特征向量v1,v2和v3表示n,有n=a v1+b v2+c v3且a2+b2+c2=1,得到:

定義q為傳感器的測量準確度(單位:±n T),Q為MGTS的測量準確度(單位:±nT/m),結合式(9)有q=Q d/2。注意,Q和q不僅反映了傳感器本底噪聲引起的測量誤差,還包括傳感器系統誤差、磁干擾誤差和未對準誤差引起的不確定性[8]。將考慮全部實測誤差因素時的最大測量偏差作為準確度Q和q的取值范圍。
在觀測點r0=(x0,y0,z0)處捕捉目標需要q滿足可靠的差分計算。因此,q要小于磁場強度在基線距離之間的實際變化值。令r為MGTS的觀測距離,結合式(8)和式(10),傳感器測量準確度q需滿足:

使得在基線距離d、該處磁異常強度B條件下的差分計算可靠。
差分測量可靠的必要條件如下:

式(12)即為差分磁梯度張量測量范圍的計算公式。由此可知,某一MGTS系統的理論探測極限不僅與傳感器的測量準確度q和系統基線距離d有關,還與磁源目標的磁矩大小M和觀測角度θ有關。這表明MGTS對同一目標在不同方向上具有不同的探測極限。
由于式(12)中與磁源相關的磁矩M和觀測角度θ為單一磁源目標參數,當測區內存在多個磁異常目標時,M應為空間某處的磁偶極子等效源模型磁矩。該等效源是由各個磁異常目標在觀測點處的磁場疊加后再反演得到的磁偶極子,θ應為該等效源磁偶極子的磁矩與觀測點間的夾角。
在MATLAB仿真環境中,建立空間笛卡爾坐標系,令x軸朝東,y軸朝北,z軸朝上。將磁偶極子放置在r1=(10 m,0,-20 m)處,磁矩M=5 000 A·m2,磁偏角為20°(東偏),磁傾角為-60°(向上)。設置MGTS的初始參數為基線距離d=0.5 m,傳感器的測量準確度q=±0.1 n T。圖2(a)中繪制了MGTS在三維空間中的理論探測極限。圖2(b)顯示了在z=-20 m平面內不同d和q時MGTS的理論探測極限。圖3(a)顯示了空間中不同θ和磁矩強度M時MGTS的理論探測極限距離。圖3(b)顯示MGTS最大探測極限距離同時隨d和q的變化而變化。
顯然,在平行于磁矩矢量的方向(即θ=180°或0°),MGTS的探測距離最大。q越小,d越長,M越大,理論探測極限距離越遠。
然而實際情況表明,當探測距離不變而d值增長到一定程度時,MGT測量值會失真。由于傳感器差分計算得到的張量值與真實值之間存在理論誤差,以分量Bxx為例說明式(3)計算得到的MGT分量值沿x軸梯度測量方向的逼近過程。間隔d的兩傳感器x軸測量值和分別自展開為泰勒級數,其中為觀測點處磁場矢量的x軸真實分量。令Bxx測量值為,真實值為,則可由泰勒級數表示:


圖2 MGTS在變參數下的空間探測極限范圍Fig.2 Detection limits of MGTS with variable parameters

圖3 MGTS在不同磁源和變參數下的理論探測極限Fig.3 Theoretical detection limits of MGTS with different sources and variable parameters
式(13)表明,在差分過程中,由于忽略了泰勒級數的三階和高階奇數項,測量值與基線距離d存在一定的正相關偏差。當d值增大且觀測距離很近時,磁場高階偏微分也變得不可忽略。此時,式(13)中第二項和后續項與首項相比不能忽略,差分測量結果失真愈發明顯。
圖4顯示了系統由窗口(10 m,0,-20 m)滑動測量至(100 m,0,-20 m),即測點由磁偶極子中心逐漸遠離的過程中Bxx和Bxy的理論測量誤差。因此,隨著觀測距離變短,由差分計算的MGT測量值可靠性會變低,且d越長這種趨勢越明顯。從這些結果來看,MGTS的理論探測極限與基線距離、測量準確度、磁矩、觀測點與磁矩矢量的夾角等有關;基線距離越長,傳感器測量準確度越高,MGTS的理論可探測距離越遠;基線距離越長,觀測距離越近,張量差分測量理論誤差越大;探測距離在平行于磁矩矢量的方向上達到最大,而在垂直于磁矩矢量的方向上急劇減小。

圖4 不同基線距離和觀測距離時MGT差分測量的理論誤差Fig.4 Theoretical error of MGT differential measurement with different baseline distances and observation distances
在工程實際中,基線距離的設置存在一個較優解,需綜合考慮精度需求、距離需求及儀器尺寸等。
本文構建了一個實際的平面十字形MGTS,基線距離為0.5 m,其中包含4個Barrington公司生產的Mag-03磁通門傳感器,以及一個由非磁性塑性樹脂材料制成的十字架,見圖5。
將圖6中4種典型的磁體標記為m1,m2,m3和m4,其中m1和m2為圓柱體,m3和m4為長方體。這些磁體的等效偶極矩未知。
實驗過程如圖7所示,以實現典型磁體的MGTS探測極限估計和驗證。磁異常測量經驗法表明,在超過物體長度約2.5倍的距離處,偶極矩會占主導作用[21],磁目標可近似于磁偶極子。將4個磁鐵視為遠離目標位置的磁偶極子,并試圖找到它們相對于幾何形狀的磁矩矢量。

圖5 平面十字形MGTS實驗裝置Fig.5 Experimental setup of planar-cross MGTS

圖6 預備的4塊典型形狀磁鐵及其尺寸信息Fig.6 Prepared four magnets with size information

圖7 MGTS探測極限估計和驗證實驗Fig.7 Detection limits estimation and verification experiments
圖7中,環境測量用以估算傳感器和MGTS的實際測量準確度,Q為MGTS的測量準確度(單位:±n T/m),有q=Q d/2。采用張量衍生不變關系定位方法[18]來估計磁鐵位置,若位置處于磁鐵的物理尺寸內,則表明磁矩估計是有效的。然后,利用放置磁鐵前后的MGT測量讀數和準確度Q判斷MGTS是否已到達其探測極限。在測量中,這里采用零相位低通濾波來消除信號干擾。
Mag-03傳感器的出廠本底噪聲幅度在±0.01~0.02 n T內[10],然而傳感器的真實讀數幾乎不可能達到此量級。實際上,實測讀數不可避免地受到測量設備的電磁干擾,地磁場的本底噪聲以及環境中其他未知渦流磁信號的干擾。通過信號處理的方法能夠過濾掉一些有規律的高頻信號,盡可能保留真實的磁場信息。
實驗地點為中國石家莊某空曠野外,劃定10 m×10 m且磁場較穩定的測區,標量質子磁強計測得該區域的磁總場強度均值為53 162 n T(±20 nT)。首先進行了環境測量,以確定濾波和降噪的截止頻率。采樣頻率為500 Hz。圖8(a)、8(b)顯示了經過約54 s靜態采樣后傳感器1的x軸分量和頻譜,并突出了前后各一秒內的讀數變化細節。靜磁場信號是一種超低頻信號,而設備電流采集卡的噪聲頻段主要集中在50,235 Hz等頻段。利用零相位低通濾波器來降低地磁信號噪聲,并將截止頻率保守地設置為5 Hz。由于地磁日變的影響以及濾波器中遺漏了部分低頻噪聲分量,實際磁場強度在不斷變化,傳感器測量準確度q難以直接獲得。然而,在MGT數據的差分計算過程中,可以等效地減去該區域地磁場的協同變化。故若已知MGTS的測量準確度Q,可以用差分原理估計q,即q=Q d/2。
圖8(c)顯示了開始采樣后一段時間4個傳感器的磁場強度分量讀數,各軸讀數明顯地隨地磁場波動而變化。由于采樣時間一分鐘內溫度相對穩定,并且傳感器的出廠偏移溫度系數(n T/°C)小于±0.1 n T[10],因此圖8(c)中的磁場讀數漂移幾乎全部來自地磁在短時間內的隨機日靜變化,此時直接估計q是不可靠的。然而,在靜態采樣期間測量的MGT應該是恒定的,因為它代表了磁場強度的空間變化率,不受地磁日變化的影響。MGT讀數隨采樣時間的絕對誤差如圖8(d)所示。Mag-03傳感器的出廠測量準確度已達到±0.02 n T,確保了此時MGT讀數的可靠性。表1列出了MGTS讀數的精確度(均方根誤差(Root Mean Square Error,RMSE)、平均誤差(Mean Error,ME)和準確度Q。

圖8 磁場低通濾波截止頻率的選擇與MGTS張量測量準確度的估計結果Fig.8 Selection of cutoff frequency of low-pass filter and estimation of MGTS tensor measurement accuracy
在約54 s的靜態采樣階段內,測得MGT分量讀數的最大絕對誤差收斂到±8.285 n T/m,作為MGTS的估計張量測量準確度Q。重復采樣后得到的MGT值也在此范圍內波動,表明局部磁場相對穩定且Q值估計有效。由于q=Q d/2,給定d=0.5 m,故q真值應收斂于±2.121 n T附近,作為現有儀器和降噪條件下估計的傳感器測量準確度。

表1 MGTS讀數的估計準確度(最大絕對誤差Q)和精確度(RMSE和ME)Tab.1 Estimated accuracy(maximum absolute error Q)and precision(RMSE and ME)of MGTS readings(n T/m)
一旦獲得MGT分量和位置矢量,可直接利用式(5)磁偶極子正演方程來反演磁矩矢量m,即:

其中:H僅與位置矢量有關,I是具有5個獨立MGT分量的列向量。請注意,式(13)中的H+是Moore-Pennrose逆。然而,現有情況不能直接判斷估計的m是否有效。環境中可能存在其他磁異常,使得由反演公式得到的疊加磁場等效磁矩代替了目標的真實磁矩。此外,前文結論表明,距離過近同樣會使張量差分測量值變得不可靠。
為了確保磁矩估計是有效的,利用張量衍生不變關系法定位磁體[18]。該方法從以下方程中提供4個可能的坐標解,其中一個是正確的。

式中v1和v2分別對應于特征值λ1和λ2的特征向量,CT是張量縮并,是僅由磁偶極子產生的磁總場強度TMI。因此,在計算ITM時應提前測量地磁場進行補償。
對坐標已知的磁體進行定位,一旦估計坐標與目標真實位置間的測量偏差控制在磁體尺寸范圍內,則反演得到的磁偶極矩是有效的。MGTS置于地面高50 cm處,以觀測點為原點,將4塊磁鐵分別放置于位置1(50 cm,50 cm,-50 cm)、位置2(100 cm,100 cm,-50 cm),如圖9所示。

圖9 磁鐵磁矩估計實驗Fig.9 Magnets magnetic moment estimation experiments
對4塊磁鐵共進行了4×2次定位實驗。采樣頻率為500 Hz,單次采樣時間為10 s左右。定位結果為每單次采樣時間內計算的平均值。設置相同的濾波條件,截止頻率均為5 Hz。定位結果列于表2。顯然,位置1處的定位結果均有偏差,而位置2處的定位結果均控制在4個磁體的幾何空間范圍內。位置1處由于探測距離過近而差分測量失真,而位置2的結果都是有效的。4個磁鐵估計的磁矩矢量如圖10所示。

表2 磁鐵定位實驗中估計的磁鐵位置和磁偶極矩Tab.2 Estimated magnet position and magnetic dipole moment in magnet positioning experiments

圖10 估計的4塊磁鐵的磁偶極矩大小和方向Fig.10 Estimated magnitude and direction of magnetic dipole moment of four magnets
一旦確定了磁矩的方向和大小、傳感器讀數準確度以及MGTS的基線距離,就可以得到該MGTS對4個磁體的空間理論探測極限范圍,如圖11(a)和11(b)所示。此外,圖11(c)顯示了4個磁鐵的探測極限隨角度θ的變化。MGTS相關參數和最大探測極限距離(rmax)的估計結果列于表3。
為了驗證估計結果,將4塊磁鐵沿MGTS的x軸滑動并連續采樣。采樣率為500 Hz,滑動窗口為(1 m,0,0)到(10 m,0,0)。規定每0.01 m采集不少于20個點時,實測MGT的準確度可信,故滑動速度應小于0.25 m/s。設置滑動速度為0.2 m/s,速度偏差為±0.05 m/s,即可滿足要求。滑動過程中始終保持角度θ為0°,即磁矩m的方向始終指向觀察點。


圖11 MGTS對4塊磁鐵的探測極限估計Fig.11 Detection limits estimation of MGTS for four magnets
通過觀察是否可以區分不同MGT分量讀數來判斷MGTS是否已達到對4種磁鐵的真實探測極限邊界,并估計出實測rmax。為了消除環境中其他未知磁異常對結果的干擾,利用背景磁場的空側數據對滑動采樣階段測量的磁場信號進行補償。一旦磁體滑動時測得的MGT分量讀數之差不再超過MGTS的測量準確度Q(±8.285 n T/m),判定磁體已脫離了系統的有效探測范圍。測量結果如圖12所示,圖中標出了每個磁體的實測探測極限rmax,并記錄在表3中,其中測量值和估計值最大偏差為0.4 m。
上述結果表明,MGTS實測探測極限與估計值吻合較好,探測極限估計準確度為±0.4 m,從而驗證了所提差分磁梯度張量測量極限估計方法的有效性。

表3 MGTS探測極限估計及其必要參數Tab.3 Detection limits estimation and its necessary parameters of MGTS

圖12 MGTS在平行于4塊磁鐵磁矩方向上測得的有效探測極限距離(rmax)Fig.12 Measured MGTS detection limit distance(rmax)of four magnets in direction parallel to magnetic moment
本文提出了差分MGT測量極限估計方法。根據差分MGT測量范圍公式,MGTS的理論探測極限與基線距離、傳感器測量準確度、目標磁偶極矩、觀測點與磁矩矢量間的夾角有關。基線距離越長,傳感器測量準確度越高,MGTS的有效探測距離越遠。探測距離在平行于磁矩矢量的方向上達到最大,而在垂直于磁矩矢量的方向上急劇減小。在實際測量中,成功估計并驗證了搭建的平面十字形MGTS的理論探測極限,該系統針對4塊典型磁鐵的探測極限估計準確度為±0.4 m。然而,本文尚未深入考慮地磁日變化、渦流磁場干擾和磁傳感器非線性誤差等因素對實測中系統探測極限的影響。