微觀能量分析氣動光學效應時變誤差的方法

楊博，于賀，*，樊子辰

1.北京航空航天大學 宇航學院，北京 100083

2.北京控制與電子技術研究所，北京 100038

高超聲速飛行器是人類今后發展的戰略運輸工具［1-2］。隨著軍事需求的快速推動，自主導航技術已經成為高超聲速飛行器發展的關鍵問題之一［3］。與其他導航系統相比，天文導航系統具有絕對的自主優勢，同時具備強抗干擾能力、高可靠性［4-5］。然而高超聲速飛行器在大氣層內飛行時，光學頭罩與來流之間形成的復雜流場會產生強烈的氣動光學效應［6-7］，由此阻礙了天文導航在高超聲速飛行器上的順利應用，因此針對高超聲速條件下氣動光學效應的研究具有重要意義。

氣動光學效應問題涉及高超聲速流場、受擾光場2 個復雜系統的耦合。光學窗口上方的復雜流場對傳輸光線產生嚴重的干擾，導致光學探測系統的目標圖像發生偏移、模糊和抖動等畸變現象［8-10］。當前氣動光學研究的關注點主要集中于2 個方向。一種是基于風洞試驗和圖像分析理論來研究目標圖像的畸變特性。通過分析不同風洞試驗條件下獲取的目標圖像研究受擾畸變圖像的移動、模糊程度，進而總結畸變規律以便對光學導航系統的測量偏差進行補償。Sutton［11］將湍流統計量和氣動光學畸變統計量進行了關聯。WyCkham［12］在此基礎上發展了邊界層光程差均方根的標度律。Jumper 和Gordeyev［13］通過統計實驗數據總結了邊界層氣動光學效應的預測模型。這種方法能夠在受擾圖像的基礎上直接獲取較為精確的畸變結果，但是難以揭示光線在復雜流場中受擾的本質過程。另一種研究方向是通過流場觀測技術，如基于背景紋影技術的波前傳感器（BOS-based Wavefront Sensor，BOSWS）［14］、納米粒子示蹤技術（NPLS）［15］，或者通過計算流體力學技術，如大渦模擬法（Large Eddy Simulation，LES）［16］、直接數值模擬法（Direct Numerical Simulation，DNS）［17］等，獲得高速流場的密度數據。當光線通過高超聲速飛行器光學窗口周邊的高速流場時，一般取其湍流邊界層（Turbulent Boundary Layer，TBL）附近區域的折射率場進行分析，利用幾何光學的光線追跡方法研究光線通過高速流場時氣動光學效應引起的相位畸變，用光程差（Optical Path Difference，OPD）、斯特列爾比（Strehl Ratio，SR）來描述受擾圖像的畸變結果［18-19］，這也是目前對氣動光學效應研究最常用的方法。這種結合流場、光路計算的研究過程能夠對光線傳輸的偏折特性進行宏觀層面上的分析，但是無法描述光線在湍流傳輸中的散射、能量耗散過程。采用波動光學對氣動光學效應進行描述時，雖然從理論上可以通過波動振幅反映能量的變化，但在針對湍流內部的光線傳輸時，一般會忽略振幅變化，只利用簡化的麥克斯韋方程組進行解算。通過弱散射理論可以證明在小尺度渦結構中，散射場功率的數值結果不能作為小值處理［20］。因此尋找一種能夠反映光在湍流中的能量傳輸過程的氣動光學效應分析方法，對于探究氣動光學效應的作用本質和完整描述具有重要意義。

光子理論是基于普朗克、愛因斯坦和玻爾等對于光的量子化概念發展出的描述光與物質相互作用的理論。根據能級躍遷假設，具有確定內部運動定態的原子集合中的每個定態的特征具有一定的總能量。當原子發射或吸收光子時，是從能量（EH）較高的定態跳到能量（EL）較低的定態（EH＞EL），其發射一個頻率為ν的光子；反之，則是吸收一個頻率為ν的光子。當處于基態的原子與能量較高的粒子碰撞，粒子將其部分能量轉移給原子去跳到較高能態，粒子則衰變至較低能級，隨之產生光子發射。在分析光的能量傳輸分析上，已有不少學者采用光子傳輸理論對介質產生的光學畸變進行分析［21-22］。而光在高速流場中的傳輸在本質上即是光子與氣流分子的相互作用過程。氣動光學效應則可認為是光子與不均勻流氣體分子相互作用后產生不同能量分布的表征。因此，本文以新的視角，提出基于光子傳輸理論從微觀能量分析的角度剖析光在高速流場中傳輸導致的氣動光學效應時變誤差的方法。通過光子統計估計氣動光學效應中的光子能量損失，建立微觀參量與宏觀特征描述的映射關系，并通過數值仿真和風洞試驗進行驗證，同時得到氣動光學效應的時變誤差，用于高超聲速飛行器天文導航系統的誤差補償。

1 氣動光學效應的微觀機制描述

氣動光學效應的微觀本質是光子在湍流中與氣體分子的相互作用，如圖1所示。以微觀物理學為基礎，研究湍流中氣體分子的吸收、散射特性，以及光子在湍流中的傳輸模型、氣動光學效應在光子能量分析層面的微觀物理表征。

圖1 氣動光學效應微觀角度的光子傳輸過程Fig.1 Photon transmission process of aero-optical effects from microscopic angle

1.1 光子與湍流氣體分子的相互作用

1.1.1 湍流氣體分子對光子的吸收作用

光子經過湍流氣體分子層傳輸時會被氣體分子吸收。令光子經過ds距離傳輸，被氣體分子吸收的概率為μads，定義物理量吸收系數為μa(r，ν，t)，描述吸收過程，是空間、時間及頻率的函數。光子與氣體分子的吸收作用伴隨著電子的極化、電子能態的改變。當光子被吸收后，氣體分子會從基態能級躍遷到激發態能級，如圖2所示。

圖2 光子被吸收產生的能級躍遷Fig.2 Energy level transition caused by photon absorption

不同的躍遷方式取決于吸收光子后的電子自由狀態，由于直接計算單一分子不同躍遷方式的吸收系數十分復雜，因此選擇復折射率的形式對光子吸收系數進行描述。根據洛倫茲的色散理論，復數折射率n?的計算公式［23］為

式中：nR為實部，即幾何光學中常用的介質折射率，用于反映湍流氣體分子的色散特性；nI為虛部，用于反映分子的吸收特性。nR、nI的關系滿足

式中：KGD為Gladstone-Dale 常數，其取值由光束的波長決定，KGD取為

式中：λ為入射光的波長。

根據湍流分子對光子的吸收系數μa與復折射率虛部的關系滿足

因此可以發現，吸收系數與空間、時間、頻率均相關。由于湍流氣體分子的分布具有時變性、不均勻性，因此不同時間光子經過不同區域的湍流氣體分子傳輸時，由于分子間距變化的不穩定，光子被氣體分子吸收作用不同，從而會造成光子能量損失程度的差異。

1.1.2 湍流氣體分子對光子的散射作用

光子經過湍流氣體分子傳輸時，氣體分子紊亂度的增加導致了光子散射作用增強，也會引起光子能量擴散及光強減弱。

假設光子在氣體中經過ds距離后被散射的概率為μs(ν)ds，定義散射系數μs(r，ν，t)為描述該過程的物理量。光子的散射與湍流分子不均勻性的尺度有關，本文主要考慮小尺度湍流氣體分子的瑞利散射，其由分子密度變化和偶極矩變化形成的光學不均勻性所致。根據光子散射理論，單分子的瑞利散射截面為

湍流分子的散射截面與λ4成反比，為了減少單一分子散射截面的計算量，實際散射系數可以表示為

通過式（8）可以看出湍流分子的散射系數也是空間、時間、頻率的函數。對于高速湍流，氣體分子轉捩程度的增加將導致散射的增強、光子能量的擴散。

1.2 光子在湍流氣體中傳輸模型

光子在湍流氣體中傳輸時與氣體分子相互作用，光子數目發生變化。設在任意時刻t，處于位置坐標r處、頻率ν的光子，沿著Ω方向傳播時，定義光子的分布函數可表示為f(r，ν，Ω，t)，單元光子數dP可表示為

式（9）即表示t時刻在位置點r附近的體積單元dV內，通過頻率處于ν附近dν單元內、方向處于Ω附近dΩ單元內的光子數。光子的分布函數f(r，ν，Ω，t)決定了湍流中光子場的強度I，關系公式為

式中：c為光速；h為普朗克常數；I(r，ν，Ω，t)為t時刻處于位置點r、在Ω方向上傳輸的頻率為ν的光子的能量，推出光子在湍流中傳輸的玻爾茲曼方程［24-25］為

式中：υ為光子速度矢量；a為加速度矢量；若忽略相對論效應，則a=0，且υ=cΩ；?r、?υ分別為幾何空間、速度空間的散度算子；Q(r，Ω，t)為光子源項，表征了光子與湍流氣體分子相互作用的過程，本文主要研究與光子能量相關的吸收、散射效應。將直角坐標系x-y-z平移獲得光子局部坐標系ex-ey-ez，單位方向向量余弦和光子方向矢量在坐標軸上的投影關系如圖3所示。

圖3 光子局部坐標系的矢量關系Fig.3 Vector relation of photon local coordinate system

令單位方向矢量Ω=[μ η ξ]T，在直角坐標系中滿足

根據式（10）、式（11）可以得到光子在湍流中的傳輸模型［26-27］為

式（13）等號右側2 項分別為光子源項Q(r，Ω，t)吸收、散射部分的表述，滿足

將式（13）改寫為散度形式［26-27］，即

其中：

式中：i、j、k分別為x、y、z坐標方向的單位矢量。光子傳輸方程不僅包含空間微分項，還包含角度積分項。針對某一湍流狀態的光子傳輸過程，采用離散坐標法對微分方程進行求解，因此，需要分別對空間、角度進行離散化。首先，將空間計算域劃分為體積單元，并通過角度分段常數求積分的方法將立體角劃分為多個部分，如將天頂角θ、圓周角φ離散為Nθ×Nφ份，因此，可以將光子傳輸方程寫為離散坐標形式，即

式中：離散的天頂角、圓周角為

對天頂角和圓周角的導數需要進行有限差分近似，可以將角分配項寫為

因此光子在湍流分子中傳輸的離散坐標方程為

在體積單元ΔVP上對式（21）進行積分，可以得到離散坐標形式

其中：

式中：NP為單元P所具有的表面數；Ai為第i個表面的面積；ei、nR，i分別為第i個表面單位外法向量、折射率的實部；ΔVP為單元P的體積；nR，P為單元P中心的折射率實部。為了使得積分方程封閉，需要設置邊界條件，這里選擇下游表面的光子強度值等于上游中心的強度值，滿足

式中：Pi為與單元P共享第i個表面的相鄰單元；IPi(r，ν，Ωm，n)為單位Pi的光子強度。上述過程是一個方向上的空間坐標離散化，其他方向相同，這樣就形成了一個離散方程組。光子強度通過求解和迭代離散方程來求解，詳細過程見文獻［26-27］。

綜上可以看出，光子在湍流中的傳輸方程為一個微分方程，變量為空間、時間坐標，由于計算時間變量會導致計算量急劇增長，且光子運動速度遠超過流場的變化速度，可以認為在一個較短的時間間隔內，光子與湍流分子作用期間不會影響湍流流動狀態，湍流結構是相對靜態的，但是為了分析氣動光學效應的時變誤差，因此后續的分析中，需要限制流場的采樣時間。

2 氣動光學效應的能量分析評估函數

湍流氣體分子運動的不確定性、分布的不均勻性，使得在其中傳輸的光子能量產生耗散現象，最終導致實際光子在接收器上的落點（星點）與成像平面中心的標準落點之間產生偏差。事實上，在成像平面上獲得的“星點”即是光子能量統計的結果，并且僅是以能夠出現在接收面上的光子的統計結果。因此，定義微觀意義下描述光子能量耗散程度的評價函數是必要的。

2.1 光子的能量散度

光子的能量散度（Photon Energy Divergence，PED）用于描述光子散射作用程度。以某一傳輸過程的接收平面為評估對象，統計矢徑在某一確定接收平面Λ上的光子分布情況。假設傳輸過程中散射并不改變光子的頻率，那么在平面Λ處的PED 為該平面上光子分布的集中度，計算公式為

式中：r0為光子直線傳輸情況下在接收平面上的矢徑。PED 為微觀統計意義下的接收平面上各處離散光子能量分布的集中度，PED 的值越大，表示光子傳輸的散度越大，光子傳輸過程中遇到的散射作用越強。

2.2 光子的能量耗散比

光子的能量耗散比（Photon Energy Loss Ratio， PELR）為光子被吸收程度的物理描述。取傳輸過程的某一接收平面Λ，PELR 為該平面上的光子總能量與初始光子能量的損耗比例，即

式中：I0為初始光源的強度。PELR 的值越大，表示光子傳輸過程中能量的耗散比例越大，由此可以說明湍流分子對光子的吸收作用。PELR 可以描述相對初始光源的光子能量的損耗。

2.3 光子的偏轉角

光子的偏轉角（Photon Deflection Angle，PDA）用于描述光子傳輸過程中速度方向偏轉程度。對于光子傳輸過程中的某一接收平面Λ，光子的偏轉角定義為，該平面上任意矢徑處光子的即時速度方向與初始傳輸方向的夾角：

式中：Ω0為初始光子的傳輸方向矢量。為了反映氣動光學畸變的宏觀特征，定義平均光子偏轉角為平面Λ上任意矢徑處光子數權重與其偏轉角乘積的和，即

PDA 的值越大，說明光子在傳輸方向上的偏移畸變越大。相對于傳統幾何光學的相位畸變分析，PDA 可以更加直接地描述光子傳輸方向受湍流結構影響的程度。

3 仿真分析與實驗數據對比

3.1 高速湍流流場的獲取

采用大渦模擬法（LES）獲取高速流場的瞬態結構［28］。高速飛行器的物理結構及尺寸如圖4所示。其中，光學窗口為凹窗設計，錐頭斜角為13.5°。xb-yb-zb、xw-yw-zw分別為本體、窗口坐標系。

圖4 高速飛行器的物理結構及尺寸Fig.4 Physical structure and size of high-speed aircraft

飛行器與光學窗口附近網格劃分如圖5所示，為了更加清晰地得到光學窗口附近的流場結構，對窗口附近網格進行加密處理。

圖5 飛行器與光學窗口附近網格分布Fig.5 Grid distribution near aircraft and optical window

飛行器的網格劃分數量為1.895×107，且已經過網格無關性驗證。LES 的仿真采用商業軟件Fluent 2021R2 及北京超算中心的高性能服務器進行湍流模擬。其中LES 的仿真條件為：飛行高度取20 km；飛行速度為Ma=3.0；來流密度為0.089 kg/m3；飛行攻角為設為0°；仿真時間步長為10-7s。流場模擬結果如圖6、圖7所示。

圖6 飛行器外流場的速度分布Fig.6 Velocity distribution of aircraft outflow field

圖7 光學窗口上方流場的密度分布Fig.7 Density distribution of flow field above noptical window

3.2 湍流渦結構對光子的能量耗散分析

從微觀分子角度看，湍流結構尺度實際上即為湍流分子的密度、速度分布不同所致，剛進入光學窗口的氣流分子由于窗口壁面黏性相對大導致分子速度變慢，較窗口之外氣體的流速變化率增大，分子與分子間的相對運動加劇，能量傳遞范圍加大，在窗口前部形成了湍流大尺度渦，而經過一段距離傳輸之后，分子動能逐漸耗散，最終湍流大尺度渦轉捩為小尺度渦［29］。這一變化過程中，湍流各種尺度渦結構對當前接受到的光束的光子的散射作用不同，接收窗口對經過湍流傳輸的光子的耗散能量也不同。采用Q準則確定湍流流場數據內部渦結構，即通過速度的旋度來進行渦結構的判定，計算公式為

式中：Vx、Vy、Vz分別為氣流在x、y、z方向上的速度分量。在飛行器高速飛行時，在高速湍流中存在不同尺度的時變渦結構，不妨將其分為大、中、小尺度結構，如圖8所示。

圖8 Q 準則下的不同尺度湍流結構Fig.8 Turbulence structures of different scales under Q criterion

湍流的渦結構與氣流分子的速度和密度分布均相關，為了探究渦結構對光子傳輸的干擾規律，分別在大、中、小尺度渦結構數據情況下，使用上述光子傳輸的能量分布來分析渦結構尺度對光子能量耗散的影響。其中渦結構的分離采用本征正交分解（Proper Orthogonal Decomposition，POD）來進行，將其能量模態按照降序分為3種進行重構［30-32］。假設光源初始設置為在邊界層流場中心高斯分布的點光源，入射角為90°，波長為572 nm，光子數目為108，光子接收成像面尺寸為50 mm×50 mm 的方形尺寸。光子在3 種尺度流場結構中傳輸模擬得到的光子能量散布，如圖9～圖11所示，其中紅色方框為光子傳輸達到接收器成像平面的最小“能量圈”。

圖9 大尺度湍流結構示意圖與光子能量散布圖Fig.9 Schematic diagram of large-scale turbulence structure and photon energy distribution

圖10 中尺度湍流結構示意圖與光子能量散布圖Fig.10 Schematic diagram of meso-scale turbulence structure and photon energy distribution

圖11 小尺度湍流結構示意圖與光子能量散布圖Fig.11 Schematic diagram of small-scale turbulence structure and photon energy distribution

通過圖9～圖11 可以看出，假定同樣數量的光子通過大尺度結構之后，可以有較多的光子達到了成像平面上，即紅框之內，并且光子集中度高，亮度（光強）大，向著成像面中心集中。同理看到通過中尺度結構次之，而通過小尺度渦結構后，散射到紅框之外的光子增加（見圖11），并且紅框內即達到成像平面的光子減少，且集中度開始下降，亮度（光強）減弱，從以上分析可看出，流場中小尺度結構是對光子能量耗散起主要作用，可造成更多的光子散射，其宏觀表象就是“星點”的模糊和抖動。

光子在湍流分子中傳輸過程，偏轉角達到一定界限時，光子便無法進入接收器的成像平面內，通過對不同尺度的渦結構影響的光子統計可以獲得散射作用與渦結構之間的聯系。表1 為根據能量分析的評價函數分別對圖9～圖11 的3 種情況進行氣動光學效應評估的結果。

表1 不同尺度渦結構下光子評價函數值Table 1 Photon evaluation function value with different scale vortex structures

由表1 可以看出光子能量散度PED 隨著渦結構尺度減小而增大，而光子能量耗散比PELR、光子傳輸偏轉角PDA 也是隨著渦結構尺度減小而增大，評價函數值的變化規律與圖9～圖11 不同渦結構尺度的光子能量散布規律相吻合。

3.3 光子的微觀能量分析與宏觀表征的統一

基于光子能量層面研究氣動光學效應是基于光子傳輸統計方法對光子能量耗散的定量分析。通過微觀能量分析的量化結果與傳統幾何光學的宏觀OPD、點擴散函數（Point Spread Function PSF）、宏觀偏角（Light Offset Angle，LOA）3 個定性分析的物理量進行比較，驗證光子微觀能量分析方法的準確性。

1） 微觀統計光子光強分布和宏觀PSF 對比

通過3.1 節中獲取的高速流場進行光子傳輸仿真，取距離光學窗口不同高度的采樣平面進行解析，通過I(r，ν，Ω，t)=chνf(r，ν，Ω，t)將光子能量分布轉化為光子場的強度。并用幾何光學的“光線追跡”法得到相應采樣平面的PSF 對應的光強，圖12 為不同傳輸距離上兩者的分布表示。

圖12 微觀統計光子場強度與平均PSF 的對比Fig.12 Comparison of micro-statistical photon field intensity and average PSF

對比圖12 中2 條曲線，接近重合，但仍然有平均1.2%的相對誤差。這里有2 點需要說明：其一，圖12 僅是對落到截取平面上的光子進行了光子微觀統計，并不包括由于湍流影響產生了能量耗散未被接收到的光子（如同圖9～圖11 中表示的接收能量邊界紅框外的光子）；其二，圖12 中PSF 是用幾何光學“光線追跡”得到的采樣點點擴散函數平均值分布，即各采樣點光場上“像點”的平均光強的分布，因此，兩者之間存在誤差是必然的。圖13、圖14 分別為微觀統計光子場強度分布和傳統宏觀PSF 的光強分布，結果也是明顯相近的，說明微觀能量分析的方法是準確的。

圖13 微觀統計光子場強度分布Fig.13 Micro statistical photon field intensity distribution

圖14 傳統宏觀PSF 的光強分布Fig.14 Light intensity distribution of traditional macro PSF

2） 光子統計光程差PS-OPD 與宏觀OPD對比

類比傳統OPD 的物理表述，在某一確定接收平面Λ中的某一矢徑r處，設當前位置光子數為NP=f(r，ν，t)，則定義光子統計光程差（Photon Statistical Optical Path Difference， PS-OPD）為

式中：r0為光子不受干擾傳輸情況下在接收平面上的矢徑。

由于瞬態流場更能體現湍流流場的動態高頻特性，因此取流場中光子傳輸距離（與圖12 采用相同的采樣間隔距離）上的不同采樣平面，統計得到各采樣平面的平均PS-OPD、通過幾何光學的“光線追跡”對各采樣平面求得的宏觀平均OPD，如圖15所示。對比兩曲線看，2 種光程差仍然存在一定的偏差，產生誤差的原因也是由于部分光子未達到接收表面而影響了統計結果，同時可以從圖16、圖17 中比較傳統幾何光學宏觀OPD 與微觀光子統計光程差PS-OPD，可以看出來其趨勢基本一致，也說明了微觀光子能量分析的方法與宏觀方法的統一性。

圖15 光子統計平均PS-OPD 與宏觀平均OPD 對比Fig.15 Comparison between photon statistical average PS-OPD and macroscopic average OPD

圖16 傳統幾何光學宏觀OPD 分布Fig.16 OPD distribution of traditional geometric optics

圖17 微觀光子統計光程差PS-OPD 分布Fig.17 PS-OPD distribution of microscopic photon statistical optical path difference

圖18 微觀光子統計平均PDA 與宏觀LOA 對比Fig.18 Comparison between microscopic photon statistical average PDA and macroscopic LOA

3.4 時變誤差仿真和微觀能量分析方法驗證

隨著高速流場的不斷變化，氣動光學效應誤差也在不斷發生變化，實際仿真過程中，由于光子與湍流分子相互作用的時間比湍流變化更短，因此采用流場的不同采樣時間間隔作為分析氣動光學效應時變誤差的步長。并將最終的仿真結果與風洞試驗數據進行對比，驗證微觀能量分析方法的有效性。

在LES 仿真40 000 代之后，飛行器光學窗口上方流場已經完全發展為湍流狀態。此時取20 ms 的流場變化數據，采樣200 個流場狀態，經過光子傳輸仿真獲得--- -----PDA 時變誤差結果如圖19所示。

圖19 基于能量分析的氣動光學效應的時變誤差Fig.19 Time-varying error of aero-optical effect based on energy analysis

為了進一步驗證本文中的微觀能量分析方法，使用文獻［33-34］中的CFD 模擬和風洞試驗結果進行比較，驗證結果如圖20所示，其中，ξ為入射角，滿足ξ∈(0，180)°，光子能量分析方法的OPDrms實際代表微觀統計PS-OPD 的均方根，這里為了表示方便，對氣動光學效應的誤差均采用OPDrms進行描述；OPDrms(90°)則表示ξ=90°時的OPDrms。最終從驗證結果可以發現氣動光學效應的微觀能量分析方法結果與實際風洞試驗的結果基本吻合。

圖20 不同分析方法的結果比較Fig.20 Comparison of results from different analysis methods

4 結 論

針對當前高超聲速飛行器光學探測中氣動光學效應的瓶頸問題，創新性地提出了基于光子微觀能量分析氣動光學效應的思想，并進行了詳細的推導和驗證，得出以下主要結論：

1） 本文推導了光子在湍流中的傳輸模型，并基于光子與湍流分子的作用機制，建立了基于光子能量分析氣動光學效應時變誤差的方法。

2） 分析了湍流不同尺度結構與光子能量耗散的關系，得出湍流小尺度結構是造成光子能量耗散的主要原因，揭示了湍流影響下的光接收成像平面上的“星點”僅是部分具備有“可達到”能量的光子散布，小尺度渦結構可使散布范圍加大，即表現為模糊、抖動加劇。

3） 從微觀光子角度建立了氣動光學評價函數，并驗證了微觀統計的量化值與宏觀定性描述物理量之間映射關系的合理性和準確性。

4）進行了氣動光學效應時變誤差仿真，可以從光子統計的角度進行氣動光學畸變規律的分析。

本文主要提出基于光子能量分析氣動光學效應的微觀思想，為高超聲速飛行器氣動光學效應的研究提供新的研究思路，后續還需針對飛行環境的適應能力進行進一步評估。