多普勒非對稱空間外差儀探測高層大氣風速的三種方法比較研究*

李文文 惠寧菊 李存霞 劉洋河 方妍 李凌青 王彥龍 唐遠河

(西安理工大學理學院,西安 710048)

1 引言

地球中高層大氣的風速、風速結構和變化特性與空間航空環(huán)境、遙感探測水平和日地物理等密切相關[1].隨著火星、土星等行星大氣風速的遙感探測熱潮,高精度探測大氣風速已成為研究熱點[2].高層大氣風速含在氣輝輻射的多普勒頻移中,通過廣角邁克耳孫干涉儀(Michelson interferometer,MI)等光學儀器的成像干涉條紋可提取視線方向上的風速.為了提高測風儀器的穩(wěn)定性和測風精度,世界范圍內的科學家不斷努力發(fā)展了多種星載和地基光學測風儀器,努力提高測風精度.WINDII(wind imaging interferometer)探測地球上空80—300 km的大氣風速是通過壓電陶瓷驅動廣角MI一臂的動鏡步進4次以實現(xiàn)“四強度法”測風,測風精度達到10 m/s[3];PAMI (polarizing atmospheric Michelson interferometer)利用廣角MI中在一個周期內的偏振態(tài)變化4次實現(xiàn)“四強度法”測風,測風精度達5 m/s[4];WAMI (waves Michelson interferometer)將MI兩臂的反射鏡都固定,把一臂反射鏡分為4個區(qū)域再鍍反射膜,分別產生0,λ/4,2λ/4,3λ/4 (λ為波長)的步進光程差,實現(xiàn)“四強度法”測風[5,6].WINDII,PAMI和WAMI等儀器均產生干涉圓條紋.多普勒非對稱空間外差光譜儀(Doppler asymmetric spatial heterodyne spectroscopy,DASH)將MI兩臂的反射鏡分別換成閃耀光柵,產生干涉直條紋,采用傅里葉變換法探測地球上空峰值高度250 km的O(1D) 630.0 nm氣輝,室內測風精度為1.6 m/s[7];MIGHTI (Michelson interferometer for global high-resolution imaging)也采用DASH思路將MI兩臂的兩個閃耀光柵的衍射光再相干疊加,采用傅里葉變換法探測高層大氣風速,測風精度達1.2—4.7 m/s[8,9].中國科學院西安光學精密機械研究所的陳潔婧等[10]分析了DASH傅里葉變換法測風中窗函數(shù)的選取對測風精度的影響,中國科學院光電技術研究所的彭翔等[11]對DSAH測風的復合光程差相移解算方法的前提也是由傅里葉變換法從干涉圖中提取位相差.寧通[12]將DASH的干涉圖用傅里葉級數(shù)擬合而提取風速,測風精度為12.6 m/s.所以目前從DASH干涉條紋中提取風速有傅里葉級數(shù)法[12]和流行的傅里葉變換法[7].

本課題組研制的地基氣輝成像干涉儀GBAII(ground based airglow imaging interferometer)成功探測地球上空90—100 km 的大氣風速、溫度、密度等信息[13-15],通過改進,利用“四強度法”測風精度能達到4—6 m/s[16].本文將GBAII的兩臂換成閃耀光柵,改造為GBAII-DASH系統(tǒng),首次提出用“四強度法”反演GBAII-DASH的大氣風速.比較研究DASH測風的3種方法——傅里葉級數(shù)法、傅里葉變換法和“四強度法”獲取風速的原理、儀器正演和野外拍攝O(1S) 557.7 nm氣輝的風速數(shù)據反演及測風精度.

2 DASH探測高層大氣風速的3種方法

2.1 DASH測風關鍵—提取干涉條紋位相差

DASH光路如圖1所示,待測氣輝經準直透鏡Lens1后入射至立方分束器(beam splitter,BS)分成兩束光強相同的相干光,經閃耀光柵G1和G2衍射后原路返回BS,再到CCD上成像干涉.軸向光以θ角入射到光柵上,如果光線以θ角沿原方向衍射回來,此時波數(shù)稱為光柵的Littrow波數(shù)(σL),此角θ稱為Littrow角.如圖2所示的平面反射式光柵的夫瑯禾費衍射光強為[17]

圖1 DASH的光路圖[7]Fig.1.Optical path diagram of DASH.

圖2 閃耀光柵結構圖Fig.2.Blazing grating structure diagram.

其中a是光柵小反射面的寬度;d是光柵常數(shù);i和i'是光線關于小反射面法線n的入射角和衍射角,φ和φ'是光線關于光柵法線N的入射角和衍射角;i,i',φ和φ'的正負規(guī)定按光線以銳角轉向法線的順時針為正、逆時針為負.πa(sini+sini')/λ是單縫衍射因子,πd(sinφ+sinφ')/λ是縫間干涉因子.縫間干涉確定主極大峰值強度、位置和數(shù)目,滿足光柵方程d(sinφ+sinφ')=kλ (k=0,±1,±2,···)時衍射產生主極大,考慮反射角與衍射角的正負,|sinφ+sinφ'|不可能大于1,這就限制了主極大數(shù)目.主極大的最大級|k| ＜d/λ,故一般選取第1級衍射.

入射到DASH光束的波數(shù)若是σL,過BS后的兩光束經光柵衍射后返回的出射波面都與光軸垂直,相位差為0;若入射光波數(shù)不是σL,兩出射波面的傳播方向與光軸都有一夾角±γ,形成斐索干涉條紋,此時φ=θ,φ'=θ -γ,將其代入光柵方程,便可求出γ.斐索干涉條紋經透鏡Lens2后成像于CCD靶面,此時斐索干涉條紋的空間頻率為

若高斯線型氣輝入射到DASH上,其干涉方程為

其中,光程差Δ=4tanθ[x+ΔL/(2tanθ)] (單位為cm),x是CCD探測器的像素點位置坐標,ΔL是兩臂的路徑差,如圖1所示,2ΔL是干涉儀的固定光程差.B(σ)=B0exp[-4ln2(σ-σ0)2/ω2] 是光源輻射譜,B0=(2/ω)/(ln2/π)1/2,零風速波數(shù)σ0與波長對應關系σ0=1/λ0,ω=[(7.16×10-7)2×]1/2是高斯線型氣輝的半高寬,其中T是熱平衡時的大氣溫度(單位K),M是發(fā)光物質的原子量.

如果視線方向的氣輝粒子對探測器有相對運動速度v,則波數(shù)變?yōu)棣?σ0(1+v/c),根據(3)式及多普勒頻移,得到含大氣風速項的相位φv:

零風速相位φ0為

由(4)式可知,DASH干涉儀測風的關鍵是探測出干涉條紋的相位差.國際上采用傅里葉級數(shù)和傅里葉變換兩種方法從DASH干涉圖中提取相位量來反演大氣風速.

2.2 傅里葉級數(shù)法測風

對單一譜線入射的氣輝光源,因為傅里葉級數(shù)可轉變成余弦函數(shù),積分變?yōu)榍蠛?由(3)式可知干涉強度I是像素點x的余弦函數(shù)[12]:

首先對DASH干涉圖進行傅里葉擬合,得出傅里葉級數(shù)方程:

根據傅里葉級數(shù)三角變換,將(7)式轉換成如下余弦函數(shù)[12]:

其中ci是權重系數(shù);φi可由(8)式中的系數(shù)求得,φi=-arctan(bi/ai);f是頻率.余弦函數(shù)(8)中包括i個不同頻率的余弦項,φi是第i項余弦項的初始相位.

當被拍攝的氣輝譜線是單一波長時,可用一級諧波作為含該主要信息的數(shù)學模型,即當i=1時,(8)式中的φ1項是(6)式中的4π(σ -σL)ΔL項.當干涉圖是零風速的光所成像時,φ1便是(5)式所示的零風速相位φ0;當干涉圖是含有多普勒頻移的光所成像時,φ1便是(4)式和(5)式中有大氣風速相位和零風速相位之和φ1=φv+φ0.傅里葉級數(shù)法分別對零風速和含風速干涉圖的數(shù)據進行傅里葉級數(shù)擬合,求出v=0的φ1及風速v≠ 0的φ1,兩者相減得到φv,結合(4)式即可求出視線方向上的風速.

2.3 傅里葉變換法測風

當DASH干涉儀接收連續(xù)譜線簡化為幾條離散光譜時,則(3)式寫成[7]:

其中j指的是第j條譜線,Sj與探測器上干涉圖條紋亮度成正比,Ej(x)依賴光譜線型和光程差的包絡函數(shù),κj=4(σj-σL)tanθL是每條譜線中心σj的外差空間條紋頻率,φj=4π(σj-σL)ΔL是固定光程差2ΔL引起的相位,δφj是每條譜線多普勒頻移導致的相位頻移量.

選擇合適的光譜帶、Littrow波數(shù)和ΔL后,(9)式的傅里葉變換產生一個復頻譜,在空間頻率κj和-κj附近有局部分離良好特性.分離出這些特性,如令j=0,通過歸零除了+κ0附近局域內的光譜元素,有效消除了所有干涉圖貢獻,僅保留一個指數(shù)項,因此在傅里葉逆變換后得到[7]

相位項可從(10)式的虛部和實部之比得到

將(11)式計算結果減去零風速相位2πκ0x+φ0,便得到相移δφ0,根據(4)式求出該譜線在視線方向上的風速,這就是國際上常用的傅里葉變換測風方法.

2.4 “四強度法”測風

WINDII用壓電陶瓷驅動MI的“動鏡”實現(xiàn)相位的4步進測風[3],而DASH結構中沒有移動部件,我們提出也可以用“四強度法”測風,闡述如下.如圖3所示的簡化GBAII-DASH光路,非Littrow波長光入射時,兩出射光的波面間會有夾角2γ,如圖4所示.兩出射光到達CCD相同像素點時,若兩者的光程差正好是一個λ的長度,或者是一個λ的整數(shù)倍時,探測器CCD靶面上會形成亮條紋.故對于GBAII-DASH來說,我們用兩出射光到達CCD相同像素點的光程差來實現(xiàn)相位的四步進法測風.

圖3 GBAII-DASH的光路圖Fig.3.Optical path diagram of GBAII-DASH.

圖4 兩出射波面夾角示意圖Fig.4.Angle of emergent wave surface.

高斯線型氣輝經廣角MI后所成干涉圖的強度為[3]

其中σ0是零風速時的波數(shù),T是大氣溫度,Q=(π2ω2)/4Tln2,是高斯線型氣輝的半高寬,M是以σ0為中心的發(fā)射線的物質的原子量.令條紋調制度V=exp(-QTΔ2),則

根據多普勒效應σ=σ0(1+v/c),將(13)式的σ0換成σ,光程差Δ分為含固定光程差Δ0、含風光程差Δv、步進光程差Δ′之和Δ=Δ0+Δv+Δ′(且Δ0?Δv,Δ′),則由(3)式可得Δ0=2ΔL,Δv=Δ0(v/c),Δ′=4xtanθ,3種光程差對應相位分別為φ0,φv,φ'.因此(13)式變?yōu)?/p>

當φ'=0,π/2,π,3π/2 (對應Δ′=0,λ/4,2λ/4,3λ/4,條紋周期很小,量級是10-5,為了防止步進太小而被吞掉,會以(π/2+2kπ)(k是整數(shù))方式步進),代入(14)式得4個強度值:

由此得到

由(18)式得到φ0+φv后,減去零風速定標所確定的φ0(定標選擇cos2πσ0Δ0=1條件得Δ0,亦即φ0=2kπ (k為整數(shù))),就得出含有風速的相位φv,再根據(4)式即可求出視線方向上的風速,這便是我們提出的GBAII-DASH的“四強度法”測風法.

3 三種測風方法的正演結果

利用計算機仿真正演GBAII-DASH干涉儀的光學成像過程,需要考慮大氣傳輸模型和儀器模型等諸多部分,這里只對大氣傳輸模型及儀器模型進行比較研究.

本系統(tǒng)以氧原子O(1S) 557.7 nm作為目標譜線,其體發(fā)射率分布為[18]

式中Ve,Vf代表地球上空E層(100—300 km)和F層(300—500 km)的峰值體發(fā)射率(單位photons·cm-3·s-1),χ是太陽天頂角(單位(°)),be=(h-He)/We,bf=(h-Hf)/Wf,h是海拔高度,He和Hf分別是E層和F層的峰值高度,We和Wf為E和F層氣體的標高.

GBAII-DASH探測系統(tǒng)使用窄帶濾波片的濾波函數(shù)與波長λ、入射角θ的關系為[19]

其中Δλ是半高寬,λ0是中心波長,n是濾光片的有效折射率.

(14)式中的φ0可由零風速定標而確定,選擇適當?shù)墓潭ü獬滩?ΔL,令φ0=2kπ (k為整數(shù)),而不出現(xiàn)在余弦函數(shù)中,則展開干涉強度函數(shù)(14)式得

令J1=I0,J2=I0Vcosφv和J3=I0Vsinφv,則CCD的第l行j列像素上獲得的模擬結果為

其中R是儀器響應度,2ΔLij是固定光程差,Nnoise是輸出信號中存在的噪聲.

儀器正演仿真時我們選取CCD是1024 × 1024面陣探測器,單像素尺寸24 μm × 24 μm,固定光程差2ΔL=7.495 cm;光柵Littrow波長550 nm,Littrow角14.3°,刻線密度900 L/mm.模擬西安上空峰值高度為98 km的O(1S) 557.7 nm氣輝.假設一個大氣風速值,得到正演仿真干涉圖后,分別用傅里葉級數(shù)法、傅里葉變換法和“四強度法”從干涉圖中提取風速.

3.1 傅里葉級數(shù)法測風的正演結果

干涉條紋進行一級傅里葉級數(shù)擬合后的結果如圖5所示,圖6是圖5的局部放大圖,從圖6可明顯觀察到,因為風速而導致的干涉圖相位的頻移.

圖5 干涉圖的傅里葉級數(shù)正演結果Fig.5.Fourier series forward results of interferograms.

圖6 圖5中正演結果的局部放大Fig.6.Local amplification of forward results in Fig.5.

假設以10 m/s為風速間隔的0—100 m/s風速,得到正演干涉結果,對風速為0 m/s干涉結果進行傅里葉級數(shù)擬合后的方程式如下:

同理可得出其他風速的傅里葉級數(shù)擬合方程.根據(7)式和(8)式可求出不同風速的相位φ,此時0風速相位φ0就是GBAII-DASH系統(tǒng)固定光程差2ΔL導致的相位,則不同風速導致的相位頻移量φv=φ -φ0,根據公式v=cφv/4πΔLσ0求出風速如表1第2列所示,平均測風相對誤差是2.98%.

表1 三種測風方法的正演結果Table 1.Forward wind speed results by three methods.

3.2 傅里葉變換法測風的正演結果

對干涉數(shù)據進行傅里葉變換后,利用(11)式分別求出不同風速的干涉圖上各像素點的相位,如圖7所示.選取第387個像素點的相位進行分析,即x=387時,不同風速的相位Φ=2πκx+φ0+δφv.而0風速的相位就是系統(tǒng)固定光程差2ΔL導致的相位Φ0=2πκx+φ0,則不同風速導致的相位頻移量δφv=Φ -Φ0,用v=cδφv/4πΔLσ0求出風速如表1第3列所示,平均測風相對誤差是4.67%.

圖7 傅里葉變換的相位分布圖Fig.7.Phase distribution diagram of Fourier transformation.

3.3 “四強度法”測風的正演結果

用“四強度法”反演風速,首先需確定如圖4所示起點I0的坐標x0,然后以兩出射光到達CCD像素點x0的光程差為起點,令光程差依次步進λ/4或者kλ+λ/4 (k是整數(shù)),故x0的選取很重要.干涉圖的強度I是關于CCD像素點坐標x的函數(shù),將GBAII-DASH干涉儀中固定光程差2ΔL所導致的相位φ0,代入到干涉圖方程中,求出相位是φ0的像素點坐標x,則此x便是步進的起始點坐標x0.

假設風速為50 m/s得到仿真干涉數(shù)據的擬合函數(shù)如圖8的綠色曲線所示,因為干涉圖的周期很小,一個周期量級為10-5,為3個像素大小,每次步進1/4個周期,即3/4個像素的距離,這個距離太小而被吞掉,故我們將周期性函數(shù)進行拉伸,得到拉伸后的“四強度法”拉伸函數(shù)結果如圖8的紅色虛線所示.經過拉伸后再步進,確定出x1,x2,x3,x4,各像素點對應的強度值I1=0.1793,I2=0.8836,I3=0.8207和I4=0.1164,利用tan(φ0+φv)=(I4-I2)/(I1-I3)和反正切求出v=50 m/s時φ=0.8745 rad.

圖8 函數(shù)拉伸后的“四強度法”Fig.8.Four steps of phase determination.

用同樣的方法處理其他風速的干涉圖數(shù)據,得出相位.將0風速的相位φ0作為GBAII-DASH固定光程差2ΔL所對應的相位,則風速導致的相位頻移量是φv=φ -φ0,利用v=cφv/4πΔLσ0求出正演風速如表1第4列所示,平均測風相對誤差是3.00%.

3.4 三種測風方法正演結果的誤差

上述正演研究中未考慮噪聲的影響,但是實際拍攝氣輝的成像干涉圖存在多種噪聲,從干涉圖提取相位之前需對原始數(shù)據去噪和平場.為了探測噪聲和平場對3種方法測風誤差的影響,我們在正演過程中對上述風速的干涉圖人為添加均值為0、標準差為0.1的高斯噪聲,對各數(shù)據進行平場處理后,利用傅里葉級數(shù)、傅里葉變化和“四強度法”進行正演,得到正演風速結果如表2所示,3種方法得到的平均相對誤差分別為2.30%,11.66%,2.27%.可見干涉圖存在噪聲時,傅里葉變換法的測風誤差相對較大,傅里葉級數(shù)法和“四強度法”測風的測風精度高,但兩種方法在測風間隔是10 m/s時的測風誤差相近,為了更好地區(qū)分出傅里葉級數(shù)法和“四強度法”的測風精度,我們繼續(xù)分析測風間隔在1 m/s和0.1 m/s時兩種方法的測風誤差.

表2 加入噪聲后的3種測風方法的正演誤差Table 2.Speed Error after adding noise by three methods.

模擬含有噪聲的以1 m/s為間隔的風速是31—39 m/s和以0.1 m/s為間隔的風速是30.1—30.9 m/s的干涉圖,用傅里葉級數(shù)法和“四強度法”從干涉圖中提取風速,求出誤差,結果如圖9所示.從圖9可以看出,1 m/s為風速間隔時“四強度法”的平均測風相對誤差是2.20%,明顯低于傅里葉級數(shù)法的平均測風相對誤差3.55%,以及0.1 m/s為風速間隔時“四強度法”的平均測風相對誤差是2.69%,也明顯低于傅里葉級數(shù)法的平均測風相對誤差4.15%.故“四強度法”的測風精度更優(yōu)于傅里葉級數(shù)法測風.

圖9 兩種方法的測風誤差Fig.9.Wind measurement error of two methods.

從正演和增添噪聲的上述結果可見,傅里葉變換法測風的誤差較大.傅里葉變換法測風時,將窗函數(shù)頻譜卷積過程中導致復原干涉圖的包絡與理想干涉圖的包絡在小光程差和大光程差區(qū)域發(fā)生明顯變形,直接影響干涉相位的計算,導致測風誤差大.減小誤差的辦法是選擇中心區(qū)域的光程差點能較準確反演出風速[20].

“四強度法”測風時,是通過(18)式的反正切求出相位,其中(I4-I2)和(I1-I3)已經把干涉圖的背景噪聲和直流部分統(tǒng)統(tǒng)減掉,故對數(shù)據平場后測風精度較高;況且“四強度”法測風計算簡便,不用考慮窗函數(shù)因子的不確定性帶來的測風誤差.

4 室外實驗3種方法的反演結果

采用圖3所示的GBAII-DASH光學系統(tǒng)(實物見圖10),于2023年4月19日凌晨1—3點在陜西西安臨潼洪慶山頂(海拔1250 m,34°19′52′′ N,109°16′56′′ E)拍攝O(1S) 557.7 nm氣輝的成像干涉圖如圖11所示,圖11(a)是GBAII-DASH的0°天頂角所拍的結果,用于零風速定標,圖11(b)是天頂角為45°時拍攝的O(1S)氣輝成像干涉圖.

圖10 GBAII-DASH 的實驗系統(tǒng)Fig.10.GBAII-DASH system in the laboratory.

圖11 GBAII-DASH拍攝O(1S) 557.7 nm氣輝的成像干涉圖 (a) 0°天頂角時拍攝的干涉圖;(b) 45°天頂角時拍攝的干涉圖Fig.11.Imaging interferogram of O(1S) 557.7 nm gas glow obtained by GBAII-DASH: (a) Interferogram taken at 0°zenith angle;(b) interferogram taken at 45° zenith angle.

通過去噪和平場等措施后,分別用傅里葉級數(shù)法、傅里葉變換法和“四強度法”提取圖11氣輝干涉圖的大氣風速如表3所列.分別用傅里葉級數(shù)法、傅里葉變換法和“四強度法”測得當晚西安上空98 km的一維風速為32.21 m/s,43.55 m/s和32.17 m/s.

表3 三種方法反演室外測風結果Table 3.Inversion wind speed outdoor experiment by three methods.

5 結論

基于DASH探測高層大氣風速,比較研究了傅里葉級數(shù)法、傅里葉變換法和“四強度法”提取風速的原理、正演、噪聲和反演等內容,3種方法探測風速均從DASH光學系統(tǒng)所得斐索干涉條紋的相位差變換而來,結論如下:

1)模擬以10 m/s為風速間隔的風速0—100 m/s的正演斐索干涉圖,用傅里葉級數(shù)法、傅里葉變換法和“四強度法”得到正演風速值,計算得到3種方法的平均測風誤差分別為2.93%,4.67%和3.00%.

2)人為添加均值為0、標準差為0.1的高斯噪聲后,假設風速是0—100 m/s,用傅里葉級數(shù)、傅里葉變換和“四強度法”分別對平場后的數(shù)據進行正演,得到平均相對誤差分別為2.30%,11.66%,2.27%.

3)以1 m/s為間隔模擬風速31—39 m/s,得到含高斯噪聲的正演斐索干涉圖,用傅里葉級數(shù)法和“四強度法”得到正演風速值,并得到平均測風誤差分別為3.55%,2.20%;以0.1 m/s為間隔的風速30.1—30.9 m/s,模擬得到含高斯噪聲的正演斐索干涉圖,用傅里葉級數(shù)法和“四強度法”的平均測風誤差分別為4.15%,2.69%;“四強度法”的測風誤差都小于傅里葉級數(shù)法的測風誤差.

4)利用GBAII-DASH拍攝西安上空98 km的O(1S) 557.7 nm氣輝,得到天頂角為0°和45°的成像干涉圖,再用傅里葉級數(shù)、傅里葉變換和“四強度法”得到風速的反演結果分別為32.21 m/s,43.55 m/s和32.17 m/s.

5)從DASH的正演、反演數(shù)據結果看,我們提出的“四強度法”探測高層大氣風速的結果較好,計算簡便且測風精度相對較高.