基于Matlab的水下自航目標體三分量磁陣列探測仿真系統

沈小兵， 李永濤

(中國地質大學(武漢) 地球物理與空間信息學院， 湖北 武漢 430074)

?

·專題研討——模擬仿真實驗(27)·

基于Matlab的水下自航目標體三分量磁陣列探測仿真系統

沈小兵， 李永濤

(中國地質大學(武漢) 地球物理與空間信息學院， 湖北 武漢 430074)

基于Matlab中GUI平臺程序開發環境，開發模擬三分量磁通門磁陣列探測仿真系統,進行實時探測水下自航磁性目標體的運動狀況，取得了如下仿真效果：通過控件實現了自航磁性目標體的三分量磁異常可視化數值模擬計算功能；運動中的磁性目標體Za、Hax、Hay三分量及總磁異常的特征分析；根據自航磁性體運動特征實現自航磁性體運動狀態及磁異常曲線的可視化顯示；根據用戶輸入的不同參數，實現磁性目標體在不同深度、不同半徑、不同長度、不同運動速度以及不同環境磁場的條件作用下的三分量磁異常變化的仿真模擬技術。經過該仿真模擬實驗，加深了對三分量磁探測水下磁性目標體方法技術的了解，為學生從事地球物理勘探工作以及軍事地球物理均具有幫助。

Matlab； 三分量； 磁陣列； 仿真系統； 自航目標體； 磁異常

0 引 言

磁法勘探研究中的磁異常是指磁性體產生的磁場疊加在地球磁場之上而引起的地磁場畸變;它是一個空間矢量場,可以通過它在3個坐標軸的分量或是正常地磁場方向的投影量來確定,也可以進一步用各分量的垂向、水平梯度來表示。所以,磁異常是一個多分量磁場。用高精度磁力儀觀測獲得磁異常多分量信息是磁法勘探的一個重要環節,可靠的磁異常多分量數據是磁法勘探的重要基礎[1-2]。在以往地面磁法勘探測量磁場的單一分量中,不管是水平分量、垂直分量還是總場,在解釋的時候通常是繪制平面等值線圖、剖面圖,通過化極、延拓等輔助方法來推斷場源中心、埋深以及范圍等。如果可以直接得到測量磁場的各個分量,繪制精確的含有矢量信息的圖件來反映異常體的各種特征,在后期的解釋過程中,可以增加反演結果的可靠程度,并且也可以根據多分量帶來的空間矢量場的方向信息,增加對磁異常的空間認識[4,8-12]。

基于三分量磁信號的目標探測技術是近年來隨著磁探測技術的不斷發展和磁探測傳感器的測量的測量精度的不斷提高而新興起的一種目標磁探測技術。它以含有鐵磁性物質的物體擾動地磁場分布，從而產生地磁場分布異常的物理現象為基礎，通過測量磁異常信號的分布，提取磁異常信號的特征量，并通過一定的數據處理最終得到目標相關信息。該項技術在軍事領域擁有極其廣闊的發展前景，比如說，隱形目標搜索，戰場監控，水雷位置探測等[3]。

Matlab是由MathWorks公司開發的一種主要用于數值計算及可視化圖形處理的工程語言[6],它將數值分析、矩陣運算、圖形圖像處理、信號處理和仿真等諸多強大功能集成在較易使用的交互式計算機環境中，為科學研究、工程應用提供了一種功能強、效率高的編程工具，它擁有強大的科學計算與可視化功能、簡單易用、開放式可拓展環境，特別是所附帶的30多種面向

不同領域的工具箱支持，使得它在許多科學領域中成為計算機輔助設計和分析、算法研究和應用開發的基本工具和首選平臺[13-20]。

基于Matlab的水下自航目標體三分量磁陣列探測仿真系統,將有助于地球物理學專業及磁法勘探領域學生得到更好的訓練和提高。

1 自航磁性體的磁異常正演

根據物理學原理，流線型物體運動時受到的阻力最小，所以水下運動磁性體的形狀多為長橢圓狀，在一定程度上，可以以一運動的磁性水平圓柱體來代替運動中的磁性體，兩者所產生的磁異常特征基本相同，因此討論水平圓柱體的磁異常特征不僅有實際的意義，也有一定的合理性[1,5]。

水平圓柱體的磁場可視為許多沿走向緊密排列的磁偶極子的磁場之和。因此,將均勻磁化球體的磁場沿走向積分,就求得斜磁化有限長水平圓柱體的磁場,其公式如下:

(1)

(2)

(3)

(4)

其中：觀測點P與圓柱體右截面圓心距離

觀測點P與圓柱體左截面圓心距離

M1=J·S為沿水平圓柱體走向單位長度上的磁矩

Mx=M1cosIcosA,My=M1cosIsinA,Mz=M1sinI

Za為斜磁化有限長水平圓柱體磁場垂直分量；Hax為斜磁化有限長水平圓柱體磁場的水平X方向分量；Hay為斜磁化有限長水平圓柱體磁場的水平Y方向分量；ΔT為斜磁化有限長水平圓柱體磁場的總磁異常；x、y、z分別為觀測點P的X、Y、Z的坐標值；R為水平圓柱磁性體的中心埋深；2L為沿走向的長度；A為剖面磁方位角(X為剖面方向)，以逆時針轉為正；I為地磁傾角。坐標系與水平圓柱磁性體的空間位置如圖1所示。

式(1)～(4)分別為斜磁化有限長水平圓柱體磁場的垂直分量,水平X分量、水平Y分量及總磁異常公式。由公式可見,它不同于無限長水平圓柱體的磁場。由于My的出現,使磁場有了端面影響,由于My和L的出現,使三度體的磁場公式比二度體復雜得多[7]。

2 自航磁性體可視化編程

2.1 仿真原理

將海底面設想為水平面，并以一定的網度布置三分量磁力測量陣列，當水中存在自航磁性目標體時，每個陣列坐標磁探頭可測量自航磁性體的Za、Hax、Hay三分量及總磁異常值的大小及其相應的矢量參數。通過空間磁異常實時變化，直觀顯示自航磁性體的運動時的磁場特征，達到水下磁性目標體在不同深度、不同半徑、不同長度、不同運動速度以及不同環境磁場條件下的仿真效果，實現良好的人機交互功能(圖2所示)。

2.2 GUI界面設計

(1) 在命令窗口輸入命令GUIDE。啟動對象編輯器，然后對要設計的圖形用戶界面進行總體布局設計，布局的要求為整體布局合理、簡潔、美觀、操作性強。

(2) 添加異常實時圖像顯示功能對象。仿真界面的磁異常圖像顯示包括磁異常的平面等值線圖(Za、Hax、Hay、ΔT)顯示及三維立體圖(Za、Hax、Hay、ΔT)顯示兩方面，此功能主要是將計算出的磁異常三分量和總磁異常以圖像的形式顯示出來，使結果呈現更加形象、具體。為實現該功能，通過GUI編輯界面，選擇添加一個坐標軸對象，調整坐標軸對象的大小及位置。

圖2 自航磁性體三分量磁陣列探測模型

該坐標軸對象區域即為異常圖像顯示區域。

(3) 添加運動磁性體運動狀態顯示功能。為使得運動磁性體的運動狀態可以形象的顯示在仿真系統主界面上，同上步一樣，通過GUI編輯界面，添加一坐標軸對象，調整坐標軸對象大小及其位置。該坐標軸對象區域為運動磁性體運動狀態顯示區域。

(4) 添加參數設置功能。利用可編輯文本框實現數據的改變輸入，實現在不同物性參數下的模型研究，體現了良好的人機交互功能。該仿真系統中參數主要包括三方面：運動磁性體的物性參數、地磁場參數及圖像顯示選項，通過GUI編輯界面，添加可編輯文本對象。可編輯文本對象框內為數據輸入區域，再添加靜態文本對象，靜態文本框內內容主要為注釋功能，使程序界面通俗易懂，便于操作、識別。

(5) 添加實時數據顯示功能。為使用戶能具體了解磁異常的大小，需要將計算出的磁異常數據在仿真系統界面上顯示。通過GUI編輯界面，添加列表框對象，調整列表框對象的大小及其位置。計算出的磁異常數據將在列表框中動態顯示。

(6) 修改、添加、調整。修改添加的各個對象屬性及名稱，調整各個對象在GUI界面上的位置，使得界面整體更加美觀、規范。

水下自航磁性體三分量磁測主界面結構框圖見圖3;根據上面步驟完成的仿真系統會顯示于界面圖。

2.3 GUI回調程序的編寫

設置好各個對象的屬性后，進行回調程序的編寫，該步驟是實現圖形用戶界面程序功能的最關鍵步驟。

2.3.1 三分量磁異常計算程序的實現

水下自航磁性體的三分量磁陣列探測,編寫回調函數時，首先要清楚研究的內容與目的，對問題進行適當分析與問題簡化，在本例中即為利用一有限長水平圓柱磁性體代替了水下自航磁性體。

將上面式(1)～(4)編輯轉換為Matlab編程語言語句，完成斜磁化有限長水平圓柱體磁場的垂直分量、水平X分量、水平Y分量及總磁異常數值模擬計算功能，這是本例中的計算基礎，其中式(1)～(4)中的計算所需的變量參數是通過可編輯文本對象進行更改。編程中需注意矩陣運算中的運算符號的變化。本例中分別定義了四個自定義函數式(5)～(8)分別計算有限長水平圓柱磁性體引起的磁異常的垂直分量、水平X分量、水平Y分量及總磁異常的數值大小，主函數調用這四個自定義函數完成磁異常各個數值的計算。

圖3 主界面結構框圖

function Daan1=Za(x,y,R,L,r1,r2,Mx,My,Mz)

(5)

function Daan2=Hax(x,y,R,L,r1,r2,Mx,My,Mz)

(6)

function Daan3=Hay(x,y,R,L,r1,r2,Mx,My,Mz)

(7)

function Daan4= Del_T(x,y,M1,I,A,R,L,r1,r2)

(8)

其中：式(5)為計算有限長水平圓柱體引起的磁異常垂直分量具體程序內容如下：

function Daan1=Za(x,y,R,L,r1,r2,Mx,My,Mz)

%UNTITLED2 此處顯示有關此函數的摘要

% 此處顯示詳細說明

Daan1=(Mx*x*R/((x^2+R^2)^2))*(((y-L)*(3*(x^2+R^2)+2*(y-L)^2))/r1^3-((y+L)*(3*(x^2+R^2)+2*(y+L)^2))/r2^3)+(Mz/((x^2+R^2)^2))*((((y-L)*(((y-L)^2)*(x^2-R^2)-(x^2+R^2)*(2*R^2-x^2)))/r1^3)-(((y+L)*(((y+L)^2)*(x^2-R^2)-(x^2+R^2)*(2*R^2-x^2)))/(r2^3)))-My*R*((r2^3-r1^3)/((r1^3)*(r2^3)));

end

式(6)～(8)中具體編程語句由上同理編輯完成。

2.3.2 磁異常圖示結果可視化的實現

實驗中將三分量磁探頭放置在海底，組成運動磁性體探測磁陣列，當運動磁性體經過三分量探測磁陣列時，自航磁性體所引起的磁異常將被海底布置的三分量磁陣列磁探頭探測(如圖2)，在上位機(PC機)中顯示出探測的磁異常。

磁異常圖示結果可視化的實現即為將式(1)～(4)式中計算的結果以圖像的形式顯示，可分為磁異常2D平面等值線圖像和磁異常3D空間立體圖像顯示。磁異常2D平面等值線圖像可直觀地分辨出磁異常等值線的密集處與稀疏處，而磁異常3D空間立體圖像可形象地觀察出測區各處的磁異常變化情況，易觀察磁異常峰值處及梯度變化處。兩種磁異常圖像各有優點，在本仿真系統可根據用戶的選擇改變繪圖類型。

(1) 平面等值線圖實現部分。利用式(5)～(8)中數值模擬系統中計算出的結果，通過Matlab中繪圖函數，其中繪制磁異常2D平面等值線圖像利用函數contour()。contour()函數是Matlab中繪制等值線函數。利用仿真系統主界面中繪圖類型控制面板，用戶可以自由選擇繪制圖像類型及圖像繪制目標：磁性目標體的Za、Hax、Hay三分量或者總磁異常ΔT。其中繪制磁性目標體Za分量平面等值線圖程序如下：

contour(x,y,Za);

其中：變量x為X坐標軸坐標值；變量y為Y坐標軸值；Za即為該坐標處磁異常的垂直分量。

完成程序編寫后，點擊運行程序，輸入水平圓柱磁性體的半徑、長度、運動速度、地磁場的大小、磁傾角、磁偏角等數據，以一距海底50 m，半徑為6 m，長15 m，運動速度為50 m/s，方向與X軸呈45度的水平圓柱磁性體為例，設置地磁場為50 mT，磁傾角45°，磁偏角45°等數據。點擊正演按鈕得到結果如圖4，同理可以繪制Hax，Hay，ΔT的磁異常等值線圖，如圖5～7。由圖4～7可以看出,磁異常等值線最密集的處即為運動磁性體的位置所在，根據磁異常等值線密集處運動狀態變化則可以判斷運動磁性體的運動方向及速度大小。

(2) 三維立體圖實現部分。利用式(5)～(8)中數值模擬系統中計算出的結果，通過Matlab中繪圖函數函數surf()繪制磁異常三維立體圖像。surf()函數是Matlab中繪制表面圖函數。利用仿真系統主界面中繪圖類型控制面板，用戶可以自由選擇繪制圖像類型及圖像繪制目標：磁性目標體的Za、Hax、Hay三分量或者總磁異常?T。其中繪制磁性目標體Za分量3D立體圖的程序如下：

surf(x,y,Za);

完成程序編寫后，點擊運行程序，輸入相應參數，選擇繪制磁異常Za3D立體圖，得到圖8。同樣可以選擇繪制磁異常水平X、Y分量及總磁異常的3D空間立體圖，如圖9～圖11。

圖4 運動磁性體Za等值線圖

圖5 運動磁性體Hax異常平面等值線

根據不同參數繪制出的磁異常平面等值線圖及磁異常3D空間立體圖的特性，研究運動磁性體的位置和運動方向及速度，總結出不同狀態下的運動磁性體的磁異常圖像特征，為實際探測到的磁異常曲線圖解釋提供理論依據，為研究并推斷運動磁性體的走向長度、半徑、運動速度及方向等物性參數理論基礎。

在該仿真系統主界面上,用戶可以非常方便地輸入運動磁性體的物性參數和幾何參數,然后直觀地觀察不同參數下的運動磁性體的磁異常在整個平面和空間的分布特征,這些特征將為實際環境中探測到的運動磁性體的磁異常特征解釋推斷提供一個有力分析工具[5]。

2.3.3 仿真系統人機交互功能的實現

通過人機交互功能實現研究不同物性參數條件下的運動磁性體的磁異常特征。仿真系統主界面處設置有數據錄入控制面板，在此處可以改變運動磁性體的半徑、走向長度、運動速度、磁場大小、磁偏角、磁傾角等參數，根據用戶輸入的參數，實現研究磁性目標體在不同深度、不同半徑、不同長度、不同運動速度以及不同環境磁場條件下磁性體三分量磁異常的變化特征。

圖7 運動磁性體ΔT異常平面等值線

圖8 運動磁性體Za3D空間立體圖

圖9 運動磁性體Hax3D空間立體圖

下式為讀取用戶設置的運動磁性體半徑函數：

r=get(handles.edit1,'string');

R=str2num(r);%圓柱半徑

其中：handles.edit1為該可編輯文本框對象的對象名，利用get()函數讀取該文本框中字符，讀取的變量r為字符型，無法直接參與數值計算，則利用str2num()函數將字符型的變量r轉變成數值型的變量R，參與數值計算，其他物性參數同上利用同樣函數讀取，仿真系統主界面中的多個可編輯文本框中讀取的物性參數為數值模擬計算的基礎。

圖10 運動磁性體Hay3D空間立體圖

圖11 運動磁性體ΔT3D空間立體圖

2.3.4 磁性體自航時動態顯示的實現

為實現對運動中的磁性體磁異常特征的研究，雖然靜止磁性體的磁異常并不能滿足研究要求，但是運動的磁性體磁異常特征可以分解為多個按一定時間順序排列的靜止磁性體磁異常特征合成。根據這一研究思想，實現磁性體的運動顯示，具體程序代碼如下：

pausetime=0.01;

pause(pausetime);

drawnow

其中：變量pausetime是屏幕圖像顯示暫停時間變量，這里設置的是0.01 s;pause()函數是暫停函數;drawnow是刷新屏幕函數，這兩個函數共同使用，則可以實現對運動磁性體磁異常特征的研究。

2.3.5 自航磁性體磁異常數據輸出的實現

磁異常平面等值線圖像和3D立體圖像可以定性的研究磁異常的特征，但是當需要定量研究某些點磁異常的特征時，磁異常的平面等值線和3D立體圖像并不能很好的滿足，則需要直接將磁異常數據輸出。本仿真系統磁異常數據輸出的方法有兩種：① 將數值模擬計算結果，直接輸出至仿真系統主界面上實時顯示;② 將數值模擬計算結果保存在數據存儲文件(txt文件、dat文件等)中。在圖4仿真系統主界面中可以選擇是否將式(5)～(8)中數值模擬系統中計算出的結果導出至存儲文件中儲存，同樣也可以選擇磁陣列中某一探頭，將該點磁異常值顯示在GUI界面上，凸顯該探頭處在磁性體運動整個過程中磁異常數值變化。下面程序代碼是將磁異常顯示在數據顯示框中：

Lstring=num2str(Za);

inform{length(inform)+1,1}=Lstring;

set(handles.listbox2,'String', inform);

其中：num2str()函數是將Za由數值類型轉換為字符型的變量Lstring；handles.listbox2是顯示數據的列表對象名稱，利用set()函數將字符型的變量Lstring顯示在數值顯示框中。

下面為將數值模擬計算結果輸出并保存在數據文件TXT文件的函數。

[fn1,pn1]=uigetfile('*.txt','Select the Data file');

fid1=fopen(strcat(pn1,fn1),'wt');

fprintf(fid1,'%f %f %f %f %f %f ',x,y,Za,Hax,Hay,Del_T);

fclose(fid1);

其中：uigetfile()是選擇指定輸出數據文件夾函數；fopen()是打開文件夾函數；fprintf()是文件夾輸出函數；fclose()是關閉文件夾函數；fid1是文件夾指針變量；x、y、Za、Hax、Hay、Del_T分別是X坐標值、Y坐標值、磁異常垂直分量、磁異常水平X分量和磁異常水平Y分量；”%f”是數據輸出格式。

某一時刻水平圓柱磁性體的三分量數據見表1。

表1 海底三分量磁陣列某一時刻的的測量數據

T0=50 μT；A=90°；I=90°；V=50 m/s；R=6 m；L=15 m

3 結 語

基于Matlab中GUI程序開發環境，設計磁通門三分量磁陣列進行運動磁性體中的探測模擬，實現運動磁性體三分量磁異常可視化數值模擬計算功能，研究運動中的磁性體磁測中Za、Hax、Hay三分量及總磁異常的變化規律。該仿真系統可根據用戶對不同深度、不同半徑、不同長度、不同運動方向和速度的磁性目標體以及不同地磁場環境進行任意參數設置，從而可以更好地反映真實環境中運動磁性體的磁異常特征。采用磁異常信號對水下運動磁性體目標探測技術能彌補聲吶探測在這方面的不足，更好地解決“隱身”運動磁性體難以探測的問題。

相比于C語言、VB語言等計算機程序語言，Matlab擁有更加強大的數學運算能力及圖形處理能力，程序語言簡便，功能十分強大，可根據用戶要求自行設計美觀合理的界面。利用強大的Matlab程序設計語言，本文實現了運動磁性體三分量磁異常的動態顯示和其運動的動態顯示及良好的人機交互能力，使得可以更加全面地定性、定量分析磁異常的特征。

[1] 張勝業,潘玉玲.應用地球物理學原理[M].武漢:中國地質大學出版社,2004:23-43.

[2] 管志寧.地磁場與磁力勘探[M].北京:地質出版社,2005:23-45.

[3] 胡祥超.基于磁異常信號的目標探測技術試驗研究[D].長沙.國防科技大學,2005.

[4] 王林松，王傳雷，西永在，等.球體磁力異常三分量數據空間成像[J].工程地球物理學報.2009,6(4):404-410.

[5] 張 劍,師學明.基于Matlab開發環境的球體重力正演[J].工程地球物理學報,2007,4(5):460-464.

[6] 劉衛國.Matlab程序設計教程[M].北京:水利水電出版社,2005:32-49.

[7] 范師美. 斜磁化有限長水平圓柱體的磁場[J]. 物化探電子計算技術,1980(1):75-84.

[8] Kim W H, Kim J C, Choi I K,etal. Improvement of the signal to noise ratio of the magnetic detection system using the geological magnetic filter[C]//Electronics, Circuits and Systems, 1999. Proceedings of ICECS'99. The 6th IEEE International Conference on. IEEE, 1999(1): 237-240.

[9] Caruso M J, Bratland T, Smith C H,etal. A new perspective on magnetic field sensing[J]. SENSORS-PETERBOROUGH-, 1998, 15: 34-47.

[10] Caruso M J, Withanawasam L S. Vehicle detection and compass applications using AMR magnetic sensors[C]//Sensors Expo Proceedings. 1999: 477.

[11] Wiegert R. Magnetic anomaly guidance system for mine countermeasures using autonomous underwater vehicles[J]OCEANS 2003. Proceedings. IEEE, 2003(4): 2002-2010.

[12] Charles L. Bernstein; Close-range sensors for small unmanned bottom vehicles[C]∥Proc. SPIE 3711, Information Systems for Navy Divers and Autonomous Underwater Vehicles Operating in Very Shallow Water and Surf Zone Regions, 92 (July 23, 1999); doi:10.1117/12.354644.

[13] 葛哲學.精通Matlab[M].北京.電子工業出版社,2008:23-65.

[14] 史旭華，俞海珍. 基于工控組態軟件及 Matlab 的計算機控制實驗平臺[J].工業控制計算機,2008,21(6):14-15.

[15] 王印松，岑 瑋. 基于Matlab/Simulink 電力系統仿真工具箱的拓展[J].電力系統保護與控制,2009，36(9):84-87.

[16] 鄭黎明，孫敬偉. 基于Matlab/Simulink的新型混合SVPWM的分析[J].電機與控制應用,2009，36(9):4-7.

[17] 尚 麗，淮文軍.基于Matlab/Simulink和GUI的運動控制系統虛擬實驗平臺設計[J].實驗室研究與探索,2010,29(6)：66-71.

[18] 李建中.強大的系統分析工具MATLA B[J].陰山學刊,2001,6(5):67-69.

[19] 陳義群,陳 華.基于Matlab的工程物探軟件快速開發[J].地球物理學進展,2004,19(4):802-806.

[20] 程 鈴,徐冬冬.Matlab仿真在通信原理教學中的應用[J].實驗室研究與探索,2010,29(2):117-119.

Detective Simulation System of the Self-propelled Object’s Three-component Magnetic Array Based on Matlab

SHENXiao-bing，LIYong-tao

(Institute of Geophysics Geomatics, China University of Geosciences, Wuhan 430074, China)

Based on GUI procedure development circumstance in Matlab, we develop a software which simulates magnetic flux gate array of three-component magnetic to survey underwater self-propelled object's movement. The simulated effects are as follows: it can make self-propelled object visible when the three-component is abnormal; it can imitate computing functions; it can also provide feature analysis of Za, Hax, Hay and the total magnet when the magnetic target is moving. According to self-propelled feature of magnetic body movement, the motion situation and magnetic anomaly curve can be observed. According to the parameter inputted by user, it can actualize simulation of magnetic object's anomalous conditions in different depths, semidiameters, lengths, speeds, and magnetic surroundings. Through the simulation experiment, it can deepen students' comprehension about the skill of detection of magnetic object underwater and give certain reference to student those who engage in geophysical exploration and military geophysical.

Matlab; three-component; magnetic array; simulation system; self-propelled target body; magnetic anomaly

2015-03-13

湖北省科技支撐計劃資助項目(2015BCE077)；中國地質大學(武漢)2014-2015學年度教學實驗室開放基金資助項目

沈小兵(1993-)，男，江西景德鎮人,本科生，就讀于地球物理學專業。Tel.：13387590418；E-mail：1136926533@qq.com

李永濤(1961-)，男，江蘇南通人，博士，教授，主要從事地球物理學及環境磁學研究。

Tel.：027-67883447；E-mail：ytli@cug.edu.cn

P 631.2

A

1006-7167(2015)11-0103-08