0 引言
原子磁力仪在现代科技应用中扮演着重要角色，其广泛应用涵盖了地球物理学、航空航天、目标探测与定位等多个行业领域[1-3]，尤其是在空间探测地球磁场和其他行星磁场中具有重要作用。目前，标量原子磁力仪主要包括光泵磁力仪、Overhauser磁力仪和相干布居囚禁（Coherent Population Trapping, CPT）原子磁力仪。光泵磁力仪灵敏度高，但机理上存在方向盲区和转向差问题，因此工程上光泵磁力仪需采用复杂的多探头组合及安装调试方法来解决方向盲区问题，多探头设计导致探头体积庞大，多探头数据切换时产生的数据干扰也会带来负面影响[4-5]。Overhauser磁力仪具有绝对精度高、无方向盲区的优点，但其采样率较低，一般不超过5 Hz，采样率较低会降低磁感应强度测量的空间磁测分辨率[6]。CPT原子磁力仪具有高灵敏度、高采样率、体积小和无方向盲区等优点，在高精度地磁图测绘和矿产资源勘探等领域具有巨大的应用前景[7-8]。
但是由于非线性塞曼效应，CPT原子磁力仪探头激光方向与外界磁场方向的夹角不同时会产生不同大小的测量误差，即转向差[9]。CPT原子磁力仪转向差曲线理论上在0~2π范围内是关于π轴对称的，如果已知激光方向与外界磁场方向夹角，则可以直接从磁场测量结果中补偿转向差，提高磁场测量精度。尽管探头内采取激光器无磁封装或导线磁屏蔽处理，但仍存在剩磁干扰，影响转向差的对称性和补偿精度。因此在转向差补偿前仍需先去除或减小剩磁产生的干扰，然后使用对称的转向差曲线进行补偿。
本文针对CPT原子磁力仪内剩磁干扰破坏理论转向差曲线的对称性，从而影响转向差的补偿精度问题，提出了一种基于粒子群算法估算剩磁参数和使用三次样条插值的转向差补偿方法。第1节介绍了CPT原子磁力仪转向差产生机理，建立了剩磁补偿模型与标定方法。第2节提出了具体的转向差补偿方法，包括粒子群算法和三次样条插值。第3节从试验的角度验证了所提出的转向差补偿方法，并对探头剩磁进行了估计，分析了转向差补偿结果。第4节为结论。
1 CPT原子磁力仪的转向差分析和标定
1.1 CPT原子磁力仪的转向差机理
CPT原子磁力仪的转向差来源主要有两个方面：一是原子能级的非线性塞曼效应引起CPT信号频率发生偏移导致的转向差，这部分为物理转向差；二是探头内部件材料中含有剩磁，其不对称性引入转向差。其影响虽然经弱磁化设计可降低，但是难以完全消除，尤其是距离原子气室比较近时，其影响不能忽略。
根据转向差产生的机理，可以建立磁感应强度测量的转向差模型，如图1 所示。待测磁感应强度为${{\boldsymbol{B}}_{\mathbf{e}}}$ ，探头轴向与磁感应强度方向的夹角为$\theta $ ，磁力仪转向差函数为${B_h}\left( \theta \right)$ ，${B_h}\left( \theta \right)$ 是与夹角$\theta $ 有关的标量。未知剩磁${{\boldsymbol{B}}_{\mathbf{r}}}$ 会影响磁感应强度测量结果，由于探头内光源和探测器的材料相对确定，光源驱动电流变化也很小，可以假定${{\boldsymbol{B}}_{\mathbf{r}}}$ 为固定磁场，即大小方向相对于探头均不变。因此，CPT原子磁力仪标量场输出$ B\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 可以表示为：
若在地磁场条件下，$\left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| \leqslant \left| {{{\boldsymbol{B}}_{\mathbf{e}}}} \right|$ ，可以忽略二次项，$\left| {{{\boldsymbol{B}}_{\mathbf{e}}} + {{\boldsymbol{B}}_{\mathbf{r}}}} \right|$ 可以写成：
式中，约等号右侧为待测磁感应强度大小，和剩磁在待测磁感应强度方向的投影。
1.2 CPT原子磁力仪的转向差标定
CPT原子磁力仪探头生产制造完成后，转向差作为传感器的固有属性，其相对于角度$\theta $ 的关系是恒定的，因此可以通过试验标定和数据处理得到$ B_{\mathrm{h}}\left(\theta\right) $ 与角度$\theta $ 之间的函数，随后用于实际磁探测中转向差补偿。试验标定过程中，设定$\left| {{{\boldsymbol{B}}_{\mathbf{e}}}} \right|$ 为已知量，由于探头存在剩磁的干扰，根据式（1）和式（2）可知，标定只能得到综合转向差函数$ B\mathrm{_H}\left(\theta,{\boldsymbol{B}}_{\mathbf{r}}\right) $ ：
从式（3）可知，${{\boldsymbol{B}}_{\mathbf{r}}}$ 的大小和方向均会影响$ B\mathrm{_H}\left(\theta,\boldsymbol{B}_{\mathbf{r}}\right) $ 的标定结果，这使转向差$ B\mathrm{_h}\left(\theta\right) $ 与$\theta $ 的关系难以确定。因此，研究剩磁估计方法，从$ B\mathrm{_H}\left(\theta,\boldsymbol{B}_{\mathbf{r}}\right) $ 数据中补偿剩磁${{\boldsymbol{B}}_{\mathbf{r}}}$ 对$ B\mathrm{_h}\left(\theta\right) $ 的干扰是后续的一项工作。
转向差的试验装置如图2 所示， CPT原子磁力仪探头和无磁码盘放置在稳定场发生装置内。该装置可以精确控制待测磁感应强度的大小和方向，电子箱与计算机放置在离磁场发生装置足够远的位置，以减小探头对磁感应强度测量的电磁干扰。设定待测磁感应强度${{\boldsymbol{B}}_{\mathbf{e}}}$ 方向为垂直方向，且与固定的无磁码盘0°刻度重合。探头可以通过机械结构在码盘中心旋转，外接计算机用于数据采集。具体标定方法如下：从探头与0°刻度重合开始，顺时针旋转探头并采集CPT原子磁力仪输出数据，探头轴向与${{\boldsymbol{B}}_{\mathbf{e}}}$ 方向之间的夹角$\theta $ 可以通过码盘直接读数。按照式（3）计算可以得到综合转向差$ {B_\mathrm {H}}\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 关于$\theta $ 的函数关系。
1.3 转向差的对称性
从转向差标定的方案可知，${B_\mathrm{h}}\left( \theta \right)$ 和$ {B_\mathrm{H}}\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 均是以2π为周期的函数。当探头处于如图3 所示的两种状态时，以顺时针旋转方向为正，状态1与0-π轴线的夹角为$\theta $ ，状态2与0-π轴线的夹角为$ {\text{-}}\theta $ 。理论上，两种状态磁力仪磁感应强度测量是相同的物理过程，转向差值也相等，即${B_{\mathrm{h}}}\left( \theta \right) = {B_{\mathrm{h}}}\left( { {\text{-}} \theta } \right)$ 。运用周期性将角度映射到$\left[ {0,2{\text π} } \right]$ 区间范围内，则
所以$ B\mathrm{_h}\left(\theta\right) $ 是关于π轴对称的函数。
但是，由于剩磁的影响，综合转向差$ {B_{\mathrm{H}}}\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 并不一定满足对称性关系，由式（3）可知，除了剩磁${{\boldsymbol{B}}_{\mathbf{r}}}$ 与待测磁感应强度${{\boldsymbol{B}}_{\mathbf{e}}}$ 相互正交或者${{\boldsymbol{B}}_{\mathbf{r}}} $ 与探头轴线重合外，其他情况均标定得到非对称的$ {B_{\mathrm{H}}}\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 。图4 所示的试验标定的一组$ {B_{\mathrm{H}}}\left( {\theta ,{{\boldsymbol{B}}_{\mathbf{r}}}} \right) $ 与$\theta $ 之间关系的散点图即为一个例子。一般情况下，转向差和剩磁很难从标定数据中区分。下面将从数值优化的角度建立剩磁补偿模型，估计剩磁的大小和方向，并从综合转向差数据中分离剩磁和物理转向差，从而建立统一的转向差补偿模型。
2 转向差的补偿方法
2.1 探头剩磁的估计
假定CPT原子磁力仪探头固有属性不变，即剩磁的大小和方向相对探头的位置均不变。建立如图5 所示的$OXYZ$ 笛卡尔坐标系，设定探头旋转面为$XOY$ 平面，待测磁感应强度${{\boldsymbol{B}}_{\mathbf{e}}}$ 与$Y$ 轴平行，${{\boldsymbol{B}}_{\mathbf{r}}}$ 为探头剩磁矢量，其大小和方向未知，设定${{\boldsymbol{B}}_{\mathbf{r}}}$ 与其在$XOY$ 平面上的投影${{\boldsymbol{B'}}_{\mathbf{r}}}$ 之间的夹角为γ，${{\boldsymbol{B'}}_{\mathbf{r}}}$ 方向与探头轴线之间的夹角为α，探头轴向与${{\boldsymbol{B}}_{\mathbf{e}}}$ 方向之间的夹角为β。$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 、α和γ为未知变量，β可以通过码盘读取或者根据磁通门磁力仪的方向数据确定。
由式（3）可知，剩磁引起的磁感应强度测量偏差为${{\boldsymbol{B}}_{\mathbf{r}}}$ 在${{\boldsymbol{B}}_{\mathbf{e}}}$ 方向上的投影，根据三余弦定理，该投影分量可表示为：
当探头沿着码盘平面转动时，采集输出数据，按式（3）计算可得出$ B\mathrm{_H}\left(\theta,\boldsymbol{B}_{\mathbf{r}}\right) $ 关于$\theta $ 的不对称散点图。$ B_{\mathrm{H}\left(i,j\right)}\left(i=1,2,\cdots,m;j=1,2,\cdots,n\right) $ 为计算得到的综合转向差数据。其中，i表示试验次数，共m次，即探头沿着旋转方向转动的圈数；j表示每次试验采集的转向差数据，共n个，即每圈采集的数据点数。去除剩磁影响的物理转向差数据表达式为:
建立优化的目标函数考虑以下三点：①剩磁引起的磁感应强度测量偏差$\Delta {B_r}$ 足够小，使得剩磁补偿后的散点曲线${B_{{\mathrm{h}}\left( {i,j} \right)}}$ 与标定得到的综合转向差数据散点曲线${B_{{\mathrm{H}}\left( {i,j} \right)}}$ 变化趋势保持一致；②剩磁补偿后，转向差散点曲线${B_{{\mathrm{h}}\left( {i,j} \right)}}$ 关于π轴对称位置数值差的平方和尽量小，以保证物理转向差的对称性；③剩磁补偿后，每次试验的转向差散点曲线与它们平均值之间的距离足够小，这反映的是剩磁补偿后不同组试验曲线的离散程度。因此，目标函数可以写成如下形式：
式中，$ N\mathrm{_{samples}} $ 为采集的数据总点数，主要起归一化的作用；$\left\lfloor {} \right\rfloor $ 为向下取整；$ mean(B_{\mathrm{h}\left(i,:\right)}) $ 为m次试验转向差数值的均值。
目标函数的约束条件主要是参数的取值范围，经弱磁化后，假定探头剩磁在10 nT以内，即$\left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| \leqslant 10$ ，α取值范围为[0,π]，γ取值范围为[0,π]。式（7）目标函数在约束范围内含有正余弦震荡，属于带约束的非凸优化问题，如果没有良好的初始条件，基于梯度或Hessian矩阵的最速下降法、高斯牛顿法、共轭梯度法等优化方法容易陷入局部最优解。这里采用粒子群算法用于目标函数优化。粒子群算法是一种随机搜索智能优化算法，通过个体间的协作与竞争，能够实现复杂空间全局最优解的搜索[13]。
合理设计种群粒子个数N，维数D，搜索次数G，学习因子${c_1}$ 、${c_2}$ ，惯性权重$w$ ，每个粒子的位置x[k] 和速度v[k]（k=1,···,N）以及边界即可运行算法。通过粒子群算法的全局寻优，理论上可得到$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 、α、γ的最优解或者它们之间的关系式，结合式（6）计算可得出多次试验的转向差$ B_{\mathrm{h}\left(i,j\right)} $ 与$\theta $ 相关的散点图。
2.2 利用三次样条插值建立转向差模型
利用三次样条插值的数值分析方法可以对给定的一组散点数据拟合出一条光滑的连续函数曲线[14]。剩磁补偿后，使用三次多项式逼近一段小区间内的转向差数据，并利用这些多项式的连接处衔接条件来保证整个曲线的光滑性。
考虑到环境干扰、电路噪声、试验误差等因素的影响，探头旋转在不同圈数、相同角度对应的转向差$ B_{\mathrm{h}\left(i,j\right)} $ 值不尽相同，单独选取一次试验数据插值的曲线也不能精确反映磁力仪转向差数值。将m次试验数据点关于π轴对称的试验点取平均值作为插值点纵坐标，插值拟合的曲线能够反映探头转向差的平均水平。插值点纵坐标的选取规则如下：
设$ \left\{ {{\theta _1},{\theta _2}, \cdots ,{\theta _n}} \right\} $ 、$ \left\{B_{\mathrm{ip}(1)},B_{\mathrm{ip}(2)},\cdots,B_{\mathrm{ip}(n)}\right\} $ 分别为n个插值点的横坐标和纵坐标的集合，在每个子区间$\left[ {{\theta _l},{\theta _{l + 1}}} \right]\left( {l = 1, \cdots ,n - 1} \right)$ 内，运用三次样条插值方法均能得到与之对应的4个三次样条曲线的系数$\left\{ {{a_l},{b_l},{c_l},{d_l}} \right\}$ 。如果无磁码盘或者磁通门磁力仪k时刻获取的探头轴向与待测磁感应强度之间的角度${\varphi _k}$ 在子区间$\left[ {{\theta _l},{\theta _{l + 1}}} \right]$ 内，则转向差补偿模型如下：
3 转向差的试验验证与分析
CPT原子磁力仪转向差补偿试验可分为三步：首先按照标定方法，获取若干组综合转向差与角度对应的离散点值；其次，用粒子群算法补偿剩磁引起的测量数据非对称部分；最后，使用三次样条插值生成转向差曲线。
实际综合转向差离散点数据获取是在中国计量科学研究院的稳定场发生装置下进行的，如图6 （a）所示。三对正交的线圈能够在中心位置产生稳定的磁感应强度，通过调节线圈电流设置待测磁感应强度大小为50000 nT，方向为垂直向上，无磁码盘固定在装置中心，且0°刻度与磁感应强度方向重合，CPT原子磁力仪探头可通过无磁机械结构绕码盘中心旋转。设置探头初始轴向角度为0°，即激光方向与磁感应强度矢量方向平行，在0°～360°范围内，顺时针每间隔15°记录一次CPT原子磁力仪稳定输出数据。重复５次试验，可以得到综合转向差值与角度之间的关系，如图6 （b）所示。从图中可以看出，５次试验散点的变化趋势基本一致，满足一定的周期性，且散点图均不是180°轴对称分布，这说明剩磁叠加在待测磁场方向上导致磁感应强度测量偏差，试验结果与上文理论分析一致。下面将用实测数据验证第３节所提出的转向差补偿方法。
设定粒子群算法的初始条件：待估计的未知参数有3个，故粒子维度设定为3。总粒子数设置为60，总搜索次数设置为100，学习因子、惯性权值为经验值，速度边界设置为[-1,1]。位置边界需要考虑参数的取值范围，对应目标函数的约束条件，$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 边界范围设定为[0,10]，α、γ边界范围均为[0,π]。试验的初始参数如表1 所示。
表 1 粒子群算法初始条件Table 1 Initial conditions of PSO |
| 变量 | 粒子数 | 粒子 维数 | 学习 因子 | 惯性 权值 | 速度 边界 | 位置 边界 | 迭代 次数 |
| 初值 | 60 | 3 | 1.2 | 0.8 | [−1,1] | [0,0,0;10,π,π] | 100 |
初始化种群以及算法中的随机值均采用MATLAB中的rand函数，待算法执行完成后可以得到剩磁$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 的大小以及α和γ的估计值。考虑到搜索算法的随机性，试验中，在算法最外层再加1000次循环，记录每次粒子群算法估计的参数最优值以及最优适应度值。绘制估计参数的概率密度直方图，如图7 （a）所示。α近似服从均值为90°的正态分布，方差为10−5量级，目标函数最优适应度值为0.0366，这两个参数均能得到稳定的估计值。$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 和γ并没有明显的确定性数值关系，而是近似服从指数分布，$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 比较小，在零附近的概率密度最大，而γ在90°附近的概率密度最大。考虑到最优适应度值稳定不变，且目标函数中$ \left|{\boldsymbol{B}}_{r}\right|·\mathrm{cos}{\text{γ}} $ 是整体出现，并且该值对应着$ {{\boldsymbol{B}}_{\mathbf{r}}} $ 在旋转面上的投影$ \left| {{{{\boldsymbol{B'}}}_{\mathbf{r}}}} \right| $ 。进一步得到算法执行1000次$ \left|{B}_{r}\right|·\mathrm{cos}{\text{γ}}$ 估计的概率密度直方图，如图7 （b）所示。该概率密度直方图也近似服从正态分布，因此，$ \left| {{{{\boldsymbol{B'}}}_{\mathbf{r}}}} \right| $ 值为0.344 nT。从上述分析可知，目标函数最优适应度值只与剩磁$ {{\boldsymbol{B}}_{\mathbf{r}}} $ 在旋转面上的投影分量$ \left| {{{{\boldsymbol{B'}}}_{\mathbf{r}}}} \right| $ 有关，而与$ \left| {{{\boldsymbol{B}}_{\mathbf{r}}}} \right| $ 和γ的单变量无关，且该投影方向垂直于探头轴线。
将估计的参数关系带入式（5）和式（6），得到５次试验剩磁补偿后的转向差离散点图，如图8 所示。相比于图6 （b）的标定原始数据，图8 中剩磁补偿后的转向差离散点图的对称度明显提高。按照式（8）计算插值点，运用三次样条插值方法可得光滑转向差曲线，任意角度的转向差值均能通过式（9）计算并补偿。
此外，转向差曲线的对称性也可以作为评估转向差补偿后剩余残差的依据，将５次试验综合转向差数据均值和剩磁补偿后的转向差数据均值分别与插值点作差，如图9 所示。该图反映了剩磁补偿前后转向差补偿残差的平均水平。从图９中可知，补偿残差关于（π，0）对称分布，这与式（8）插值点纵坐标计算规则有关。剩磁补偿前，转向差补偿后平均残差均在±0.175~±0.410 nT之间，并且0°~180°综合转向差值小于插值点纵坐标，180°~360°综合转向差值大于插值点纵坐标，这种不对称关系是由剩磁引起的。剩磁补偿后，大部分角度对应的平均残差在±0.182 nT之间，但在±165°达到0.321 nT，这主要是由于在180°附近探头光传输方向与待测磁感应强度的角度关系从$\theta $ 转化为$2{\text π} - \theta $ ，也可以解释为探头换向导致转向差曲线局部的曲率较大，微小角度的变化可能引起较大测量偏差。为了综合比较剩磁补偿前后对转向差补偿残差指标的影响，本节使用图9 中补偿残差的标准差来评估，计算可得标准差由剩磁补偿前的0.288 nT减小到剩磁补偿后的0.132 nT，标准差减小了54%。
以上试验验证了转向差补偿方法的有效性，使用粒子群算法从有剩磁干扰的非对称综合转向差散点数据中提取物理转向差数据，随后在此基础上使用三次样条插值获得转向差曲线，剩磁补偿后转向差补偿的平均残差大幅减小。在实际磁探测任务中，可以根据磁通门磁力仪提供的角度信息和预先标定以及剩磁补偿获取的转向差曲线补偿CPT原子磁力仪转向差。
4 结论
在非线性塞曼效应的影响下，CPT原子磁力仪在0°~360°范围内的转向差曲线理论上是轴对称的。虽然探头尽量避免使用磁性材料，但探头微弱剩磁依旧存在，并且物理上破坏了转向差的对称性，影响转向差的补偿精度。本文设计了转向差标定方法，建立了剩磁补偿模型，并用粒子群算法估计了剩磁参数。试验数据验证了模型和算法的有效性，探头剩磁在旋转面投影分量的大小为0.344 nT。剩磁补偿后，运用三次样条插值获得了对应的转向差曲线，并实现了转向差补偿。转向差补偿后的残差的标准差由剩磁补偿前的0.288 nT减小到0.132 nT，标准差减小了54%。