相位檢測誤差對核磁共振陀螺輸出頻率的影響

江奇淵，羅 暉，楊開勇，趙洪常，汪之國，夏 濤，張 燚

(國防科技大學前沿交叉學科學院，長沙 410073)

0 引言

核磁共振陀螺通過測量原子系綜核自旋的Larmor進動頻率在慣性參考系下隨載體轉動的變化量來實現轉動測量，而原子系綜核自旋進動的測量通常是通過內嵌的磁力儀測量其磁化強度矢量來實現的。原子氣室作為貯存其工作物質的容器，通常包含堿金屬原子(Rb，Cs)、惰性氣體原子(Xe，He)及緩沖氣體[7-8]。堿金屬原子不僅作為其內嵌的磁力儀的工作物質，同時通過光泵浦和自旋交換作用實現惰性氣體原子核自旋的超極化[9]。由于核自旋產生的磁場很弱，單核自旋通常無法通過堿金屬磁力儀實現探測，因此，需要對其原子系綜進行探測。原子系綜內各原子核自旋由于相位是隨機的，因此，其橫向宏觀磁化強度矢量為零，無法探測到其Larmor進動頻率，需要在橫向施加一個與惰性氣體原子Larmor進動頻率相等的激勵射頻場來統一各原子核自旋的相位，實現其原子系綜宏觀的進動效應。最終通過記錄這個實時跟蹤惰性氣體原子Larmor進動頻率的激勵射頻場，實現轉動測量。

因此，如何精確、實時地跟蹤探測到的核自旋進動磁場，并施加反饋激勵磁場，就成了核磁共振陀螺的關鍵技術之一。通常用于其信號跟蹤的方法有反正切法、快速傅里葉變換法(Fast Fourier Transform，FFT)[10]、自激振蕩以及鎖相環方案[11]等，這些頻率跟蹤方案大多數都是通過探測信號的相位量進而得到頻率值的。理想工作的頻率跟蹤方案能夠實時跟蹤探測信號的相位及幅度，并對時變的相位作微分處理得到信號的頻率，但是由于相位延遲等因素導致的額外相位所引入的頻率，無法與Larmor進動頻率區分開，而這樣的頻率差值就可能會引入額外的測量誤差。

核磁共振陀螺原子系綜的輸出頻率測量精度決定了其作為陀螺儀的性能，因此，建立起它的頻率誤差方程，找出引起頻率誤差的因素，具有十分重要的意義。T.G.Walker等從理論上給出了核磁共振陀螺中惰性氣體原子系綜的輸出頻率的誤差方程[12]，并從各個方面分析了各影響因素帶來的頻率誤差量級。然而，此方程沒有考慮頻率跟蹤方案可能引入的誤差，目前也未見有將其考慮在內的頻率誤差方程。本文在假設頻率跟蹤方案理想工作的前提下，著重分析了通過測量相位跟蹤核磁共振陀螺Larmor進動頻率所引入的額外頻率誤差，并建立相應的誤差方程。通過使用建立的誤差方程進行仿真分析，得到各因素對額外頻率誤差的影響，最終找到抑制此頻率誤差的方案。

1 理論

考慮系統設置如圖1所示，X軸施加激勵射頻場Bx=2B1cos(ω1t+β1)，其中B1為激勵射頻場有效幅值，ω1和β1為對應的圓頻率和初始相位；Z軸磁場條件為Bz=B0+Bccos(ωct)，B0為Z軸靜磁場強度，Bc為Z軸高頻載波場強度，ωc為對應的圓頻率。泵浦光極化方向沿Z軸正向，探測光沿X軸。PBS為偏振分光棱鏡，PD為光電二極管，BPD為由2個PD組成的平衡探測器結構。原子氣室內含有堿金屬原子87Rb，惰性氣體原子129Xe和131Xe以及緩沖氣體。

圖1 核磁共振陀螺系統設置框圖Fig.1 Systematic diagram of the nuclear magnetic resonance gyroscope

在這樣的條件下，對原子系綜進行基于Bloch方程的理論分析如下，系統的矢量Bloch方程為[13]：

(1)

(2)

可得Rb的橫向磁化強度如下[14]：

(3)

其中，Jn(·)表示n階第一類Bessel函數，n為整數，p為正整數。根據共振條件γRbB0+nωc=0，設置載波場的頻率,使其滿足ωc≈γRbB0，可得觀測的是n=-1條件下的磁化強度分量，則

(4)

(5)

使用鎖相放大器對其進行基于一次諧波的解調低通，有

(6)

(7)

若進行更為細致的考慮，實際應用中通常只能設置ωc≈γRbB0，即γRbB0-ωc=0并不精確成立，而存在一個失諧量Δωz=γRbB0-ωc，則式(4)可變為

(8)

這里m表示修改的。同樣式(5)也對應地修改為

(9)

得到基于一次諧波的解調低通信號

(10)

此時鎖放解調低通后的X及Y輸出通道的信號可表示為

FastXm=

(11)

FastYm=

(12)

考慮頻率跟蹤方案對探測到的信號(以X通道為例)進行估計跟蹤，可假設其擬合輸出的信號為：

FastXm=Acos(ωt+φ)+D

(13)

其中，A、ω、φ分別表示頻率跟蹤方案估計的幅度，圓頻率和初始相位，D為估計的直流分量。為了方便表達，令

b1=ΔωzT2

c1=[J-2(γRbBc/ωc)-J0(γRbBc/ωc)]sinφ

d1=[J-2(γRbBc/ωc)+J0(γRbBc/ωc)]cosφ

則式(11)又可表示為

FastXm=a1[(c1+b1d1)Bx+(b1c1-d1)By]

(14)

根據所加的橫向激勵場和Xe本身在橫向產生的磁場，可假設

Bx=2B1cos(ω1t+β1)+BXecos(ωXet+βXe)+Bx0

(15)

By=BXesin(ωXet+βXe)+By0

(16)

其中，BXe、ωXe和βXe分別表示Xe磁場的強度、圓頻率和初始相位，Bx0和By0分別表示系統感受到的磁場在X、Y兩軸上的分量直流部分。若假設頻率跟蹤方案能理想工作，則ω1=ω，β1=φ。根據式(13)～式(16)，有

Acos(ωt+φ)+D=a1{(c1+b1d1)[2B1cos(ω1t+

β1)+BXecos(ωXet+βXe)+Bx0]+

(b1c1-d1)[BXesin(ωXet+βXe)+By0]}

(17)

根據式(17)可得

D=a1[(c1+b1d1)Bx0+(b1c1-d1)By0]

(18)

[A-2a1(c1+b1d1)B1]cos(ωt+φ)=

a1BXe[(b1c1-d1)sin(ωXet+βXe)+

(c1+b1d1)cos(ωXet+βXe)]

(19)

sin(ωXet+βXe+βad)

(20)

其中，βad=arctan[(c1+b1d1)/(b1c1-d1)]，由式(20)可得

(21)

由于輸出的圓頻率值為所得相位在單位時間內的變化率，即ωoutput=d(ωt+φ)/dt，則有

(22)

可以看到，式(22)中等號右邊第一項為Xe自旋的本征圓頻率ωXe；第二項為Xe自旋初始相位的變化率ωβXe=dβXe/dt=0；第三項為鎖相環跟蹤頻率所額外引入的相位差的變化率ωβad=dβad/dt。前兩項的變化導致的頻率測量誤差在文獻[12]中已有討論，本文在不考慮前兩項變化的基礎上，著重考慮第三項，即鎖相環探測陀螺頻率所引入的誤差。

由于通常使用信號發生器或是鎖相放大器給線圈施加高頻載波場并進行解調，而儀器所發生的信號的穩定性遠高于系統內的其他參數。因此，可以假設Bc、ωc、φ為基本不隨時間變化的常數，僅有Z軸的靜磁場強度認為是時變量(橫向馳豫時間T2此處也看作常數)。則c1、d1為常數，根據βad的定義容易求得

(23)

又

(24)

則式(23)可變為

(25)

考慮ωβad的變化引起的頻率誤差，對式(25)求微分，又由f=ω/(2π)，得到微擾帶來的頻率誤差的表達式

(26)

通過式(26)即可預測B0不同的變化率及其二階導數，以及失諧量的函數b1，對于鎖相環引入頻率誤差的影響。

對于核磁共振陀螺來說，探測得到的Xe的Larmor進動圓頻率ωXe=γXeB+ωr,與其所感受到的靜磁場強度正相關，其中ωr為系統所感受到的轉動圓頻率。因此，在系統相對于慣性系靜止的情況下，若需要使探測到的頻率穩定，就需要使Xe感受到的靜磁場穩定。然而在實際應用中，雖然地磁場等環境磁場的波動可以通過設置高屏蔽效率的磁屏蔽來進行被動抑制約6～7個數量級[15]，同時也可通過主動的磁補償來進一步主動抑制約3～4個數量級[2]，磁場的影響仍然無法被忽略。為了抵消磁場帶來的影響，通常使用Xe的2個同位素129Xe和131Xe來同時探測他們的Larmor進動圓頻率[12,16]，如下所示：

ωXe129=γXe129B+ωr

(27)

ωXe131=γXe131B+ωr

(28)

其中，γXe129和γXe131分別為129Xe和131Xe的旋磁比。由式(27)～式(28)可以得到不受磁場影響的轉動圓頻率的表達式，即

(29)

通過雙同位素方案可以從理論上消除磁場對轉動頻率的影響，但是式(29)并沒有將相位引入的誤差考慮在內，若將其考慮在內，則2個同位素的輸出頻率分別為

ωout129=γXe129B+ωr+ωβad

(30)

ωout131=γXe131B+ωr+ωβad

(31)

可以看到，此時2個同位素有一個相同的額外頻率ωβad，若通過與式(29)相同的方式將磁場抵消，可得此時的轉動圓頻率的表達式

(32)

此時的轉動頻率就不可避免地受到式(26)的頻率誤差帶來的影響，因此，在使用頻率跟蹤方案探測陀螺頻率時，如何抑制這個誤差就顯得尤為重要。

2 仿真設置

考慮與實驗情況接近的參數設置，對頻率跟蹤方案探測陀螺頻率所引入的誤差進行具體的數值仿真分析,其中Rb的橫向馳豫時間采用文獻中的經典值[14]，實際實驗中可能會有不同的數值，但并不影響本文的分析。具體的基本參數設置如表1所示，其余的參數設置根據具體的仿真分析而定。

表1 參數設置Tab.1 Parameter settings

3 仿真分析

根據式(26)考慮不同的b1、dB0/dt、δB0和δ(dB0/dt)對δfβad的絕對值的影響，仿真結果如圖2～圖6所示，注意仿真結果圖的坐標軸均采用雙對數軸表示。這里磁場變化率及變化量的取值考慮目前電路水平一般能夠達到的范圍，并適當加大來考慮各因素的影響。

圖2 頻率誤差在不同dB0/dt條件下隨b1的變化曲線Fig.2 Curve of frequency error vs.b1 with different dB0/dt

圖3 EPR歸一化強度隨b1的變化曲線Fig.3 Curve of normalized EPR amplitude vs.b1

圖2和圖3在給定微擾值δB0=1nT和δ(dB0/dt)=0.01nT/s的情況下，考慮不同的b1和dB0/dt條件下|δfβad|的變化規律。其中圖2為|δfβad|隨b1的變化曲線，圖中各不同顏色的曲線分別表示磁場變化率dB0/dt=10-4,10-3,10-2,10-1,1(nT/s)。可以看到，在磁場變化率較小時(dB0/dt≤10-1nT/s)，|δfβad|隨b1的增加而呈現下降趨勢；而當磁場變化率較大時(dB0/dt>10-1nT/s)，|δfβad|隨b1的增加則總體呈現先增大后減小的趨勢。值得注意的是，當dB0/dt>10-1nT/s時，曲線在b1的不同取值處出現了局部的極小值，最小值甚至達到了10-10Hz量級。雖然這樣的結果確實給出了大幅抑制甚至消除額外頻率誤差δfβad的可能性，但是從式(26)中不難發現，這樣的極小值的出現是由于二階微擾項δ(dB0/dt)的存在，而實際上要精確控制磁場的二階微擾項是極其困難的。因此，通過這樣的極值點來抑制或消除δfβad并不可靠。另一方面，注意到當b1≥1時，所有曲線均呈下降趨勢，因此，通過選取一個相對較大的b1值，可以至少把δfβad抑制1個數量級。但是b1值并不是越大越好，因為它越大，就表示與n=-1的共振點越偏離，從式(11)可知，過大的b1值會造成信號衰減過大。因此，需要根據信號強度和頻率誤差選擇一個合適的值。由式(11)對電子順磁共振(Electron Paramagnetic Resonance,EPR)強度隨b1值的變化進行了仿真，仿真結果如圖3所示，其中對強度進行了歸一化處理。

從圖3中可以看出，隨著b1的增大，EPR歸一化強度開始基本不變，直到b1≈1時，則開始迅速下降。若將EPR歸一化強度下降1個量級作為信號探測的邊界條件，即超過1個量級即無法探測，則要求b1≤5。另一方面，從圖2中可以看出，當b1≈5時，|δfβad|相對于b1≈0時約有1個數量級的降低；b1≈10時更是有約2個數量級的降低；考慮到信號強度衰減b1取值在1～5之間較為合適，b1=5則能在實現信號探測的前提下最大程度地降低頻率誤差，其通過設置合適的ωc或B0則很容易實現，同時相對于b1≈0時的抗干擾能力也強很多。在此基礎上，通過降低磁場的擾動可以進一步降低頻率誤差δfβad。

圖4所示為|δfβad|隨dB0/dt的變化曲線，各不同顏色的曲線分別表示為b1=10-3,10-2,10-1,1,10。除去與圖2類似的局部極小值之外，不同于圖2的是，|δfβad|隨dB0/dt的增長基本呈現平穩的趨勢，僅在磁場變化率較大時(dB0/dt≥0.01nT/s)才開始變化，其下降的趨勢是由于局部最小值的存在，不考慮局部最小值則其基本保持不變或呈現上升趨勢(b1=1)。因此，只要將磁場變化率控制在較小的量級，其對額外頻率誤差的影響基本是不變的。對比圖4中不同的曲線可以清楚地看到，b1=10比b1=10-3時|δfβad|降低了約2個數量級，這與圖2的結果也完全符合。

圖4 頻率誤差在不同b1條件下隨dB0/dt的變化曲線Fig.4 Curve of frequency error vs.dB0/dt with different b1

圖5和圖6則在給定條件b1=5和dB0/dt=0.01nT/s的情況下，考慮不同的微擾值δB0和δ(dB0/dt)對 |δfβad|的影響。其中圖5為|δfβad|隨δB0的變化曲線，圖中各不同顏色的曲線分別表示靜磁場二階微擾量δ(dB0/dt)=10-6,10-5,10-4,10-3,10-2(nT/s)。與圖4類似，除去局部極小值外，在二階微擾量較小時(δ(dB0/dt)≤0.01pT/s)，|δfβad|隨δB0的增加而基本不變；在二階微擾量較大時(δ(dB0/dt)>0.01pT/s)，則呈現先平穩后單調遞增的趨勢。目前的電路水平一般至少能夠將電流導致的磁場變化抑制在0.01nT/s左右，而變化率的微擾量則更低，至少可以達到1pT/s的水平。注意到隨著磁場的二階微擾量δ(dB0/dt)降低1個數量級，|δfβad|也對應降低1個數量級，這從式(26)中也容易得到。因此，通過降低磁場的二階擾動，可以進一步地有效降低額外頻率誤差δfβad。

圖5 頻率誤差在不同δ(dB0/dt)條件下隨δB0的變化曲線Fig.5 Curve of frequency error vs.δB0 with different δ(dB0/dt)

圖6 頻率誤差在不同δB0條件下隨δ(dB0/dt)的變化曲線Fig.6 Curve of frequency error vs.δ(dB0/dt) with different δB0

圖6所示為|δfβad|隨δ(dB0/dt)的變化曲線，各不同顏色的曲線分別表示靜磁場微擾量δB0=10-4,10-3,10-2,10-1,1(nT)。同樣不考慮局部極小值的情況下，|δfβad|隨δB0的增加而呈現先平穩后單調遞增的趨勢，這是由于式(26)中δfβad與磁場的一階微擾量δB0和二階微擾量δ(dB0/dt)均呈線性關系。因此，可以看到，在較小的靜磁場二階微擾量的條件下，通過抑制磁場的一階擾動，同樣能起到降低額外頻率誤差δfβad的作用。

通過圖2～圖6的分析，可以做一個簡單的計算，首先在b1=5和dB0/dt=0.01nT/s的條件下，假設一個通過簡單控制就易于達到的磁場擾動值δB0=1nT及靜磁場二階微擾量δ(dB0/dt)=1pT/s，可得|δfβad|=0.224μHz，對于高精度陀螺(導航級)，這樣的誤差是無法忽略的；而在同樣的條件下，如果通過抑制磁場擾動值使其達到δB0=1pT以及相同的靜磁場二階微擾量δ(dB0/dt)=1fT/s，可以得到|δfβad|=0.224nHz，這樣的頻率誤差換算成轉速為2.9×10-4(°)/h，已經遠小于導航級的轉速誤差量級，可以認為對于導航級的陀螺其誤差可以忽略。

4 結論

1)本文基于Bloch方程推導了核磁共振陀螺輸出頻率的表達式，考慮了頻率跟蹤方案的探測相位引入的額外頻率誤差，并通過頻率誤差方程著重對引入的額外頻率誤差進行了仿真分析。

2)仿真結果表明，頻率跟蹤方案的探測相位引入的額外頻率誤差主要受到Z軸靜磁場強度及其一階變化率的影響。一方面通過選擇合適的Z軸靜磁場強度或載波場頻率使得其產生一定的失諧量，可以至少將引入的額外頻率誤差抑制1個數量級；另一方面，通過減小Z軸靜磁場的一階及二階擾動，其也會得到對應的抑制，結合兩者可以將額外頻率誤差抑制4個數量級至亞nHz量級，在導航級陀螺精度時已可以忽略。

3)本文首次將頻率跟蹤方案的影響考慮在核磁共振陀螺的頻率誤差方程內，目前并未見有考慮了測頻方案影響的頻率誤差方程的報道。

4)本文并未考慮堿金屬原子的橫向馳豫時間隨時間變化的情況，而這一項的變化也會對額外頻率誤差產生貢獻，下一步將對其進行深入分析。