丁学振，李卓轩，李予国,2,3*，刘浩
1 中国海洋大学海洋地球科学学院，山东青岛 266100 2 中国海洋大学海底科学与探测技术教育部重点实验室，山东青岛 266100 3 青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266237
磁法探测在地质调查、资源探测、军事、环境等领域有重要的应用价值（管志宁，1997；
林君等，2017）.磁法探测的发展可分为总磁场测量、总场梯度测量、矢量场和磁梯度张量场测量等3个阶段（张昌达,2006）.相对于磁场数据，磁梯度张量数据具有异常信号更敏感、空间分辨率高、降低地磁噪声、减少对基站的需求等优势（Pedersen and Rasmussen,1990;Schmidt and Clark,2006;Clark,2012）.
磁梯度张量是磁场矢量在3个正交方向上的空间变化率（Schmidt et al.,2004）.在实际测量中，可以采用差分近似法、旋转调制法（Tilbrook,2009）和直线测量法（Sunderland et al.,2009）测得磁梯度张量数据.差分近似法易于实现，并且能够获得全张量信息，因此得到了广泛应用（苗红松,2017）.目前，基于差分近似原理的磁梯度张量数据采集系统大致可分为基于超导量子干涉仪的磁梯度张量系统（Schmidt et al.,2004;Gamey,2008;Keenan et al.,2010;汪瀛等,2015）和基于磁通门传感器的磁梯度张量系统（Koch et al.,1996;Sulzberger et al.,2006;Sui et al.,2014;Yin et al.,2014）.磁梯度张量超导量子干涉系统具有很高的灵敏度，但量程范围较小、成本高，并且对系统安装工艺和测量环境要求也较高.而基于磁通门传感器搭建的磁梯度张量系统成本较低，对安装工艺要求相对较低，利于大批量生产制造（李青竹等,2017）.由于磁通门传感器制作工艺和磁梯度张量系统安装精度的限制，张量测量精度会受到单个传感器误差（零偏、三轴非正交误差、三轴灵敏度不一致）和多传感器非对准误差的影响，磁梯度张量系统输出误差有时可达上千nT/m （李青竹等,2018），因此必须对其进行校正.
目前，大部分磁梯度张量系统校正方法为两步校正方法，该方法包括两个步骤（Pang et al.,2013a,b;Yin et al.,2014,2015a;迟铖等,2017）.第一步是对单传感器的误差进行校正.该校正可分为矢量校正方法（Pang et al.,2014）和标量校正方法（Pang et al.,2013b;Yin et al.,2014,2015a;迟铖等,2017）.矢量校正方法以高精度磁场矢量为标准对磁通门传感器进行校正，该方法需要测量准确的磁场矢量，对试验条件要求较高.标量校正方法以高精度磁场总场为标准对磁通门传感器进行校正，该方法因较为方便且成本低而得到广泛应用.第二步是对磁梯度张量系统的非对准误差进行校正.以磁梯度张量系统中的一个传感器或搭载平台坐标系（Yin et al.,2014;Li et al.,2018）为基准，将多个传感器校正至同一坐标系中.两步校正方法的缺点在于第二步的校正结果依赖于第一步校正结果的精度，第一步校正结果的误差直接影响着第二步校正效果.为了解决两步校正方法存在的问题，Yin等（2015b）和Li等（2018）提出了基于最小二乘拟合方法的一步校正方法.一步校正方法将磁梯度张量系统的理想输出校正至搭载平台坐标系中，需要令张量系统绕三个正交的坐标轴旋转采集数据，该方法依赖于高精度的三轴无磁旋转台且试验难度较大.上述两步校正方法和一步校正方法大多基于最小二乘拟合方法估计磁梯度张量系统的误差参数，而最小二乘拟合反演方法对初始参数比较敏感，为获得较好的校正结果，需对初始参数进行仔细调整（Pang et al.,2013c）.
为了进一步解决一步和两步校正方法存在的问题，本文提出了一种基于改进差分进化算法的磁梯度张量系统校正方法.差分进化算法具有全局寻优、不受初始模型选取限制等优点，可以有效解决基于最小二乘拟合的校正方法依赖于初始参数的问题.我们将张量分量和张量不变量引入反演目标函数中，只需一步反演即可实现单传感器误差和张量系统非对准误差的校正，有效避免了两步校正方法中第一步校正精度对第二步校正结果的影响问题，提高了校正精度.
1.1 张量数据采集系统
磁梯度张量系统包括十字形结构、正方形结构、三角形结构、直角四面体和正四面体结构等5种基本结构（刘丽敏,2012）.磁梯度张量系统采用差分近似的原理测量张量数据，不同结构的张量系统会带来测量误差，其中，十字形结构的磁梯度张量系统测量误差最小.图1为十字形磁梯度张量系统示意图，它由4个三轴磁通门传感器组成，其中1号和3号传感器沿x轴分布，2号和4号传感器沿y轴分布，基线距离为d.根据差分近似原理，磁梯度张量系统中心位置o处的张量矩阵为（Li et al.,2018;Yin et al.,2015b）
（1）
式中，Bij表示第i（i=1,2,3,4）个传感器测得的j（j=x,y,z）方向上的磁场分量.在磁法测量中，通常将磁性物体产生的异常场看作为无源的静磁场，异常场的散度和旋度都为零，故张量矩阵G中只有Bxx,Bxy,Bxz,Byy,Byz5个独立分量（尹刚等,2016）.
由张量矩阵G的特征值可以计算得到磁梯度张量不变量，三个张量不变量的表达式如下（Clark,2012;Mu et al.,2019）：
（2）
式中，λ1,λ2,λ3为张量矩阵G的3个特征值，且满足λ1≥λ2≥λ3,|λ1|≥|λ2|,|λ3|≥|λ2|.
图1 十字形磁梯度张量系统示意图Fig.1 Schematic diagram of the cross-shaped magnetic gradient tensor system
1.2 误差校正
在理想情况下，三轴磁通门传感器的三个轴保持两两完全正交，但是由于制作工艺的限制，传感器普遍存在非正交误差.建立传感器非正交模型如图2所示，O-XYZ为标准正交坐标系，O-X1Y1Z1为传感器实际坐标系，O-X2Y2Z2为理想正交坐标系.假设OZ1和OZ2共轴，平面Y1OZ1和平面Y2OZ2共面，则OY1和OY2之间的夹角为ψ，OX1与平面X2OY2的夹角为φ，OX1在平面XOY的投影OX′1与OX2的夹角为θ.一旦非正交角φ,θ,ψ确定后，传感器理想正交坐标系O-X2Y2Z2即可唯一确定.
图2 单传感器非正交误差示意图Fig.2 Schematic diagram of the non-orthogonal error for single sensor
在理想情况下，磁通门传感器在无磁环境下输出为零，但是由于磁芯存在剩磁或者电路部分存在零位偏移导致传感器三轴输出不为零.一般情况下，传感器的零偏可达几十nT.另外，由于传感器磁敏元件制作工艺的限制，可能会导致三轴灵敏度不一致，从而使得三轴在相同磁场环境下的测量值不同.传感器三轴实际测量值与真实值的比值称为灵敏度标度因子（李青竹等,2017）.假定传感器三轴零偏为b=（bx,by,bz），三轴灵敏度标度因子为kx,ky,kz.综合考虑零偏、灵敏度不一致误差和非正交误差，建立的单传感器误差模型为（Li et al.,2018）：
Bc=KCB+b，
（3）
在实际情况下，传感器的实际输出Bc是已知的，传感器误差的校正就是由其实际输出求取理想输出B，即
B=（KC）-1（Bc-b）.
（4）
由于传感器安装精度的限制，构成磁梯度张量系统4个矢量传感器的测量轴指向可能不一致，张量系统非对准误差可通过旋转方式进行校正（Pang et al.,2013b;Li et al.,2018）.定义传感器绕x轴旋转为横倾角α，绕y轴旋转为俯仰角β，绕z轴旋转为方位角γ.如图3所示，O-X1Y1Z1表示1号传感器的正交坐标系，O-XiYiZi表示i（i=2,3,4）号传感器的正交坐标系.假定i号传感器先绕y轴旋转，再绕x轴旋转，最后绕z轴旋转，可以令i号传感器坐标系与1号传感器坐标系重合，即可实现非对准误差校正.
图3 传感器非对准误差校正示意图Fig.3 Schematic diagram of misalignment error calibration
由式（4）可知，单传感器误差校正后的输出为B，经非对准误差校正后得到的输出Br可表示为（迟铖等,2017）
Br=TγTαTβB，
（5）
综合考虑单传感器误差和多传感器非对准误差，可得
Br=TγTαTβ（KC）-1（Bc-b）.
（6）
由上式可知，如果得到单传感器误差参数和多传感器非对准误差参数即可实现张量系统误差校正.
差分进化算法是一种模拟自然界生物遗传和进化机理发展起来的基于种群的全局优化算法，它是由Storn和Price于1995年提出的（Storn and Price,1995）.差分进化算法具有原理简单，控制参数较少，鲁棒性较好的优点（Storn and Price,1995,1997），已得到较为广泛的应用（Das and Suganthan,2011）.差分进化算法主要包括参数初始化、变异、交叉和选择4个步骤.参数初始化就是在预先设定的搜索空间中随机产生NP个维度为P的参数向量，在该步骤中随机产生的参数向量应尽可能覆盖整个搜索空间.变异操作的目的是产生新的个体以增加种群多样性，当前种群向量称为目标向量，通过使用某种变异策略对目标向量进行变异操作产生新的变异向量.在交叉操作中，令目标向量和变异向量随机交换参数，产生新的试验向量.选择的目的是保持种群数量不变，使用“贪婪性”选择策略将优良个体保留至下一代种群中.下面，结合磁梯度张量系统误差校正的目的介绍差分进化算法的实现步骤.
（1）参数初始化
在参数初始化这一步骤中，随机产生NP个维度为P的参数向量.针对磁梯度张量系统误差校正这一应用，每个传感器有12个误差参数，于是参数向量的维度为48.假设当前进化代数为g，第m个体可以表示为
m=1,2,…,NPi=1,2,3,4
（7）
式中，φ,θ,ψ表示传感器的非正交角（图2），kx,ky,kz表示三轴灵敏度标度因子，bx,by,bz表示零偏误差，α,β,γ表示非对准角（图3），i表示传感器序号.
在参数初始化过程中，应尽量保证随机产生的参数向量能够覆盖整个搜索范围.假设参数向量的取值上下界为Xmax=（x1,max,x2,max,…,xP,max）和Xmin=（x1,min,x2,min,…,xP,min），则初始化过程可表示为
（2）变异操作
（9）
式中，F为处于[0,2]之间变异因子（Storn and Price,1997），r1,r2,r3为[1,NP]之间随机选择的不同于m的整数.变异策略可以表示为DE/X/Y/Z的形式，其中DE表示差分进化算法，X表示指定的变异向量，Y表示变异过程中使用的差分向量个数，Z表示交叉方式（Das and Suganthan,2011）.
（3）交叉操作
（10）
（4）选择操作
为了使下一代中种群数量保持不变，采用“贪婪性”策略将试验向量与目标向量进行对比，优良个体将会保留至下一代种群中，选择操作可以表示为
（11）
式中，f（x）为目标函数，其值越小表明个体表现越好.
在传统的差分进化算法中，变异因子F和交叉概率因子CR均为固定值，且只采用DE/rand/1/bin变异策略.为了提高算法的搜索效率和精度，避免陷入局部最优，本文使用DE/rand/1/bin、DE/target-to-best/2/bin变异策略（陈亮,2012）、时变变异因子（颜学峰等，2006）和时变交叉概率因子（王天意,2015）对差分进化算法进行了改进.
改进差分进化算法迭代过程中变异策略选取方法为
（12）
式中，rand（0,1）为处于[0,1]之间的随机数，φ为随迭代次数变化的阈值，可以表示为
（13）
其中，φmax=1,φmin=0.1分别为φ取值上下界（陈亮,2012），itermax和iter分别为最大迭代次数和当前迭代次数.
时变变异因子可以表示为（颜学峰等，2006）
（14）
式中，F0为初始变异因子.
时变交叉概率因子可以表示为（王天意,2015）
（15）
式中，CR0为初始交叉概率因子.
如果磁梯度张量系统不存在误差，在匀强磁场环境下旋转其姿态测量数据，由磁通门传感器三轴读数计算得到的总磁场为一定值，且不同传感器相同方向的测量轴读数一致.磁梯度张量系统误差校正标准为
（16）
式中，Bxi,Byi,Bzi（i=1,2,3,4）分别为由每个个体向量计算得到的第i个传感器的三分量数据，Bti为由三分量数据计算得到的总磁场，Bc为参考总磁场.
本文针对磁梯度张量系统校正目的，设计如下目标函数：
（17）
由式（17）可知，目标函数包括两项：第一项为由4个磁通门传感器3分量数据构成的总磁场与质子磁力仪测得的参考总磁场之间的偏差，它主要用于校正单传感器误差（零偏、三轴灵敏度不一致和三轴非正交误差）.第二项为传感器三分量数据的偏差、张量分量和张量不变量，主要用于校正磁梯度张量系统的三轴非对准误差.
3.1 仿真试验1
为验证本文方法的有效性，设计了仿真模拟试验.如表1所示，磁梯度张量系统零偏误差在-200～200 nT之间随机产生，非正交误差和非对准误差在-2°～2°之间随机产生，灵敏度参数在0.9～1.1之间随机产生.假定磁梯度张量系统的基线距离为0.4 m，试验地点处总磁场为52000 nT，磁偏角和磁倾角分别为-7°和54°，使用磁梯度张量系统在三维空间内随机旋转采集200组数据.
在实际应用中，磁梯度张量系统采集数据可能会存在一定的测量误差，且不可避免地受到非均匀背景场的影响.因此，在仿真试验中加入随机噪声，在4个磁通门传感器的3个测量轴中分别加入了最大幅值为0 nT、2 nT、5 nT和10 nT的随机噪声.加入不同大小噪声时，反演得到的磁梯度张量系统误差参数如表1所示，校正后的总磁场和张量分量如图4所示.当仿真模拟数据不含随机噪声时，反演得到的误差参数和预设参数完全一致，校正后的总磁场和张量分量分别为52000 nT和0 nT/m.当仿真模拟数据加入2 nT、5 nT和10 nT的随机噪声时，反演得到的零偏误差参数与预设参数的最大偏差分别为0.367 nT、0.909 nT和1.825 nT，非正交误差与预设参数的最大偏差为0.001°、0.003°和0.006°，对灵敏度不一致误差和传感器非对准误差6个误差参数的反演结果没有影响.由图4可知，当加入最大幅值为10 nT的随机噪声时，校正后的张量分量仍小于0.2 nT/m.
表1 预设误差参数与反演得到的误差参数Table 1 Preset and estimated error parameters
图4 加入不同大小随机噪声得到的总磁场和张量分量校正结果Fig.4 Calibration results of total magnetic intensity and tensor component
图5 不同种群大小得到的校正数据RMS随迭代次数的变化情况Fig.5 RMS error of calibration results with different population size
差分进化算法的反演精度与设定的种群数量NP、变异因子F和交叉概率因子CR有关（Balkaya et al.,2017）.本文分析了这些参数对反演结果的影响.将种群数量分别设定为10、15、20、50、100和200，分别反演200组不含噪声的磁场数据.令最大迭代次数为10000，利用每次迭代中最佳个体校正磁场数据，并计算校正数据的RMS.RMS计算公式如式（18）所示.图5为不同种群数量时校正数据RMS随迭代次数的变化情况，可以看出随着种群数量增大，差分进化算法陷入局部最优解的可能性降低.当种群数量为50时即可实现精度较高的反演.
RMS=
（18）
式中，N为磁场数据量.
假定变异因子和交叉概率因子的取值范围分别为[0.2,0.9]和[0.2,1]，我们分别采用传统差分进化算法和改进差分进化算法估计张量系统误差参数，并计算校正后总磁场与张量分量RMS的和，结果如图6所示.假定校正后总磁场与张量分量RMS之和小于0.01时校正精度是可以接受的.由图6可知，与传统差分进化算法相比，改进差分进化算法对变异因子和交叉概率因子依赖性较小.本文算例中F0取值为0.3，CR0取值为0.8.
参与反演的数据量可能也会影响误差参数的反演精度，为分析数据量对反演结果的影响，令磁梯度张量系统在三维空间中随机旋转采集10组、20组、50组、100组、200组和500组磁场数据，并将其输入改进差分进化算法反演误差参数.令种群数量为200，最大迭代次数为10000，利用每次迭代中最佳个体校正磁场数据，并计算校正数据的RMS.由图7可知，当数据量为10组时，无法得到良好的校正效果，而当数据量为50组时即可实现误差参数的精确反演.
将本文方法所得结果与遗传算法、最小二乘拟合方法的校正结果进行了对比.在反演中使用了200组磁场数据，最大迭代次数为10000，种群数量分别设置为50和200，利用每次迭代中最佳个体校正磁场数据，并计算校正数据的RMS.改进差分进化算法和遗传算法反演过程中RMS随迭代次数变化情况如图8所示.由图8可知，改进差分进化算法的搜索效率与反演精度明显高于遗传算法.当种群数量为50时，改进差分进化算法迭代次数约为2000时反演得到的RMS收敛于10-4量级，而遗传算法迭代10000次得到的RMS仍为17.68.
校正前后磁场数据的均方根误差如表2所示.由表2可知，当不存在随机噪声时，本文方法和最小二乘方法的校正结果相同，均能准确反演出磁梯度张量系统的误差参数.添加随机噪声后，本文方法反演结果的RMS明显小于最小二乘方法，本文方法具有较强的抗噪能力.与遗传算法相比，本文方法反演结果的RMS小于遗传算法反演结果的RMS，这意味着本文方法的反演精度较高.
表2 磁场数据校正前后的均方根误差Table 2 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration
图6 两种差分进化算法校正后磁场数据的均方根误差（a） 传统差分进化算法；
（b） 改进差分进化算法.Fig.6 RMS error of calibration results with differential evolution algorithm（a） Traditional differential evolution algorithm；
（b） Improved differential evolution algorithm.
图7 不同数据量得到的校正数据RMS随迭代次数的变化情况Fig.7 RMS error of calibration results with different magnetic data
图8 改进差分进化算法和遗传算法反演结果对比Fig.8 Comparison of RMS error of calibration results with improved differential evolution algorithm and genetic algorithm
3.2 仿真试验2
为进一步测试本文方法的实用性和可靠性，依据某国产磁通门传感器测定的误差大致取值范围，设计了仿真模拟试验.如表3所示，磁梯度张量系统的零偏误差在-100～100 nT之间随机产生，非正交误差在-0.1～0.1度之间随机产生，灵敏度参数在0.9999～1.0001之间随机产生，非对准误差在-1～1度之间随机产生.同样，假定磁梯度张量系统的基线距离为0.4 m，试验地点处总磁场为52000 nT，磁偏角和磁倾角分别为-7°和54°，使用磁梯度张量系统在三维空间内随机选择采集200组数据，在磁场数据中加入了最大幅值为0 nT、2 nT、5 nT和10 nT的随机噪声.
由表4可知，当仿真模拟数据不含随机噪声时，磁梯度张量系统的误差参数可以得到精确的反演.当仿真模拟数据加入2 nT、5 nT和10 nT的随机噪声时，反演得到的零偏误差参数与预设参数的最大偏差分别为0.339 nT、0.844 nT和1.715 nT，非正交误差与预设参数的最大偏差为0.001°、0.003°和0.006°，灵敏度误差与预设参数的最大偏差为2.2×10-5、5.6×10-5和1.13×10-4.由表5可知，经过校正后，磁梯度张量系统的测量误差得到了较大程度地压制，当仿真模拟数据加入10 nT的随机噪声时，校正后总磁场、张量分量Bxx,Bxy,Bxz,Byy,Byz的均方根误差分别为7.327、0.379、0.211、0.151、0.044和0.075.
表3 预设误差参数与反演得到的误差参数Table 3 Preset and estimated error parameters
表4 磁场数据校正前后的均方根误差Table 4 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration
为验证本文方法的实用性，在地磁干扰比较小的夜间开展了试验工作.试验中将十字形磁梯度张量系统安装于三轴无磁旋转台上（图9a）随机旋转采集数据，采用质子磁力仪（图9c）测量试验场地的地磁场.磁梯度张量系统基线距离为0.16 m，质子磁力仪分辨率为0.1 nT、测量精度为1 nT，三轴无磁旋转台旋转角度可读分辨率为1°.为降低试验中地磁场的天然变化对试验结果的影响，本次试验中只采集了60组数据，以缩短试验所需时间.
图9 试验所用仪器（a） 张量系统与无磁三轴旋转台；
（b） 张量系统记录仪；
（c） 质子磁力仪.Fig.9 Test equipment（a） Magnetic gradient tensor system and three-axis turntable；
（b） Data acquisition system；
（c） Proton magnetometer.
利用本文提出的校正方法对实测数据进行了校正，校正前和校正后的总磁场和张量分量如图10和图11所示.磁梯度张量系统误差未校正前测量值波动较大，地磁场总磁场的波动值可达数百 nT，磁梯度张量分量的最大偏差可达1000 nT/m，这说明磁梯度张量系统的误差会对测量结果产生较为严重的影响，必须对其进行校正.从图10和图11中可以看出，经过本文方法校正后，总磁场均收敛至试验地点处地磁场（51600 nT）附近，张量分量均收敛至0 nT/m附近.
图10 总磁场校正前后对比Fig.10 Comparison of total magnetic intensity before and after calibration
图11 张量分量校正前后对比Fig.11 Comparison of tensor components before and after calibration
校正前后总磁场与张量分量的均方根误差如表5所示.由表5可以看出，经过校正后，磁场数据的均方根误差明显减小，且本文方法的校正结果明显优于基于遗传算法和最小二乘方法的误差校正方法.以总磁场为例，校正前总磁场均方根误差为179.566 nT，基于本文方法、遗传算法和最小二乘方法的总磁场均方根误差分别为3.398 nT、18.197 nT和14.210 nT.
表5 磁场数据校正前后的均方根误差Table 5 Comparison of RMS error of total magnetic intensity and tensor components before and after calibration
磁梯度张量系统的测量精度受到磁场传感器零偏、灵敏度不一致、非正交和传感器非对准误差的影响.本文提出了一种基于改进差分进化算法的磁梯度张量系统误差校正方法，只需一步反演即可实现磁梯度张量系统48个误差参数的估计.用两个仿真模拟算例验证了本文方法的有效性，仿真结果表明本文方法与遗传算法相比具有更高的求解精度和搜索效率，本文方法较最小二乘拟合方法具有更高的抗噪能力.改进的差分进化算法能有效避免传统差分进化算法搜索效率低、易陷入局部最优解、对变异因子和交叉概率因子取值依赖大等问题.实测数据总磁场、张量分量Bxx,Bxy,Bxz,Byy,Byz校正后RMS改善比分别为52.845、54.530、25.415、40.057、38.330和25.586，本文方法能够有效降低磁梯度张量系统误差引起的测量误差，提高测量精度.本论文的不足在于实际试验场地磁场环境不够理想，导致实测数据总磁场校正后的均方根误差仍为3.398 nT.在下一步的工作中需寻找环境干扰更小的场地进行试验，以更好地估计磁梯度张量系统的误差参数.
致谢三位审稿人提出了建设性意见，对本文的质量提升帮助很大，在此表示诚挚的感谢.
猜你喜欢 张量差分反演 RLW-KdV方程的紧致有限差分格式数学杂志(2022年5期)2022-12-02反演对称变换在解决平面几何问题中的应用中等数学(2022年5期)2022-08-29符合差分隐私的流数据统计直方图发布湘潭大学自然科学学报(2022年2期)2022-07-28一类张量方程的可解性及其最佳逼近问题 ①西南师范大学学报（自然科学版）(2022年1期)2022-03-02基于ADS-B的风场反演与异常值影响研究成都信息工程大学学报(2021年5期)2021-12-30利用锥模型反演CME三维参数空间科学学报(2021年4期)2021-08-30数列与差分新世纪智能(数学备考)(2021年5期)2021-07-28严格对角占优张量的子直和华南师范大学学报（自然科学版）(2021年3期)2021-07-03四元数张量方程A*NX=B 的通解五邑大学学报（自然科学版）(2020年4期)2020-12-09一类麦比乌斯反演问题及其应用中等数学(2020年2期)2020-08-24