彭向凯 吉经纬 李琳 任伟 项静峰 刘亢亢 程鹤楠 张镇 屈求智 李唐 刘亮† 吕德胜‡

1)(中国科学院上海光学精密机械研究所,量子光学重点实验室,上海 201800)

2)(中国科学院大学材料与光电子研究中心,北京 100049)

1 引 言

由于原子在外磁场作用下的塞曼效应导致原子能级分裂与频移,使原子钟、原子干涉仪以及冷原子重力仪等以原子作为传感器的精密设备都会受环境磁场影响[1−5],在地面实验室环境中,主要是地磁场的影响.解决此问题的通常做法是对待测原子样品做多层磁屏蔽以降低环境磁场对屏蔽内部原子的影响,磁屏蔽材料一般选用磁导率较高的坡莫合金铁磁材料[6,7].但是对于绕地球飞行的星上设备来说,比如星载原子钟与空间冷原子钟,感受到的地磁场和地面实验室的很大不同之处在于卫星会经历一个周期变化的磁场[8−10].例如对于搭载在“天宫二号”空间实验室的空间冷原子钟来说,卫星在近地约400 km的斜45°轨道运行,经历磁场变化受卫星绕地周期与地球自转两方面的影响,磁场随时间变化近似于被20 d周期调幅的90 min周期正弦曲线[11].在这样的环境磁场影响下,星上原子钟磁屏蔽内的原子样品感受到的磁场不仅有被屏蔽衰减后的剩余磁场,更严重的是有铁磁材料在变化外场下被磁化产生的感应交变磁场.

如果能够精确测量磁屏蔽内感应磁场,可以通过洛伦兹线圈产生动态补偿磁场来有效消除变化外场的影响,由于磁场传感器自身磁性对原子样品会产生干扰,在靠近原子样品的屏蔽内放置磁场传感器的方法并不适用,需要通过探测环境磁场来精确预测屏蔽内感应磁场.磁屏蔽内感应磁场随外磁场呈磁滞曲线变化,可以通过线性拟合或多项式拟合的方法近似预测[8,9],更精确的方法是使用Jiles-Atherton(J-A)磁滞模型[12−14],该模型是目前最为常用的描述铁磁性材料磁滞特性的经典数学模型,是一种具有物理基础的现象学模型,能够真实地描述磁屏蔽内感应磁场与外磁场的非线性关系,通过求解J-A磁滞模型方程便能够得到较为准确的磁滞回线,从而根据探测到的环境磁场精确预测磁屏蔽内感应磁场

在J-A磁滞模型方程中,代表磁性材料特性的一组J-A参数是整个问题的基础,目前还没有一种快速、便捷、精确的方法对这一问题进行完整的解决,通常是在经验值的基础上根据特性磁场环境实验测试进行拟合.机器学习技术,尤其是依赖人工神经网络的深度学习,有极其强大的拟合能力[15−17].同时使用在线学习的技术,可以自动寻找到最优化的参数值[18−20].这个过程往往远快于人工调节,并且结果往往优于人工调节.机器学习在线优化技术已经应用于相当多的量子实验中[21−25],而深度神经网络也已经应用于冷原子实验的磁光阱(magnetooptical trap,MOT)参数的优化[26].之前也有研究利用神经网络结合遗传算法通过离线学习对硅钢片J-A参数进行优化拟合,其结果优于利用解析模型的拟合[27,28].神经网络从数学模型角度可以看作对复杂函数的拟合,多参数拟合对人工计算是比较困难的工作,对于现在电脑来说却是非常简单的工作,参数越多这种差异越明显.使用神经网络,得益于电脑的运算速度,即便数百个参数的拟合仍然可以快速准确地完成,而对于人工来说几乎是不可能完成的.本文使用神经网络模型在线学习对JA参数进行最优化拟合,可以用更少的时间找到最优参数,并且应用于大规模参数时具有更好的鲁棒性[26].根据拟合的参数,在模拟卫星载荷感受地磁变化的情况下对磁屏蔽内剩余磁场进行补偿,得到了优于手动参数拟合的结果.

2 将J-A模型应用于磁屏蔽

为了将J-A模型应用于磁屏蔽,我们将屏蔽内的剩余磁场分为两部分: 被磁屏蔽衰减后的环境磁场BSE以及磁屏蔽本身被磁化后产生的感应磁场BM[29],因此有

第一项以磁屏蔽系数SE随环境磁场线性变化;第二项与磁屏蔽的磁化强度成比例,可表示为BM=C×M,C为比例系数,有和真空磁导率相同的量纲,单位为nT × m/A,在不同的位置有不同的取值,可以通过有限元仿真得到各个位置的取值.屏蔽内剩余磁场由于BM的存在产生了磁滞回线的形状.

使用J-A模型求磁屏蔽被环境磁场磁化后的磁化强度M,J-A模型由下式表示[30]:

其中He=µ0(H+αM)是J-A模型中的有效场,µ0是真空磁导率,H为环境磁场;Ms,k,a,a,c为J-A参数,需要通过实测结果来确定.当 dH<0 时d=1,当 dH⩾0 时d=–1.

为了简化计算以及将来更方便地用于磁场补偿系统,我们对模型进行了合理的简化.首先,磁畴耦合参数a在软磁材料中影响很小[31],我们将软磁材料a经典值4 × 10–4与0值分别代入J-A模型计算屏蔽内剩余磁场,发现磁屏蔽内剩余磁场计算值差别小于0.5 nT,并且a取0时,积分方程可以解析求解展开为多项式形式,这样可以在结果影响很小的情况下极大地提高计算效率,便于软件编程,这对于实现磁场的实时补偿是很必要的.另外,当He/a<0.1 时,可以近似使用在空间磁场环境下,这个条件总是满足的.基于上述简化,公式可以改写为

其中当x> 0时 [x]+=x,x≤ 0时 [x]+=0.

此外,我们对环境磁场采集是离散值而非连续值,因此为了方便计算可以将公式中M(t)改写为M(Hn),即用当前环境磁场状态下的磁化强度表示当前时刻的磁化强度,最终可以得到整个剩余磁场的表达式

在选择最优参数值的过程中,首先根据经验值对参数进行一些简化和范围的限定,限定范围首先可以确保得到的参数值具有合理的物理含义,并且可以更加迅速准确地拟合出参数值,同时又可以防止得到的参数值陷于局域最优值而非全局最优值.饱和磁化强度Ms可以根据坡莫合金的参数表得知在5.4× 105附近,我们将范围限定于5 × 105—6 × 105.扎钉参数k近似等于矫顽力,可以根据参数表将范围限定于1—50.可逆参数c的范围根据定义限定于0—1.Langevin梯度参数a根据实测磁滞图形可以限定于1.5 × 104—2.5 × 104.

在本次实验中,测试对象是单层磁屏蔽,对于更广泛应用的多层磁屏蔽来说,只需要将多层磁屏蔽当作一个系统,对该系统利用上述J-A方程分析,求解J-A参数、衰减系数SE和比例系数C,然后对屏蔽内关心的位置展开计算分析,就可以准确计算屏蔽内剩余磁场.

3 实验测试装置

为了模拟近地轨道的磁场环境,我们搭建了一个准亥姆霍兹线圈,通过控制线圈电流来控制中心轴向磁场.该装置由三组八边形线圈构成(便于使用铝型材搭建),线圈直径约为1.3 m,整体高度约为1.5 m.我们将待测屏蔽筒置于线圈内部处,屏蔽筒直径为30 cm,高度为80 cm,并在屏蔽筒内部中心位置放置一个磁通门用于采集屏蔽内磁场,在线圈轴线远离屏蔽筒处放置另一个磁通门用于采集环境磁场.图1(a)是整个测试系统的示意图,F1和F2分别为探测环境磁场和屏蔽内部磁场的磁通门.图1(b)是用于模拟近地轨道磁场的环境磁场变化图.实验过程中,首先让线圈电流产生正弦变化的环境磁场,根据得到的磁滞回线作为神经网络训练与反馈的值来获得屏蔽筒J-A参数,并将这组参数用于预测近地轨道环境磁场条件下的屏蔽内磁场,将预测值和实测结果进行对比,验证神经网络算法预测J-A参数的准确性.

图1 (a)实验装置示意图,同时也是有限元仿真时所用的模型图;(b)对模型进行测试所用的近地轨道环境磁场Fig.1.(a)The schematic diagram of the experimental device and also the model diagram used in our finite element simulation;(b)low Earth orbit environmental magnetic field used to test the model.

在利用神经网络对J-A参数优化之前,除了预先设定参数范围外,预测屏蔽内磁场还需要知道屏蔽筒内各处衰减系数SE和磁化强度与磁场之间的比例系数C,我们使用有限元分析方法(finite element method,FEM)[31]仿真屏蔽筒内空间各处的衰减系数SE和比例系数C优化前的初始值.图2展示了屏蔽内部一个纵向截面上有限元计算所得的SE和C.我们探测所用磁通门置于屏蔽最中心处,在此处的衰减系数约为31,比例系数约为–1.8.使用其他方法获得J-A参数时,这两项参数一般都仅仅使用仿真值,虽然有限元模型可能存在的各种误差,但是往往因多参数之间关系复杂不对这两个参数进行调节.而使用神经网络自动调节参数时,由于神经网络具有强大的仿真计算能力,我们将这两个参数也作为代求参数.对SE进行有限元仿真时,选择了厂家给出的相对磁导率,实际上磁导率可能无法严格等于厂家给出的参数,并且屏蔽效果可能由于形变等原因与仿真结果有一定偏差,因此我们根据仿真结果31,将范围限定于20—40.对C进行有限元仿真时,我们将磁屏蔽看作永磁体,设置磁化强度,能够得到空间中不同位置的感应磁场,从而计算比例系数.同样,由于形变、位置偏差等一些原因,实际结果可能也会和仿真结果有偏差,因此根据有限元仿真结果–1.8,因此我们将范围限定于–1.5—–2.5.

图2 (a)数据所在磁屏蔽筒截面示意图;(b)有限元仿真所得到屏蔽内一个截面上的屏蔽效率分布图;(c)有限元仿真得到的屏蔽内一个截面上的磁化强度比例系数图Fig.2.(a)A schematic view of the section of the data;(b)the shielding efficiency distribution on a section inside the shield calculated by FEM software;(c)the magnetization intensity ratio coefficient on a section of the shield calculated by FEM software.

4 神经网络模型的建立

本文需要优化的是J-A模型中的5个参数Ms,k,a,a,c和衰减系数SE以及比例系数C共7个参数.我们建立如图3的参数优化神经网络.

在图3中,神经网络作为一个复杂的函数,我们给定输入参数值,图中x1,x2,···,x7对应J-A模型中的5个参数Ms,k,a,a,c和衰减系数SE以及比例系数C,每一组参数可通过神经网络对应一个输出值y,同样的参数值,在给定环境磁场的情况下,代入(4)和(5)式,可以算出相应的剩余磁场,并和实测的剩余磁场相减后求标准差y0,我们使用神经网络这样一个通用的函数形式来拟合这个参数值到标准差的复杂过程,一旦完成拟合,使用scipy库[32]可以方便求出令输出值即标准差最小的参数值.此外我们使用了在线学习的方法,视情况将预测的参数值重新代到神经网络中进行训练,这样可以充分利用每组参数,用尽可能少的轮次拟合出最准确的神经网络.

神经网络模型建立过程中,网络层数与每层的神经元数是首先需要考虑的问题[15].神经网络层数以及神经元数量会影响神经网络的性能,如果层数较少,或者神经元数量较少,可能会导致无法拟合出复杂的对应关系(欠拟合),如果层数较多,或者神经元数量较多,有可能会将原本简单的关系拟合复杂(过拟合).神经元数量依惯例通常选用2的次方个,我们尝试过64个神经元,以及3层或4层网络,结果显示出明显的过拟合,即预测最优参数的标准差起伏很大且不理想,如图4(b)所示.我们同时也尝试过16个神经元,2层网络,结果显示出明显的欠拟合(所需轮次较多,且调参结果没有32个神经元结果好).综合上述实验,我们认为选用2层,每层32个神经元,对于我们这个只有7个馈入节点的神经网络来说是比较合适的,如图4(a)所示.随着参数的增加,我们需要增加神经元以及神经网络层数,具体所需数量还要根据实验来确定.

首先随机选择12组参数并求得相应的标准差对神经网络进行训练,梯度下降方法使用Adam方法,每一层的激活函数选用Relu函数,训练轮次为500次,完全在线学习,批量选为1.根据训练得到的神经网络,使用L-BFGS-B方法预测当前最优参数值[33],为了防止参数陷于可能的局部最优值而非全局最优值,我们对每一组预测最优值增加一个很小的随机偏置量,并将求得的标准差代入神经网络继续进行100轮次的训练.此外,若连续5次预测得到的标准差均大于之前得到的最优值的标准差,将重新完全随机选择下一组参数值,以求跳出局部最优值.训练过程持续到得到一个优于预期最低标准差的值或者到达最大循环次数或者连续20次预测值没有优于之前的最优值.流程图如图5所示,图中只显示了一个结束条件.

图3 神经网络调参原理示意图 使用2个隐藏层的全连接神经网络,每层32个神经元,一旦神经网络训练完成,就能调用scipy库预测最优参数Fig.3.Principle of optimizing with neural network.We use 2 hidden layers of full connected neural network,32 neurons per layer.Once the neural network training is completed,we can use scipy to predict the optimal parameters.

由于注入值较少,选用了一个两层的的神经网络,因此网络训练时间很短,预测最优参数也很快,每组参数求标准差每组需要约2 s,求得一组理想参数大约仅需要3 min.我们使用keras[34]建立神经网络,使用scipy调用L-BFGS-B方法进行最优参数的预测.相比神经网络对参数拟合,我们利用人工仅对J-A模型中的5个参数进行优化拟合,一般情况也需要几个小时.

图4 (a)当选用2层网络各32个神经元时,预测标准差收敛性好且最小值仅为4.9;(b)当选用3层网络各64个神经元时,由于过拟合预测标准差振荡且最小值为27Fig.4.(a)When 32 neurons of the 2 layers network are selected,the convergence of the prediction standard deviation is good and the minimum value is only 4.9;(b)when 64 neurons of the 3 layers network are selected,the convergence of the prediction standard deviation is poor because of over-fit and the minimum value is 27.

5 参数优化结果与验证

使用神经网络调参的结果如图6所示,首先用一个正弦变化的环境磁场进行调参.图6(a)横坐标为总循环次数,纵坐标为我们设定的标准差,蓝色点为初始随机地用于训练神经元的12组参数,可见在第20轮左右,即初始训练后预测的第8组参数就已经开始明显收敛于最优值,最终得到的JA参数为Ms=540074,c=0.58,a=19817,k=9.50,SE=36.37,C=－2.37.标准差值为4.90,在该组参数下计算得到的磁滞回线和实测值如图6(b)所示,红线为预测值黑线为实测值.

相应地,当我们手动调节参数时,根据参数表确定Ms=540000,并根据有限元仿真结果确定SE=31,C=–1.8.在此基础上调整k,c和a,根据经验观察,a和幅值关联较大,k与c和磁滞程度关联较大,经历长时间反复调节之后,得到的最优的标准差值为7.09,在此时实测和预测的磁滞回线的形状仍然高度重合,很难进一步调节,然而手工调参和神经网络调参的标准差相差了接近一倍,意味着仍有很大的可调空间.显然使用神经网络自动调参,节省了工作量,节约了大量时间,同时得到了更加准确的参数.另外,最优参数得到的SE和C的取值与有限元结果的差异,提示了可能存在于有限元仿真过程中的误差,因素可能是材料磁导率的理论值和实际值的偏差、模型和实际屏蔽结构误差等.表1显示了手动和自动方法获得参数的对比.

图5 自动调参的流程图Fig.5.The flow chart of optimizing magnetic shielding characteristic parameters.

图6 (a)一个典型的求参过程图,预测参数值的标准差随实验轮次逐渐降低;(b)相应参数计算得到的磁滞回线与实测回线的比较Fig.6.(a)A typical continuation process graph,the standard deviation of the predicted parameter values is gradually reduced with the experimental round;(b)comparison of hysteresis loop and measured loop calculated by corresponding parameters.

另外,初始随机参数的不同可能导致不同的预测速度,但通常都能够收敛到理想的参数点.图7展示了另两组初始随机参数进行预测的标准差和循环轮次图,同样蓝色点为初始随机的用于训练神经元的12组参数,可以看出由于初始随机值不同,在本次训练中更快达到了标准差小于5的期望值.

表1 手动调参和自动调参得到的参数值Table 1.Parameter values obtained by manual tuning and automatic tuning.

接着将该组参数应用于接近实际卫星在轨经历的环境磁场,得到的结果如图8所示.图8(a)和图8(b)分别是在环境磁场随时间变化情况下,利用使用神经网络调参后得到的参数来计算屏蔽内磁场和实测屏蔽内磁场的对比图,横坐标分别为时间和环境磁场; 图8(c)表示计算值和实测值的差值随时间变化曲线; 图8(d)表示计算值和实测值的差值随环境磁场变化曲线.图中标准差越小说明计算越准确.其中黑线是仅考虑衰减项BSE不考虑磁化项BM的计算结果差值,另外两组则分别用手动得到的参数和神经网络得到的参数来计算,其中手动找到的参数标准差为9.33,神经网络调参标准差为7.01.可见使用J-A模型计算屏蔽内参数是非常必要的,仅考虑衰减项而忽略磁滞项将会导致很大误差,另外利用神经网络获取J-A模型的参数比手动参数更加快速准确,为我们对空间钟磁屏蔽补偿参数的选择提供了有效的帮助.

图7 另外两组不同随机初始训练参数下,标准差随实验轮次的收敛图Fig.7.Convergence graph of standard deviation with experimental rounds under two different random initial training parameters.

图8 在近地轨道磁场环境下实测磁场与计算磁场的对比图(a)横坐标为时间;(b)横坐标为环境磁场;(c)使用手动调节参数,自动求得的参数,以及仅考虑衰减项不考虑磁化项,计算磁场与实测磁场的差值,横坐标为时间;(d)横坐标为环境磁场Fig.8.A comparison of the measured magnetic field with the calculated magnetic field in a near-Earth orbit magnetic field environment;(a)The x axis is time;(b)the x-axis is the ambient magnetic field;(c)use manual adjustment parameters,automatically obtained parameters,and simple BSE to calculate the difference between the magnetic field and the measured magnetic field,the abscissa is time;(d)the x-axis is the ambient magnetic field.

6 结 论

基于描述铁磁材料磁滞特性的J-A模型,我们对坡莫合金磁屏蔽装置在变化外磁场情况下屏蔽筒内磁场进行预测.结合有限元模型对磁屏蔽特性的初步分析,使用神经网络在线学习的方法以获得磁屏蔽的J-A参数与屏蔽特性的最优值,并模拟通过低轨卫星经历的地磁场情况下验证拟合参数的准确性.拟合结果与实验测试表明,使用简单的两层神经网络就可以高效、快速地获得超过手工调节精度的J-A参数.一方面通过该方法对磁屏蔽特性的精确测定,可以应用于在轨运行的高精度磁敏感卫星载荷,通过在轨磁场补偿的措施,获得更好的磁屏蔽效果; 另一方面验证了神经网络应用于空间冷原子钟相关系统参数调整的可行性.该方法从理论上扩展了目前对于多参数拟合优化的工具,对于实际应用中通常会遇到的多参数物理实验寻找系统最优解会有很大帮助.例如可以应用于类似复杂系统的调参工作,对冷原子实验中与多个光功率、频率、时序相关的冷原子数目和冷原子温度进行优化[35,36],相关参数往往有数百个之多,受限于在轨运行时资源,人工调节几乎没有可能找到最优参数.我们希望神经网络的应用能够帮助优化空间冷原子钟的性能,并且实际最优参数与理论最优参数的差别有可能帮助我们进一步反思微重力条件下原子冷却的物理机制.