星光成像的大氣影響研究（Ⅱ）：大氣湍流

陶志煒，戴聰明，武鵬飛，任益充，梅海平，童杰，5，駱傳凱，青春，馮云松，魏合理，饒瑞中

（1 中國科學院合肥物質科學研究院 安徽光學精密機械研究所 中國科學院大氣光學重點實驗室， 合肥 230031）

（2 國防科學技術大學 電子對抗學院 紅外與低溫等離子體安徽省重點實驗室， 合肥 230037）

（3 先進激光技術安徽省實驗室， 合肥 230037）

（4 合肥工業(yè)大學 物理學院， 合肥 230601）

（5 中國科學技術大學研究生院科學島分院， 合肥 230026）

0 引言

人類利用恒星來進行導航，最早可追溯到古代人們通過北極星來確定方位。直到20 世紀50年代，星敏感器的橫空出世，大大提升了恒星導航的精度。星敏感器是一種高精度的姿態(tài)敏感測量儀器，它通過成像系統(tǒng)對星空成像，測量恒星矢量在星敏感器坐標系中的分量，利用已知的恒星精確位置來確定載體相對于慣性坐標系的三軸姿態(tài)［1］。星敏感器導航技術的高精度、抗干擾性強、可不依賴其它系統(tǒng)進行獨立導航等優(yōu)點在近地空間各類機載、艦載、車載平臺上有著廣泛的應用［2］。

早期，星敏感器應用于衛(wèi)星平臺，承擔了衛(wèi)星姿態(tài)測量的任務，是衛(wèi)星平臺不可或缺的測量設備。由于大氣層外可近似看成是真空環(huán)境，因此，星載星敏感器在工作時幾乎不受大氣的影響，其測姿精度通?？蛇_到角秒級甚至亞角秒量級。隨著空天觀測平臺的發(fā)展以及大氣層內(nèi)星敏感器觀測高度的降低，大氣對恒星探測的影響愈發(fā)顯著，星敏感器在成像時不可避免的會受到地球大氣背景輻射、湍流、折射等因素的影響。白天大氣分子與氣溶膠粒子會對太陽光產(chǎn)生散射，使得天空產(chǎn)生復雜的背景光，對星敏感器在白天觀星時產(chǎn)生極大的干擾，嚴重降低成像的信噪比，使得恒星目標湮沒在背景中，無法識別［3］。另一方面，大氣分子與氣溶膠粒子的存在會使得星光傳輸產(chǎn)生衰減，再加上大氣本身存在折射率起伏，還會影響光波的相位，并對振幅產(chǎn)生調制，引發(fā)星光的閃爍和抖動［4-10］。由于大氣折射率隨高度存在一定變化，因此當星光經(jīng)過大氣傳輸時，大氣折射會延長星光的傳輸路徑，改變星光的傳輸方向，嚴重影響星光導航的精度。

為了定量評估大氣對星光成像的影響，本系列文章之一著重研究了星光成像中恒星“看的見”這一重要問題。其次，為了進一步獲取恒星的位置采用星敏感器進行定位和導航，還需著重研究恒星“看的準”這一問題，即如何評估星光成像的大氣湍流和大氣折射影響。實際上，關于星光大氣湍流影響的研究最早可追溯至對星光圖像結構和運動的研究。這種運動在幾何光學角度看來最直接的影響就是使得星光產(chǎn)生一定程度的角向位移［4］。隨著光波隨機介質傳輸理論的建立［11］，研究人員才逐步弄清星光抖動產(chǎn)生的原因，即由于大氣湍流使得波前產(chǎn)生傾斜導致的。除此之外，大氣湍流會使星光產(chǎn)生另一個高階效應，即對星光的強度產(chǎn)生調制［12］（折射率的隨機起伏會對星光的能量產(chǎn)生耗散和起伏），進而使得星光發(fā)生“忽亮忽暗”的現(xiàn)象，這便是人們常說的閃爍效應［4-10］。

本文是星光成像大氣影響三部分研究（背景輻射、湍流及折射）的第二部分，目標在于構建星光大氣湍流傳輸?shù)臄?shù)理模型，并針對我國幅員遼闊特點，選取典型地區(qū)，獲取測試數(shù)據(jù)，優(yōu)化完善我國典型區(qū)域星光閃爍和抖動的數(shù)據(jù)模型。首先，結合大氣湍流廓線數(shù)據(jù)，獲得典型地區(qū)星光傳輸鏈路的相位屏分布特性和湍流特性。其次，結合分步式光傳輸方法，研究星光經(jīng)過大氣湍流傳輸后的閃爍量、到達角以及到達角起伏等物理量的變化規(guī)律，定量評估大氣湍流對星光成像的影響。本文的結果對于指導星敏感器的設計、星像質心算法的開發(fā)研究具有一定的指導價值。

1 星光大氣湍流傳輸模型

1.1 分步式光傳輸模型

星光成像大氣湍流仿真的核心問題在于研究大氣湍流對星光傳輸?shù)挠绊?。為了從?shù)值角度給上述問題一個定量的結果，首先介紹星光大氣湍流傳輸?shù)姆植绞焦鈧鬏斈Ｐ停渲饕枷胧菍⒋髿庑诺揽闯梢幌盗械耐牧飨辔黄恋慕M合，每個相位屏會在光束上引入一個滿足大氣湍流特性的隨機相位?，并將兩個相位屏之間的傳輸過程看成真空環(huán)境下的Fresnel 衍射［13-14］。值得說明的是，每個相位屏不會導致光束振幅發(fā)生改變，而是通過真空衍射不斷改變光束強度的分布情況，如圖1 所示。對于每一段的傳輸過程，可描述為［14］

圖1 星光大氣湍流傳輸?shù)脑韴D和星光傳輸?shù)姆植绞酱髿鈧鬏斈Ｐ虵ig.1 Schematic diagram of the split-step beam propagation model

1.1.1 相位屏生成

針對第一個問題，目前常用的構造相位屏的方法是基于傅里葉變換得到的，其基本原理是利用湍流折射率譜和復高斯隨機數(shù)矩陣產(chǎn)生相空間復隨機場，然后進行逆傅里葉變換獲得二維相位的空間分布，其過程可表示為［13，15-17］

需要說明的是，雖然采用傅里葉變換產(chǎn)生湍流相位屏的方法較為簡單，但使用該方法生成的相位屏明顯具有低頻不足的缺點，這是因為采用該方法生成的相位屏的最小頻率為Δfx，其中不包含和這部分頻率的貢獻。為了補償相位屏中低頻成分的貢獻，國內(nèi)外相關學者提出了一種次諧波補償?shù)南辔黄辽煞椒?，其?shù)學形式為［13，18］

式中，p為次諧波級數(shù)，且次諧波頻率間隔為Δfp=1/(3p L)。本文采用次諧波補償?shù)母道锶~變換法來改善相位屏模擬大尺度的統(tǒng)計特征。

1.1.2 湍流描述

針對本文研究內(nèi)容，首先需得到我國典型地區(qū)不同時段的湍流廓線，其次根據(jù)湍流廓線設置相應的相位屏分布。眾所周知，大氣湍流的分布具有一定的垂直結構，其通道內(nèi)的光學湍流強度變化可描述為折射率結構常數(shù)C2n隨高度h的函數(shù)。常見的C2n垂直廓線的參數(shù)化模式有HV 模式［19］，HMN 模式［20］以及Dewan 模式［21］等。然而，由于光學湍流的變化往往具有一定的地域和季節(jié)等特征，因此若盲目采用上述模式計算星光大氣湍流傳輸后的結果，會得到不符合我國局部的實際地理和大氣狀況的結果。因此，本文的計算和分析并未采用經(jīng)典的參數(shù)化模式，而是基于ERA5 數(shù)據(jù)1ERA5 是歐洲中期天氣預報中心（European Centre for Medium-Range Weather Forecasts，ECMWF）第五代全球大氣再分析數(shù)據(jù)集，提供的變量達到240 個，包含不同的大氣參量，如溫度、相對濕度、風速等，使用戶能精確地分析過去長時間范圍大氣的狀態(tài)。數(shù)據(jù)下載地址詳見：ERA5-Land hourly data from 1950 to present (copernicus.eu)。并結合光學湍流預報方法［22-23］，在得到的湍流廓線數(shù)據(jù)基礎上進行的。為了計算典型區(qū)域不同時刻的星光閃爍情況，本文選取了2008年1月1日三沙（地理位置：112.5°E， 17°N，溫度：21.6 ℃，濕度：73.2%，壓強：1 017.69 hPa）、大柴旦（地理位置：95.25°E， 37.75°N，溫度：?12.6 ℃，濕度：41.8%，壓強：665.79 hPa）以及合肥（地理位置：117.25°E， 31.75°N，溫度：0.81 ℃，濕度：50.9%，壓強：1 030.27 hPa）三地一天四個時間段的湍流情況進行計算，相應的湍流廓線如圖2 所示。

圖2 2008年1月1日三沙、大柴旦、合肥三地一天四個時間段的湍流廓線圖Fig.2 Turbulence profiles at Sansha， Dachaidan and Hefei for four time periods of one day on January 1， 2008

1.1.3 相位屏分布

根據(jù)以上得到的典型區(qū)域不同時刻的湍流垂直分布情況，便可設置星光大氣傳輸信道相應的相位屏分布。一般而言，對于折射率起伏均勻的光傳播路徑，一般設置等間隔的相位屏；而對于非均勻路徑下的光傳播，如地球大氣中的上行或下行傳播問題，由于傳播路徑上的折射率起伏強度不相等，若仍設置等間隔的相位屏，則會導致強起伏區(qū)不能被充分采樣，而弱起伏區(qū)被過度采樣。這不僅使得數(shù)值模擬不能真實地反映光傳播的實際過程，又增大了不必要的計算量。因此本文采用等Rytov 指數(shù)間隔相位屏分布方式來進行數(shù)值仿真［24］。具體原理可分為以下幾個步驟（詳見附件流程圖）。

1）設定常數(shù)c，使得相鄰兩層相位屏的Rytov 指數(shù)σ2R(Δzi)≡c。

2）根據(jù)近地面的C2n(0)以及，計算第一層相位屏對應的高度，用于設定初始值。

3）對于第i個相位屏到i+1 個相位屏的間距Δzi，采用第i層相位屏對應高度的C2n(zi)和恒等式進行計算。

4）迭代第3）步，直到相位屏間距總和等于總的傳輸距離，停止迭代。

上述采用等Rytov 指數(shù)間隔的相位屏分布是假定天頂角為零的條件下進行，若天頂角為θZ且上述生成的相位屏的高度分別為h1，h2，…，hn，則應調整相位屏分布的高度為h1secθZ，h2secθZ，…，hnsecθZ。

對于天頂角為θZ的各層相位相位屏的大氣相干長度r0i，采用式（4）進行計算。

1.2 星光閃爍

上一節(jié)介紹了星光大氣湍流傳輸?shù)姆植际焦鈧鬏斈Ｐ?，結合典型地區(qū)湍流廓線設置相應的相位屏分布，便可得到星光大氣湍流傳輸后的結果。根據(jù)已知星光大氣湍流傳輸后的電場，結合式（5），可計算星光的閃爍指數(shù)，它通常刻畫了接收端星光的強度起伏情況［25］

式中，符號I和分別表示光場強度和系綜平均值。一般而言，在任意起伏條件下，閃爍指數(shù)一般不能獲得一個明確的解析表達式，但對于Kolmogorov 譜情況下的平面波而言，任意湍流強度下垂直信道的閃爍指數(shù)可表述為［26-27］

式中，σ2R為Rytov 指數(shù)有

式中，h為觀測高度，θz表示天頂角的大小。需要說明的是，式（6）給出的結果是平面波垂直信道的點閃爍公式，而式（5）給出的計算公式具有一定的普適性，可計算任意孔徑區(qū)域內(nèi)的閃爍指數(shù)。一般而言，采用數(shù)值接收電場計算閃爍指數(shù)時總是會得到一定區(qū)域范圍內(nèi)閃爍指數(shù)平均的效果（這是由于網(wǎng)格的空間分辨率決定的）。眾所周知，閃爍指數(shù)一般與接收端望遠鏡接收孔徑大小成反比，因為較大的望遠鏡將平均孔徑上的波動。孔徑平均閃爍應小于用點接收器測量的閃爍。因此，為了將數(shù)值結果與理論結果進行對比，需對數(shù)據(jù)結果除以一個孔徑平均因子，才可與理論結果進行對比。對于平面波而言，孔徑平均因子可表示為

式中，D代表望遠鏡的接收口徑大小，有限面積上的任意起伏條件下閃爍方差σ2I(D)為［28］

1.3 遠場成像

根據(jù)接收電場不僅可以計算星光經(jīng)過大氣湍流后的閃爍值，還可計算星光經(jīng)望遠鏡遠場成像后的光斑的抖動量、到達角以及到達角起伏，后者通常是星光到達像面上傳播方向的度量。一般而言，星光經(jīng)過大氣湍流傳輸后被接收端望遠鏡接收，而后在遠場成像，其成像模型可等效為物體緊貼棱鏡的成像方式［13］。因此，星光的遠場光斑可表述為二維遠場衍射積分［13］

式中，uout(x1，y1)為星光經(jīng)大氣湍流傳輸后的電場，P(x1，y1)為望遠鏡的光瞳函數(shù)，f代表望遠鏡的焦距。根據(jù)式（10），可計算遠場光斑的抖動量。值得注意的是，由式（10）還可以看出，如若需要統(tǒng)計遠場光斑的抖動量，則需要事先知道成像系統(tǒng)的焦距f。換句話說，遠場光斑的抖動量與成像系統(tǒng)的焦距有關，但這一點是我們不希望的。因此，為了研究大氣湍流本身對于遠場光斑抖動的影響，本文采用到達角來描述星光的抖動，其可表述為

式中，Δr代表遠場光斑的抖動量。例如：若接收端望遠鏡的焦距為1 m 且到達角為2 μrad，則遠場光斑相應的抖動量為2 μm。另一方面，根據(jù)式（11）還可計算遠場光斑的到達角起伏，這一點也是本文所關注的。

2 星光大氣湍流影響分析

2.1 閃爍指數(shù)的變化

圖3 給出了不同觀測高度下（觀測高度為0 km、1.5 km、3 km、4.5 km、6 km、7.5 km）三沙、大柴旦、合肥一天四個時間段的星光閃爍隨天頂角的變化情況，其中不同形狀的曲線代表不同觀測高度下的閃爍結果，每幅圖中的插圖代表數(shù)值仿真結果與根據(jù)式（6）理論計算結果的對比情況。

圖3 不同觀測高度下三沙、大柴旦、合肥三地一天四個時間段星光閃爍指數(shù)隨天頂角的變化情況Fig.3 Scintillation index of starlight as a function of zenith angle in Sansha， Dachaidan and Hefei under different observation altitudes

從圖3 中插圖可以看出典型區(qū)域不同時間段數(shù)值仿真計算得到的星光閃爍指數(shù)基本與式（6）給出的理論結果相吻合，說明采用上述數(shù)值仿真方法計算星光抖動相關指標具有一定的可行性。除此之外，可以看到，當天頂角逐漸增加時，數(shù)值仿真與理論值誤差相對較大，這是因為隨著天頂角增大，傳輸距離逐漸變化，如果采用與天頂角較小時一樣的網(wǎng)格尺寸，便會產(chǎn)生一定的計算誤差。

對比圖3 中不同子圖結果可以發(fā)現(xiàn)：不同區(qū)域、不同時間段以及不同觀測高度下，星光閃爍指數(shù)均隨著天頂角的增加不斷增加，這是因為隨著天頂角的增加，傳輸距離逐漸變大，此時星光大氣傳輸會受到更多的湍流影響。除此之外，對比不同觀測高度下的星光閃爍指數(shù)，可以看到不同區(qū)域不同時間段指定天頂角情況下，星光閃爍指數(shù)隨著觀測高度的增加不斷較小，這是因為近地面湍流較強，湍流強度隨高度呈減小趨勢（詳見圖2 中的大氣湍流廓線）。對比同一地區(qū)不同時間段的星光閃爍隨天頂角的變化，可以看出不同時間段星光閃爍指數(shù)之間的差異較小，且對于不同地區(qū)而言，星光閃爍最大值出現(xiàn)的時刻也并不相同（例如：對于三沙地區(qū)而言，0 點和6 點的閃爍指數(shù)要略高于12 點和18 點；對于大柴旦地區(qū)而言，正午時刻星光的閃爍指數(shù)略高于其他時間段）。最后，對比不同地區(qū)的閃爍指數(shù)，可以發(fā)現(xiàn)不同區(qū)域的星光閃爍值之間的差異較小，但相比于沿海三沙地區(qū)，內(nèi)陸地區(qū)（大柴旦、合肥）星光閃爍指數(shù)要略大于三沙地區(qū)。

2.2 到達角和到達角起伏的變化

圖4 和圖5 給出了觀測高度為0 km、1.5 km、3 km、4.5 km、6 km、7.5 km 情況下三沙、大柴旦、合肥一天四個時間段的星光遠場光斑到達角及到達角起伏隨天頂角的變化情況，其中不同形狀的曲線代表不同觀測高度下的結果。

圖4 不同觀測高度下三沙、大柴旦、合肥三地一天四個時間段遠場光斑的到達角隨天頂角的變化情況Fig.4 Arrival angle of the far-field optical spot as a function of the zenith angle at four times of the day in Sansha， Dachaidan and Hefei under different observation altitudes

圖5 不同觀測高度下三沙、大柴旦、合肥三地一天四個時間段遠場光斑的到達角起伏隨天頂角的變化情況Fig.5 Arrival angle fluctuation of the far-field optical spot as a function of the zenith angle at four times of the day in Sansha，Dachaidan and Hefei under different observation altitudes

由圖4、5 中可看出典型區(qū)域，不同時刻以及不同觀測高度的到達角、到達角起伏隨天頂角的變化整體呈現(xiàn)相同的趨勢，到達角和到達角起伏整體與星光閃爍指數(shù)呈正相關關系，到達角主要在μrad 量級變化，到達角起伏主要在10?12量級上變化。從對比圖中的到達角以及到達角起伏的變化曲線，可以發(fā)現(xiàn)：典型區(qū)域，不同時刻以及不同觀測高度下的到達角、到達角起伏整體隨天頂角的增大而逐漸增大。隨著觀測高度的增加，這種增加逐漸減弱，特別是當觀測高度大于3 km 時，到達角、到達角起伏幾乎不隨天頂角的增加而變大，這是因為隨著觀測高度的增加，高空湍流相對較弱，這時改變觀測仰角，星光大氣傳輸并不會受到更多的湍流影響。另一方面，對于特定天頂角情況下，不同區(qū)域四個時間段的到達角和到達角起伏隨著觀測高度的增加逐漸減小。這種減小趨勢在觀測高度在0～3 km 變化時表現(xiàn)最為明顯，這是因為隨著觀測高度的增大，該高度層的湍流急劇減小的原因導致的。

3 結論

本文構建了星光大氣湍流傳輸模型，著重研究了大氣湍流對星光成像的影響。基于ERA5 數(shù)據(jù)和光學湍流預報方法得到大氣湍流參數(shù)廓線，采用等Rytov 指數(shù)間隔方法給出了典型地區(qū)不同時刻的星光傳輸鏈路的相位屏分布情況，結合分步式光傳輸模型，計算了典型地區(qū)不同時刻及觀測條件下星光的閃爍指數(shù)，通過對比星光閃爍的理論值，驗證了數(shù)值計算的可靠性。與此同時，本文研究了湍流大氣中星光遠場光斑質心的抖動特性，得到了典型觀測高度及觀測天頂角情況下的恒星到達角及到達角起伏。研究發(fā)現(xiàn)：星光抖動的到達角和到達角起伏整體呈現(xiàn)相同的變化趨勢，并且與星光的閃爍指數(shù)呈正相關關系。當觀測高度小于3 km 時，隨著星敏感器的觀測高度增加及其觀測天頂角的減小，星光抖動逐漸減小，當觀測高度大于3 km 時，由于此時大氣湍流強度較弱，星光抖動幾乎不受觀測高度、天頂角的變化而變化。綜上所述，提高星敏感器的觀測高度和觀測仰角能一定程度上減輕星光成像的大氣湍流影響。

附錄

圖A1 等Rytov 間隔的相位屏高度分布的流程Fig.A1 Flow chart of the phase screen height distribution with equal Rytov intervals