基于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