张笑楠，孙晓洁，寇 军，魏宗康，任 章

(1.北京航空航天大学自动化科学与电气工程学院，北京 100191；2.北京航天控制仪器研究所，北京 100039)

原子干涉磁力仪信号鉴频优化算法设计

张笑楠1，2，孙晓洁2，寇 军2，魏宗康2，任 章1

(1.北京航空航天大学自动化科学与电气工程学院，北京 100191；2.北京航天控制仪器研究所，北京 100039)

原子干涉磁力仪是通过确定CPT鉴频信号曲线峰值所对应的频率值来计算磁场强度值的.信号峰值定位锁定算法的搜索速度、准确度和稳定性直接影响磁力仪的分辨率和响应速度.常见的斜率算法只适用于小步长的搜索，定位速度慢，锁定状态的稳定性差.针对此问题，通过建立磁力仪鉴频信号曲线数学模型，设计了PI跟踪控制算法和快速锁定算法两种优化算法，并经仿真实验对比了性能.研究结果表明:3种算法中快速锁定算法不仅受初始搜索步长的影响小，而且峰值定位速度最快，稳定性好,对于原子干涉磁力仪系统的性能提高有很大的应用价值.

原子干涉磁力仪；CPT鉴频信号曲线；峰值定位锁定；优化算法

原子干涉磁力仪是基于相干布居囚禁(coherent population trapping，CPT)原理工作的.CPT是激光光场与原子相互作用时发生量子干涉的结果，是一种量子干涉现象.原子的磁子能级在磁场为零时都是简并的.当存在外磁场的情况下原子能级会发生塞曼分裂.当两个与原子作用的激光场的频率差与原子基态超精细子能级的劈裂间距精确匹配时，原子被束缚在这两个能级之间的无相互作用的暗态上.这个现象在激光强度变化上的表现为吸收光谱中出现非常窄的透射信号.通过对该光谱信号对应频率的测量和计算可以得到外磁场的强度值[1-3].

传统的光谱信号鉴频方法是通过求取信号的一次微分来实现的，也就是将信号的幅值信息转化成斜率信息进行鉴别.由于单个信号峰的斜率是由小到大连续变化的，峰值点对应的斜率为零值，所以通过数值计算定位零值所对应频率就可以达到鉴频的目的.目前国内外CPT测磁实验中最常用的方案是使用锁相放大器和滤波电路模块获得光功率信号的一次微分曲线[4].该方案如果使用小型电路和程序来替代精密仪器嵌入小型样机中，就会受计算速度的限制，并且对数据存储量的要求也较大，稳定性不高.所以本文从提高曲线鉴频计算速度、数据量和稳定性方面来进行优化算法的设计.

1 原子干涉磁力仪信号数学模型

由于原子干涉磁力仪是通过对原子干涉信号进行精确分析来计算相应的磁场强度的，为了深入分析信号的处理方法，需要首先建立鉴频信号的数学模型.原子干涉过程中原子系统的演化可以用原子密度矩阵的形式表达[5-6].在此基础上将光量子传感器感应到的与原子作用的激光光强信号分为3个组分，分别是激光器射入原子气室的入射激光功率Pin，原子气室吸收的光功率Pa和CPT信号.传感器内的光电探测器感应到光功率信号表示为Pout，建立它们之间的关系模型如下.

Pout=Pin-Pa-Pcpt，

(1)

式中，

(2)

(3)

(4)

(5)

式(2)中的w1和ωa分别为激光器实际输出频率和设定频率；式(3)中的Г为原子激发态到基态退相干；λ为激光器输出激光光波长；n为原子气室内充原子数；m为原子质量；Tcell为原子气室温度；kb为玻尔兹曼常量；式(4)中的c为光速；式(5)中的H为CPT信号幅值系数；wset为微波源设定频率；dw为微波源设定频率与实际输出频率的差值；γ为原子旋磁比；B为外磁场强度；υ为CPT信号线宽.为了便于计算，将以上所建立的信号数学模型在峰值附近的小范围内进行简化，得到数学表达式如下：

(6)

其中，ω为磁场对应的谱线频率；ω0为碱金属原子基态两超精细能级的半频差，二者的频率差dω=ω-ω0.根据该数学模型编写Matlab程序来对原始鉴频信号数据进行拟合，结果如图1所示，可以看出该模型在峰值附近的数据拟合效果良好，故以下使用该简化模型来进行峰值寻找锁定优化算法的设计和仿真.

图1 采集到的系统鉴频信号和数据拟合效果Fig.1 Acquired system frequency discrimination signal and data fitting result

2 峰值定位锁定方案

磁力仪通常采用斜率法作为曲线寻峰算法.斜率法寻峰是根据CPT信号峰模型曲线的斜率变化规律来确定峰值位置的[7-9].图2显示的是CPT模型中选取频率差dω为-100～100时计算得到的相应斜率值.

图2 CPT模型中频率差值所对应的斜率值Fig.2 The slope value of the frequency difference in the CPT model

CPT信号峰值附近小范围的频率差值所对应的斜率值是逐渐增大的，当频率差值正好为0时，也就是到达峰值处时得到的斜率值等于0，且峰值左侧和右侧的频率值符号相反.根据信号曲线斜率的这种特点，斜率算法峰值锁定的实现步骤如下：

(1) 设定曲线上的起始点dω0和固定步长d_ω；

(2) 确定寻找方向，比较前后两点的斜率值D(H0)和D(H1)大小，如果D(H0)

(3) 迭代计算，按固定步长连续计算频率差值并进行方向判断，如果当前计算的斜率值与初始斜率值同号，则继续正向寻找，否则反向寻找.

斜率寻峰法的优势在于程序结构简单，可以根据实际情况设置合适的步长并快速定位峰值.但是，在到达目标之后如果继续迭代计算，则会在峰值处以步长值保持振荡，搜索步长越大则振荡越剧烈.

3 鉴频信号算法设计

本文设计了两种新的峰值定位锁定算法，即PI跟踪控制算法和快速锁定算法.

3.1PI跟踪控制法

PI跟踪控制法是基于PI控制器对误差进行实时跟踪消除的一种算法.由式(6)可知，在理想情况下，也就是当频差为0时，CPT曲线的光功率值达到最低点.在峰值对应频率附近进行小范围的频率扫描，当扫描范围的中心频率值与ω0的差值dw刚好为0时，从扫描曲线中可以看到该扫描中心频率值对应着CPT曲线的极值点，并且得到的功率谱线完全对称.当dw不为0时，扫描范围中心与峰值频率发生偏移，功率谱线不再对称.例如，取扫描范围为200 Hz，扫描范围中心频率偏移值为20，得到的光功率谱线如图3所示，纵坐标表示光功率.所以根据扫描曲线的对称性可以实现对频率的调整和锁定.

(a) 中心偏移值为0 Hz时的CPT信号曲线 (b) 中心偏移值为20 Hz时的CPT信号曲线图3 扫描范围为200 Hz时得到的CPT信号曲线Fig 3 The CPT signal curve obtained when the scanning range is 200 Hz(a) CPT signal curve when the center offset is 0 Hz; (b) CPT signal curve when the center offset is 20 Hz

PI跟踪控制法的具体实现步骤如下：

(1)定义算法误差量.定义一个误差量e来衡量扫描曲线的对称性.把扫描频率范围以中心为界分为两组频率值，分别记为数组f1和f2，两组频率对应的功率值分别为数组HL和HR，对两个功率值数组内的数据求和，分别得到ER1=sum(HL)和ER2=sum(HR).这两组功率和值的差值可以作为寻找峰值的误差量，表示为e=ER1-ER2.误差量e和dw的关系为线性变化.

(2)设定PI控制器参数.为了使频率偏移值保持为0，需要消除扫描模型计算得到的误差量.本系统使用PI控制器来消除误差，它的表达式为[10]：

其中，u为系统输出量；e为随时间变化的误差量；Kp为比例系数；Ki为积分系数.将上式离散化成差分方程后，第k拍的输出表示为：

(3)频率锁定.反复迭代(1)和(2)的计算步骤，直到误差逐渐被减小到0时，进入峰值频率锁定状态.

3.2 快速锁定算法

图4 快速锁定算法加速确定搜寻点方法示意图Fig.4 Search point determination acceleration method diagram of fast locking algorithm

快速锁定算法的设计是基于原子鉴频信号峰谱线模型，如图4所示，以峰值点为虚拟横坐标轴，通过两个搜寻点dω0和dω1来推断其所在直线与该虚拟轴的交点，从而加速寻找下一个最接近峰值的搜寻点，通过调整频率偏移值dω=dω1-dω0来实现峰值定位.

快速算法具体的实现步骤如下.

(1) 由式(7)得到其一阶近似模型

当Pin-Pa-H为已知量时，取系统的观测量为

(2) 给定初始值dw0,求P(dw0).

(3) 给定增量.设置一个正向增量d_w，此时，有dw1=dw0+d_w，求出相应的P(dw1).

(4) 确定方向.如果P(dw0)>P(dw1)，表明当前频率差值位于锁定频率的右侧，与希望的锁频方向相反，则需重新选取dw1=dw0-d_w，并求取P(dw1).

(5) 构造直线.由两点(dw0,P(dw0))和(dw1,P(dw1))计算出所在直线的斜率值：

(6) 求取下一个点.设在下一个频率点，在直线方程中有P(dwm)=0，则该点的表达式为：

此点对应的观测量为P(dwm)：

(7) 数据更新.取dw1=dwm，P(dw1)=P(dwm)和dw0=dw1,P(dw0)=P(dw1)，代入步骤(6)，进入持续迭代过程，当搜索点到达峰值位置时，锁定在峰值点处.

(8) 在Pin-Pa为未知量的情况下，仍然使用观测量P，此时需根据观测量不断调整H0值.执行步骤(2)到(7).

(9) 求取下一个点.设在下一个频率点，在直线方程中有,则该点的表达式为

此时判断dωm和dω1的符号，如果二者同号，则此点对应的观测量为

否则，对Pp进行更新，有

然后更新Pp和P(dωm).继续执行步骤(7)，然后重复执行步骤(9).

为了更直观地进行分析，以初始值为-100，步长为3的设定情况为例来观察算法迭代和稳定的表现，仿真结果如图5所示，左图的纵坐标表示光功率.

图5 锁频过程和收敛过程Fig.5 The process of locking process and convergence

图5展示了利用快速锁定算法实现的锁频过程和收敛过程，左图的纵坐标表示光功率.可以看出在程序运行能够实现目标加速搜索，在快速定位峰值后稳定状态良好.

4 算法仿真结果分析

图6 3种算法的仿真结果(dω=-100 Hz，d_ω=3 Hz)Fig.6 Simulation results of the three algorithms

为了观察两种优化算法与传统斜率算法的计算速度和收敛稳定性，分别对3种算法编写程序进行仿真并对结果进行比较.3种算法的频率差值dω均设置为-100 Hz.仿真结果如图6所示.

其中，斜率法的步长选择为d_ω=3 Hz，仿真结果在图中用三角形数据点表示.当迭代计算到33步时，dω=0 Hz；PI跟踪控制法以dω=-100 Hz为扫描范围的中心频率值，扫描范围为dω-100 Hz～dω+100 Hz，同时设定控制器比例放大系数kp=0.01，积分放大系数ki=0.000 01，仿真结果在图中用圆圈数据点表示，当迭代步数为41时，dω被调整为0，之后系统锁定在峰值处；快速锁定算法的仿真结果在图中用星型数据点表示，搜索点在开始时快速逼近峰值点，随后减慢速度，当迭代到第11步时找到峰值，之后持续迭代计算并锁定在峰值位置.

从仿真结果可以看出，斜率算法始终按照固定步长更新数据，步长设置越大，峰值定位所需步数越少，但是在到达峰值处之后在dw=0附近保持振荡，振幅等于设置的步长值，这也意味着该算法不能使用大步长搜索；PI跟踪控制法搜索峰值所需的迭代次数最多，但是与斜率算法相比，锁定峰值点的稳定性更好；快速锁定算法在进行搜索点初始决策时，能快速逼近峰值点，所用的迭代次数最少.另一方面，该算法越接近峰值点时计算步长越小，可以实现峰值频率的准确定位，并保持很好的稳定状态.

综上，快速算法的综合性能与其他两种算法相比更为出色.一方面它的峰值定位速度快，在初始步长要求设置较小的情况下具有绝对的优势，另一方面，该算法在峰值频率锁定之后的稳定性高，不受初始步长值大小的影响.根据原子干涉磁力仪系统鉴频信号的控制需求，最终选用快速锁定算法作为CPT信号曲线锁定算法.

[1] STAHLER M,KNAPPE S,AFFOLDERBACH C,etal.Picotesla magnetometry with coherent dark state[J].Europhys Lett,2007,54(3):323-328.

[2] NAGEL A,GRAF L,NAUMOV A,etal.Experimental realization of coherent dark-state magnetometers[J].Europhys Lett,1998,44(1):31-36.

[3] 刘国宾,孙献平,顾思洪,等.高灵敏度原子磁力计研究进展[J].物理,2012,41(12):803-810.

[4] 黄宇翔.CPT原子磁力仪的初步研究[D].杭州：浙江大学,2013.

[5] 王 丰,刘 强,曾宪金.Cs原子磁力仪共振谱线宽度的研究[J].光电子·激光,2010,21(6):846-847.

[6] 满文庆，杨世琪，钟旭滨，等.以原子谱线作参考的半导体激光器的频率锁定[J].激光技术,1998,22(1):9-10.

[7] 邹鹏毅,罗深荣,顾建松，等.两型光泵磁力仪比对试验结果及分析[J].声学与电子工程,2008,41(2):35-45.

[8] 李 芳.磁通门磁力仪的设计与研究[D].长春：吉林大学,2014.

[9] 柯超凡，彭 翔，徐盛晨，等.基于虚拟仪器的光泵磁力仪测控系统设计[J].红外与激光工程,2013,42(12):3316-3319.

[10] 任 昊，王学锋，杨学礼，等.输入光功率对探测器组件带宽的影响分析[J].导航与控制,2012,11(1):63-67.

(编辑 CXM)

Frequency Discrimination Optimization Algorithm Design for Atomic Interference Magnetometer

ZHANGXiao-nan1，2,SUNXiao-jie2,KOUJun2,WEIZong-kang2,RENZhang1*

(1.School of Automation Science and Electrical Engineering,Beihang University,Beijing 100191,China;2.Beijing Institute of Aerospace Control Instrument,Beijing 100039,China)

Atomic interference magnetometer is able to calculate the value of the magnetic field strength by collecting the CPT frequency discrimination signal curve peak which corresponding to the frequency value though the light quantum sensor.The resolution and response speed of magnetometer are directly affected by the search speed,accuracy and stability of signal peak positioning and locking algorithm.Slope algorithm is widely used in magnetometer system and only suitable for small step length search with slow location speed and low stability of the locked state.In order to improve the algorithm,two kinds of optimization algorithms are designed based on the mathematical model of magnetometer signal,which are PI tracking control algorithm and fast locking algorithm.Simulation experiments are also conducted to compare the performance of the three algorithms.The results show that the fast locking algorithm is nearly not affected by the initial search step size,which also has the fastest peak position speed,and good locking stability.It is of great application potential for the performance improvement of atomic interference magnetometer system.

atomic interference magnetometer; CPT frequency discrimination signal; signal peak positioning and locking; optimization algorithm

10.7612/j.issn.1000-2537.2017.01.011

2016-09-22

国家自然科学基金资助项目(11204011)

* 通讯作者,E-mail：meiye100@hotmail.com

U652

A

1000-2537(2017)01-0071-06