一種高精度氣動光學星光偏折預測方法

劉宇鑫, 王新龍*, 丁 偉, 胡曉東

(1. 北京航空航天大學 宇航學院, 北京 100083; 2. 北京控制與電子技術研究所, 北京 100038; 3. 航空工業 西安飛行自動控制研究所, 西安 710065)

0 引 言

高超聲速飛行器具有飛行速度快、 突防能力強等優點, 是當前世界各國軍事博弈的重要戰略武器[1-2]。 受其任務需求影響, 高超聲速飛行器必須具備更強的導航自主性和可靠性[3]。 基于星光折射的新型天文導航方法有自主性強、 導航精度高的特點, 因此將天文導航技術拓展到高超聲速環境中進行研究, 對高超聲速飛行器實現全自主導航有重要的價值和意義。

然而, 飛行器在大氣層內高速飛行時與來流發生強烈的相互作用, 會在機體周圍形成激波層、 邊界層、 剪切層等復雜的非均勻高速氣體流場[4]。 星光穿過復雜非均勻流場時會發生氣動光學效應, 導致其傳播方向發生偏折, 引起星光矢量測量誤差。 當飛行器以馬赫數3以上的速度在臨近空間飛行時, 氣動光學效應會使觀測到的星光矢量發生幾至十幾角秒的偏折[5-6], 嚴重影響了星敏感器的測量精度, 因此對星光偏折進行預測校正的意義重大。

目前對氣動光學效應星光偏折的預測方法可分為兩類: 理論模型預測法與數據擬合模型預測法。 理論模型預測法是基于空氣動力學理論, 建立高超聲速流場的解析模型, 進而基于折射定律建立理論模型用于預測[7-8]。 但其建模過程中僅考慮了高超聲速飛行器上楔面激波所導致的光線偏折, 忽略了激波后非均勻流場對光線的偏折作用, 因此其預測精度有限。 由于星光偏折角與飛行高度、 飛行馬赫數、 攻角以及視線角等因素間具有復雜的非線性映射關系[9], 數據擬合模型預測法是基于實驗或計算得到的氣動光學偏折數據, 利用神經網絡、 支持向量機等智能優化算法[9-12]擬合得到相應的光線偏折角模型, 進而實現對偏折角的實時在線預測。 但由于偏折角與其影響因素間的映射關系復雜, 具有極強的非線性特點, 這種方法需利用光線追跡法計算大量數據進行擬合, 預測精度受訓練樣本的影響大。 此外高超聲速流場數據的網格量大且形狀不規則, 光線追跡法需耗費大量的計算資源, 制約了其應用的便捷性。

本文通過分析星光在高超聲速飛行器外流場中發生偏折的機理, 綜合考慮激波以及激波后流場的影響, 建立了一種精確的星光偏折角模型, 進而基于星光偏折角模型設計了一種高精度的氣動光學星光偏折預測方法。

1 氣動光學星光偏折角模型

1.1 星光在非均勻流場中的傳播路徑分析

當超聲速氣流流過飛行器表面時, 會形成復雜的非均勻流場結構, 根據流場中氣體的流動情況, 可將流場劃分為如圖1所示的激波與激波后流場, 進而分別分析其對星光的偏折作用。

圖1 星光穿過非均勻流場發生偏折Fig.1 Starlight deflects through an inhomogeneous flow field

圖1中的激波是由諸多微弱壓縮波堆疊在一起形成的突躍壓縮波, 是高度壓縮的流場形態, 會在連續介質中產生一個幾乎無厚度的、 狀態劇烈變化的曲面, 導致激波前后氣體的密度等物理參數出現“斷面”式顯著變化。 超聲速氣流經過激波以后, 其流向折轉一個較小角度, 然后氣流在激波后連續地進行等熵壓縮[13], 導致激波后流場氣體的密度等物理參數連續變化。

由洛倫茲—洛倫茨公式可知, 氣體折射率n與密度ρ間具有如下的線性關系[14]:

n=1+KGD·ρ

(1)

式中:KGD為格拉斯通-戴爾常數, 對于恒星星光所處的可見光波段, 可將其視為常值。

氣體密度分布不均會導致其折射率分布不均。 星光穿過折射率分布不均的流場時會發生折射, 導致其傳播方向發生偏折。 如圖1所示, 星光在穿過激波面時, 由于流場氣體折射率發生突變, 會發生明顯的折射現象, 星光的傳播方向產生較大幅度的偏折; 在穿過激波后流場時, 由于流場氣體折射率連續變化, 星光會發生持續的折射, 其傳播方向發生連續的偏折。 因此, 可根據星光在流場中的傳播過程, 分別建立激波以及激波后流場導致的星光偏折角模型。

1.2 激波導致的星光偏折角模型

在星光傳播過程中, 當其以入射角η0射入激波時, 發生偏折的示意圖如圖2所示。

圖2 星光在激波處發生偏折Fig.2 Starlight deflects at the shock wave

根據Snell定律[14]可計算得到星光在激波界面處的折射角ηd1:

(2)

式中:n1、n2分別為激波前后流場氣體的折射率。

由此, 可得激波所導致的星光偏折角Δηs為

(3)

將式(1)的洛倫茲—洛倫茨公式代入式(3), 可得

(4)

由式(4)可知, 激波導致的星光偏折角受激波前后的氣體密度ρ1、ρ2與星光入射角η0的影響。

激波前氣體的密度ρ1為來流氣體的密度, 與高度h近似成指數關系:

(5)

式中:ρ0為高度h0處的大氣密度;H為密度標尺高度。

而激波后的氣體密度ρ2又可根據激波關系式計算得到[13]:

(6)

式中:γ為氣體的比熱比, 對于空氣取γ為常值1.4;M1為來流馬赫數;β為激波角, 即激波與氣流速度間的夾角, 對于錐形體而言, 激波角可通過求解Taylor-Maccoll方程計算得到[13], 計算所得激波角β與半錐角δ和馬赫數M1的關系曲線如圖3所示。

圖3 錐面激波的半錐角、 激波角及馬赫數關系曲線Fig.3 Relationship between half cone angle, shock angle and mach number of conical shock wave

將式(5)～(6)代入式(4)的偏折角計算公式中, 可得激波導致的星光偏折角模型為

(7)

由式(7)可知, 激波導致的星光偏折角受飛行高度h, 飛行馬赫數M1, 激波角β以及星光入射角η0的影響。

1.3 激波后流場導致的星光偏折角模型

超聲速氣流經過激波以后, 在激波后流場中連續等熵壓縮[13], 使得激波后流場的折射率連續變化, 進而引起星光的連續偏折, 如圖4所示。

圖4 星光在激波后流場發生偏折Fig.4 Starlight deflects in the flow field behind shock wave

根據錐形流理論, 繞軸線旋轉而成的任意一個中間錐面上, 所有氣流參數都均勻分布[13], 因此可以Δθ為角距離將激波與壁面間的流場劃分為m個子層, 使得折射率在每個子層內保持恒定。

星光在激波處發生折射后以ηd1角度射出, 遇到第1個子層。 在該界面處的入射角為ηd1i, 折射角為ηd2, 在該界面應用Snell折射定律可得

(8)

由幾何關系可知, 星光在第2個子層的入射角ηd2i為

(9)

在連續的中間界面處應用Snell定律, 可得星光在第m個子層發生折射后最終的出射角η1為

(10)

由此可得激波后流場所導致的偏折角Δηf為

Δηf=η1-ηd1-β=η1-ηd(m-1)i+ηd(m-1)-

ηd(m-2)i+…+ηd2-ηd1i=

(11)

將式(1)代入式(11), 可得

(12)

由式(12)可知, 激波后流場導致的星光偏折角受各激波子層密度ρd(j)的影響。

各激波子層密度可通過對式(13)的球坐標系下的Taylor-Maccoll方程組進行數值積分, 并結合氣體的等熵壓縮方程計算得到[13]:

(13)

式中:vr為徑向速度;vθ為軸向速度。

該方程組的求解可以使用Runge-Kutta方法, 以式(14)所示的緊靠激波后的流場速度作為初始條件進行積分求解:

(14)

式中:T1為來流氣體溫度;h*為滯止焓;R為理想氣體常數。

各子層的氣體密度受飛行高度h, 飛行馬赫數M1, 激波角β以及子層位置的影響, 而激波角又可根據飛行馬赫數, 飛行攻角以及半錐角計算得到, 因此可將各子層的氣體密度ρd(j)表示為飛行高度h, 飛行馬赫數M1, 飛行攻角α, 飛行器半錐角δ以及子層數j的函數, 即

ρd(j)=f(h,M1,α,δ,j)

(15)

將式(15)代入式(12)的偏折角計算公式, 可得激波后流場導致的星光偏折角模型為

(16)

綜合考慮式(7)的激波導致的星光偏折角模型及式(16)的激波后流場導致的星光偏折角模型, 可得氣動光學效應所導致的星光偏折角為

Δη=Δηs+Δηf

(17)

式(17)推導的為星光矢量與折射點法線所構成平面內的偏折角, 為便于對星光偏折進行補償, 可以星敏感器坐標系(s)為基準, 將星光偏折角分解為沿Xs和Ys方向的分量Δηx和Δηy。 定義如圖5所示的星光入射角。

圖5 星光入射角示意圖Fig.5 Diagram of starlight incidence angle

圖中星光矢量傾角ψ為星光矢量與XsOYs平面間的夾角, 星光矢量滾轉角φ為星光矢量在XsOYs平面內的投影與Xs軸間的夾角。 根據幾何關系可得星光偏折角在星敏感器坐標系下Xs和Ys方向的分量Δηx、 Δηy為

(18)

將式(17)代入式(18)后, 可得氣動光學效應所導致的星光偏折角模型, 進一步可基于該模型對氣動光學效應星光偏折角進行實時預測。

2 基于星光偏折角模型的偏折預測方法

通過建立的星光偏折角模型式(18)可知, 為計算星光偏折角, 需已知飛行器的半錐角δ、 飛行高度h、 飛行馬赫數M1、 攻角α, 激波后流場的密度分布ρd(j)以及星光矢量入射角ψ、φ。 其中飛行器的半錐角δ在設計時已知; 飛行高度h、 飛行馬赫數M1、 攻角α等飛行參數可由嵌入式大氣數據傳感系統(FADS)實時測量所得; 星光矢量入射角可由星敏感器測量所得, 而激波后流場的密度需通過對式(13)的Taylor-Maccoll方程組進行數值積分, 并結合氣體的等熵壓縮方程計算得到。 但Taylor-Maccoll方程組的形式復雜, 對其求解的計算量大, 難以用于氣動光學星光偏折角的實時在線預測。 因此, 可預先離線建立如式(15)所示的流場密度代理模型, 用于星光偏折角的實時在線預測。

2.1 密度代理模型

激波后流場各點氣體密度的分布規律與飛行參數相關, 因此可利用神經網絡能夠充分逼近任意復雜的非線性關系的特點, 基于流場計算數據對飛行參數與激波后密度的分布關系進行擬合, 建立其密度代理模型。

BP神經網絡是目前應用最廣泛的網絡, 其模型結構簡單、 運行速度快, 適合于求解內部機制復雜的問題, 適用于本文的建模環境。 BP神經網絡進行訓練時, 通常使用梯度下降法來最小化損失函數, 進而更新網絡的閾值和權值。 但其全局搜索能力一般, 易使網絡權值與閾值困于局部最優值, 進而影響預測效果。

因此, 采用全局搜索能力強, 尋優速度快的粒子群算法(PSO)對經典的BP網絡進行優化, 進而利用PSO-BP網絡建立密度代理模型。

2.1.1 粒子群優化算法

粒子群優化算法是一種群體智能算法, 該算法參考了鳥群覓食的自然現象, 通過模擬鳥群間的協同合作和知識共享, 利用自身的經驗和同伴中最好的經驗來決定下一步的運動, 使粒子群以最短的時間尋得全局最優解[15-16]。

粒子群中各粒子依據式(19)與式(20)對自身的速度和位置分別進行更新:

(19)

(20)

2.1.2 基于PSO-BP的密度代理模型

利用PSO對BP網絡進行優化, 將神經網絡各層的連接權值編碼成粒子, 適應度值為使用該組權值時的網絡輸出均方誤差。 利用PSO算法, 在預設的迭代次數內搜索最優的網絡權值, 進而獲得神經網絡中的參數。

以飛行高度h、 飛行馬赫數M1、 攻角α, 以及流場中各錐面與飛行器軸線間的夾角為輸入變量, 以各錐面上的氣體密度為輸出變量, 利用PSO-BP網絡建立其密度代理模型, 流程如圖6所示。

圖6 PSO-BP網絡算法流程圖Fig.6 PSO-BP neural network algorithm flow chart

具體步驟如下:

(1) 確定BP神經網絡的拓撲結構。

(2) 根據BP網絡的拓撲結構, 設定種群粒子數以及粒子維數。 另設置慣性權重w、 學習因子c1和c2, 設定收斂精度和最大迭代次數, 并對每個粒子的位置和速度進行隨機初始化。

(3) 將各粒子作為BP網絡的參數進行訓練, 將訓練誤差作為各粒子的適應度值, 并更新當前的局部最優解和全局最優解。

(4) 根據式(19)與式(20)對空間中的粒子的速度和位置進行不斷的迭代更新。

(5) 基于迭代次數與粒子位置判斷是否滿足迭代終止條件。 若滿足則停止迭代, 輸出全局最優位置, 否則返回步驟(3), 繼續迭代直至滿足預設要求。

(6) 滿足最大迭代次數或達到誤差精度要求后終止尋優算法, 將當前的最優參數賦值給BP神經網絡用于流場氣體密度分布的預測。

2.2 星光偏折角預測方法

這樣, 基于PSO-BP網絡建立的密度代理模型并結合式(18)的星光偏折角模型, 可對星光偏折角進行實時在線預測, 該方法如圖7所示。

圖7 星光偏折預測方法框圖Fig.7 Block diagram of starlight deviation prediction method

在預測過程中, 利用FADS系統測量得到當前的飛行高度、 飛行馬赫數以及攻角, 利用星敏感器測量得到星光矢量在星敏感器坐標系下的入射角。 以飛行參數作為密度代理模型的輸入, 進而通過密度代理模型計算得到當前飛行狀態下流場氣體密度的分布, 并結合所建星光偏折角模型計算得到偏折角的預測值用于后續的校正補償。

3 仿真驗證與分析

以典型的帶凹窗結構的高速飛行器為例, 對所提氣動光學星光偏折預測方法的有效性進行驗證。

3.1 數據集的構建

仿真中計算得到飛行速度馬赫數3～9, 高度從20～50 km、 攻角從-15°～15°變化時, 不同飛行條件下的飛行器外流場。 以馬赫數3、 高度20 km、 攻角0°的飛行條件為例, 其流場密度分布如圖8所示。

圖8 流場密度分布圖Fig.8 Flow field density distribution diagram

光線追跡法是通過求解光線傳播方程, 從而得到光線在非均勻介質中的傳輸路徑的一種數值計算方法。 該方法具有較高的理論精度[6], 因此仿真過程中利用四階龍格庫塔光線追跡法計算星光穿過流場后的偏折角, 并將其作為基準, 對所提星光偏折預測方法的性能進行驗證與分析。

仿真中計算了星光矢量傾角以及滾轉角從0°～70°變化時的星光偏折角, 部分計算結果如表1所示。

表1 星光偏折角部分計算結果Table 1 Calculation results of starlight deviationangle

由計算結果可知, 氣動光學效應導致的星光偏折為幾至十幾角秒量級, 嚴重影響了星敏感器的測量精度。

3.2 密度代理模型預測結果與分析

利用計算所得的196組不同飛行條件下的流場密度分布數據, 對數據進行歸一化處理后, 隨機選取其中85%的數據作為模型訓練樣本集, 剩余數據作為模型測試樣本集。

密度代理模型采用3層BP網絡, 輸入層節點數為4, 隱藏層節點數為7, 輸出層節點數為1, 激活函數采用 sigmoid 函數, 學習速率為 0.01。 粒子群算法中的種群規模為40, 最大迭代次數為300, 學習因子c1=c2=2, 最大更新速度為5。

在利用粒子群算法尋找全局最優的權值與閾值過程中, 粒子群的適應度變化曲線如圖9所示。

圖9 粒子群算法適應度變化曲線Fig.9 Fitness curve of particle swarm optimization

由圖9可知, 粒子群算法在迭代20次左右趨于平穩, 優化得到了全局最優的神經網絡權值與閾值。

以馬赫數5、 高度20 km, 攻角0°飛行工況為例, 由PSO-BP密度代理模型預測所得的密度分布與實際的密度分布如圖10所示。 通過對比預測值與實際值可知, PSO-BP密度代理模型預測所得的密度分布規律與實際分布規律相符, 且貼合程度較高, 預測效果良好。

圖10 PSO-BP密度代理模型預測密度分布Fig.10 Density distribution predicted by PSO-BP density surrogate model

進一步, 對PSO-BP密度代理模型在測試集上的預測精度進行統計, 得到其預測的均方根誤差為0.073 9 kg/m3, 擬合優度R2為0.997 8, 能有效預測不同飛行工況下的飛行器外流場密度分布情況。

3.3 星光偏折角預測結果與分析

以四階龍格庫塔光線追跡法計算所得星光穿過流場后的偏折角為基準, 與所提方法的預測結果進行對比。 以其中100組數據為例, 預測所得的各向偏折角與偏折角的真實值如圖11所示。 通過對比仿真結果可知, 所提方法的預測值與實際值保持一致, 誤差相對較小。

圖11 星光偏折角預測結果Fig.11 Prediction results of starlight deviation angle

為進一步驗證所提方法的性能, 將所提方法的預測誤差與文獻[9]的基于IPSO-BP網絡的數據擬合模型預測法的誤差進行對比, 相應的偏折角預測誤差如圖12所示, 其誤差統計結果如表2所示。

表2 星光偏折角預測精度統計結果Table 2 Statistical results of prediction accuracy

圖12 星光偏折角預測誤差Fig.12 Prediction error of starlight deviation angle

由表2可知, 所提方法能實現亞角秒級別的星光偏折角預測精度, 預測誤差約為0.4″, 能補償80%以上的氣動光學效應所導致的星光偏折。 而雖然IPSO-BP模型預測法利用改進的PSO算法對BP網絡進行優化并用于預測模型的訓練, 具有更強的全局尋優和非線性擬合能力。 但受訓練樣本的精度以及數量限制, 其在擬合星光偏折角與攻角、 馬赫數、 飛行高度、 視線傾角和視線滾轉角間復雜的非線性關系時所能取得的精度有限。 所提預測方法由于綜合利用了星光偏折角模型進行預測, 緩解了傳統數據擬合模型預測法過分依賴訓練數據的問題, 能取得較好的預測效果。 通過對比預測誤差以及擬合優度可知, 所提預測方法的預測精度優于基于IPSO-BP網絡的模型預測法。

此外, 基于IPSO-BP網絡的模型預測法在建立預測模型時, 需利用光線追跡法計算不同入射角下的星光偏折角, 而由于高超聲速飛行器流場數據的網格量大且形狀分布不規則, 進行光線追跡計算需耗費大量計算資源。 與此同時, 由于星光偏折角與其影響因素間具有復雜的非線性關系, 需利用大量數據進行擬合, 在仿真中使用了2 564組數據進行訓練。 而所提方法在建立星光偏折預測模型過程中, 利用所建立的偏折角解析模型替代了利用光線追跡法的數值計算過程, 僅需基于流場數據建立密度代理模型, 能大幅簡化數據集的構建工作。 另外, 由于激波后流場氣體密度的變化較為規律, 在訓練時僅使用了167組數據即可實現比IPSO-BP網絡模型預測法更高的預測精度, 能大幅度簡化建立預測模型的過程, 具有較高的適用性。

4 結 論

當星光穿過高超聲速飛行器外的復雜非均勻流場時會發生氣動光學效應, 導致其傳播方向發生偏折, 引起星光矢量測量誤差進而影響天文導航的精度。 為消除氣動光學效應所導致的星光偏折對天文導航的影響, 通過分析高超聲速飛行器外流場的密度分布特性以及星光在流場中的傳播機理, 提出了一種基于解析模型的高精度氣動光學星光偏折預測方法。 與傳統方法相比, 所提方法由于綜合了偏折角解析模型, 緩解了傳統數據擬合模型預測法過分依賴訓練數據樣本的問題, 能取得更好的預測效果; 另外, 所提方法避免了利用光線追跡法計算偏折角并構建數據集的過程, 解決了傳統方法構建數據集時計算量大、 步驟繁瑣并且會引入額外計算誤差的問題, 能大幅簡化預測模型的構建過程。 仿真驗證表明, 所提方法在建立預測模型時能大幅減少所需的訓練樣本數量, 并且具有比傳統數據擬合模型法更高的預測精度, 可實現亞角秒級別的預測精度, 具有較高的適用性。

另外, 高超聲速環境下, 星敏感器在觀星過程中還面臨著高溫氣體輻射對星光的遮掩問題以及大氣層折射效應所導致的星光矢量測量誤差, 也會影響天文導航的精度。 因此, 后續將進一步研究高溫氣體輻射以及大氣折射對星光矢量測量精度的影響情況, 以進一步提高天文導航系統在高超聲速飛行器中的應用精度。