楊 博，李嘉興，楊照華，劉玉雙

(1. 北京航空航天大學(xué)宇航學(xué)院，北京 100191；2. 北京航空航天大學(xué)儀器科學(xué)與光電工程學(xué)院，北京 100191； 3. 北京機電工程總體設(shè)計部，北京 100854)

0 引 言

高超聲速飛行器近些年來不斷朝著射程更遠、精度更高的方向發(fā)展[1]，其導(dǎo)航技術(shù)是重要的研究課題之一。目前國外的高超聲速飛行器以GPS和慣性組合導(dǎo)航為主，而GPS具備的“三定一弱”(特定信標(biāo)以特定頻率發(fā)播特定格式的導(dǎo)航電文，導(dǎo)航信號微弱)特性再加上軍事用途的約束，因此，探索適合于軍用的自主導(dǎo)航技術(shù)成為了目前我國對該領(lǐng)域攻關(guān)研究的主要方向。天文導(dǎo)航是一種完全的自主導(dǎo)航方法，將其與慣性組合導(dǎo)航在特性上取長補短大大提高了慣性系統(tǒng)的精度和適應(yīng)度。而將其應(yīng)用于高超聲速飛行器上仍處于研究階段，氣動光學(xué)效應(yīng)是首當(dāng)其沖的問題。飛行器周圍的大氣被加熱，流場發(fā)生巨大改變，產(chǎn)生激波、附面層等復(fù)雜流場，光線通過湍流后致使星圖產(chǎn)生偏移、抖動和模糊等負面效果[1]，嚴重影響導(dǎo)航能力。對此類問題的研究通常分為兩個步驟：首先研究高超聲速飛行器周圍高速流場的氣動光學(xué)效應(yīng)，然后對湍流和運動模糊星圖進行補償算法研究。

當(dāng)前研究高速流場中星光氣動光學(xué)效應(yīng)的方法主要是通過數(shù)值模擬流場或?qū)嶒灉y量真實流場的方式進行。數(shù)值模擬的方法通過幾何光學(xué)、物理光學(xué)、統(tǒng)計光學(xué)的方法模擬光線穿過高速流場后的波前；而實驗測量的方法則是直接測量穿過真實流場后的光線波前數(shù)據(jù)，再由波前計算出點擴散函數(shù)，進而分析因流場帶來的光線偏折、像模糊等影響。例如Pond和Sutton[2]采用標(biāo)準(zhǔn)k-ε湍流模型對三維凸臺周圍的流場建立了數(shù)值仿真，利用相位差、斯特列爾比等參數(shù)對氣動光學(xué)效應(yīng)進行了評價；Gao等[3]利用流場可視化方法重構(gòu)氣體密度場，對超聲速湍流邊界層的氣動光學(xué)效應(yīng)進行了時間演化方面的研究。

近幾年圖像復(fù)原技術(shù)作為一種解決高超聲速飛行器氣動光學(xué)效應(yīng)問題的方法得到相當(dāng)?shù)闹匾?。洪漢玉等[4]提出了一系列針對紅外圖像的氣動光學(xué)復(fù)原算法，主要研究了基于極大似然估計準(zhǔn)則的正則化復(fù)原算法，文獻[4]提出了針對圖像與點擴散函數(shù)特點的非線性各向異性正則化函數(shù)，文獻[5]提出了基于泊松模型的多幀圖像復(fù)原算法。正則化復(fù)原方法近年來在國際上發(fā)展迅速，F(xiàn)ergus等[6]提出基于變分貝葉斯的圖像盲復(fù)原方法，Krishnan和Fergus[7]用超拉普拉斯分布來擬合圖像梯度的重尾分布，Ohkoshi等[8]提出結(jié)合全變差模型和Shock濾波器的盲復(fù)原算法，Pan等[9]提出基于暗通道稀疏先驗的盲復(fù)原算法。這些復(fù)原算法的研究思路主要是針對灰度、梯度、模糊核等特征的分布特性設(shè)計正則函數(shù)，通過最優(yōu)化求解得到清晰圖像。

上述的正則化復(fù)原算法研究對象均以自然圖像為主，缺乏針對星圖的研究。同時圖像盲復(fù)原算法的關(guān)鍵是要根據(jù)圖像豐富的邊緣信息估計出準(zhǔn)確的模糊核。對于高超流場影響下的星圖基本失去了邊緣信息，因而利用一般方法是無法建立準(zhǔn)確的模糊核用于高精度復(fù)原。為此本文針對高超星圖特點提出一種利用先驗知識提取準(zhǔn)確模糊核建立的半盲復(fù)原算法，并進行了計算機半物理仿真校驗，獲得了理想的結(jié)果，為高超聲速飛行器利用星光進行自主導(dǎo)航提供了理論依據(jù)。

1 高超湍流場的運動星圖模糊核

基于高超湍流的大渦模擬流場數(shù)據(jù)，利用幾何光學(xué)和物理光學(xué)，統(tǒng)計計算湍流場影響下的星圖模糊核；并根據(jù)湍流疊加模型，引入隨機運動模糊，建立具有氣動光學(xué)效應(yīng)和運動模糊的高超湍流的運動模糊核。

1.1 高超湍流的大渦模擬

星光經(jīng)過激波、湍流邊界層、彈體冷卻層等之后，光學(xué)窗口外復(fù)雜高速流場產(chǎn)生的氣動光學(xué)效應(yīng)使星敏感器接收的圖像存在嚴重的偏移、抖動和模糊。由激波引起的穩(wěn)定偏折可由文獻[1]中的經(jīng)驗公式計算。而本文只重點研究湍流模糊與運動模糊問題。

為了盡可能準(zhǔn)確地揭示高超湍流場中星光光線傳播的特性，使用大渦模型對湍流流場進行模擬，其流場密度分布如圖1所示。

1.2 湍流場中的星光傳輸效應(yīng)

基于上述大渦流場模擬結(jié)果，利用物理光學(xué)方法計算流場的光學(xué)傳輸效應(yīng)。根據(jù)Gladstone-Dale定律將密度分布轉(zhuǎn)換為折射率分布，利用光線追擊法計算平面波的光程差(Optical path difference,OPD)，如圖2所示。

圖2 光程差示意圖Fig.2 Optical path difference

w(x,y)=kΔOPD(x,y)

(1)

(2)

式中：x和y為像面坐標(biāo)系中的橫縱坐標(biāo)，k=2π/λ，λ為星光波長，a為波面的振幅分布，j為虛數(shù)單位。

1.3 湍流影響的退化星圖建模

(3)

(4)

(5)

(6)

式中:F(·)表示快速傅里葉變換，上標(biāo)*表示復(fù)共軛，° 表示矩陣點乘(即對應(yīng)位置元素相乘)。

設(shè)高超聲速飛行器飛行過程中光軸抖動產(chǎn)生運動模糊效應(yīng)，并有隨機運動模糊的模糊核為hM，將其與高超湍流模糊核卷積即可得到高超湍流場運動模糊核，如圖3所示。

h=hH?hM

(7)

其中，?為空間卷積運算。根據(jù)湍流凍結(jié)理論，認為短曝光的湍流退化圖像的模糊核線性位移不變，可將湍流退化過程簡單表示為目標(biāo)圖像與模糊核的卷積過程，同時考慮噪聲有：

g(x,y)=f(x,y)?h(x,y)+n(x,y)

(8)

式中：g(x,y)為退化圖像，f(x,y)為清晰的目標(biāo)圖像，n(x,y)為加性噪聲。為方便表示，下文將寫作g,f,h,n。

圖3 高超湍流運動模糊核Fig.3 Kernel of hypersonic turbulence and motion

2 高超湍流星圖的復(fù)原預(yù)處理

對星光通過大渦模擬的湍流場時，受到衰減后的星圖點擴散函數(shù)圖3，提出運用三種復(fù)原預(yù)處理技術(shù)：低通去噪、星點粗提取、模糊核估計，以此提高模糊核的估計準(zhǔn)確度，并保證獲得最佳的復(fù)原效果。

2.1 低通濾波

圖4 大渦模擬流場星圖Fig.4 Astronomical image simulation from LES

由式(8)隨機生成模擬流場星圖，如圖4所示。從圖4可以看出，星點有大量的分布信息淹沒于隨機噪聲中，直接采用閾值分割，會丟失過多的星點能量，造成星點亮度降低質(zhì)心提取不準(zhǔn)確，甚至可能將抹去低亮度星點的全部信息，使得星點跟蹤失敗。而從圖3可以看出，湍流模糊核分布比較平滑，信息變化緩慢，說明模糊圖像具有較強的低頻特性。因此采用參數(shù)適當(dāng)?shù)母咚沟屯V波器先對原始星圖進行濾波處理，得到的去噪圖像如圖5所示。

圖5 低通濾波后的高超星圖Fig.5 Image after LPF

比較圖5與圖4可以看出，經(jīng)過合適的低通濾波器預(yù)先處理，圖像去除了大量高頻噪聲，而大部分星點信息被分離保留下來。但由圖5可知，星點能量處于分散狀態(tài)，并且鄰星之間存在能量疊加區(qū)域(左上部的兩顆星)，為了獲得較高的質(zhì)心提取精度，需對圖像進行復(fù)原，集中星點能量。

分割域值選取，是通過精確地分析統(tǒng)計噪聲的分布P，得到噪聲的均值與方差，同時根據(jù)大量濾波后星圖灰度的統(tǒng)計特性，設(shè)計灰度分布擬合函數(shù)：

(9)

圖6 灰度分布擬合示意圖Fig.6 Fitting of gray distribution

2.2 星點初提取

圖7 星點提取方法對比Fig.7 Comparison of star point extraction methods

從圖7中上下兩個白色方框內(nèi)的分割結(jié)果可以分別看出，改進方法在平滑多峰和避免重疊方面的效果，利用此提取方法可以大大避免由氣動光學(xué)效應(yīng)造成的星點提取錯誤。盡管如此，初步提取的星點坐標(biāo)并非是由該星點的全部圖像信息提供，因而坐標(biāo)準(zhǔn)確度較低，并且能量較弱或相鄰很近的星可能會無法提取。

2.3 模糊核提取

天文圖像的模糊核提取與光學(xué)系統(tǒng)點擴散函數(shù)測量問題中的點脈沖法十分相近，但在氣動光學(xué)作用下，需解決模糊核重疊、圖像邊界處不完整及成像噪聲等問題。

假設(shè)星圖中有M顆星，坐標(biāo)是(xm,ym)，(m=1,2,…,M)，能量為Im的星點圖像為Imδ(x-xm,y-ym)，式中δ(x,y)為脈沖函數(shù)，則原始星圖可表示為：

(10)

經(jīng)過點擴散函數(shù)h的退化，并加入噪聲后的退化圖像可表示為：

g(x,y)=f?h+G(μn,σn)=

(11)

假設(shè)1. 閾值分割階段可以估計出準(zhǔn)確的噪聲均值。

當(dāng)假設(shè)1成立時，在退化圖像中去除噪聲均值，即g-μn后，邊界延拓L/2，以(xm,ym)為中心，提取L×L的圖像為：

Ωm(x,y)=Imh(x,y)+G(0,σn)+

(12)

式中：L為估計模糊核的支持域邊長，Δxmm′=xm′-xm，Δymm′=ym′-ym，D(sm,sm′)為m星到m′星之間的距離。

由式(12)可知，Ωm(x,y)中包含m星對應(yīng)的點擴散函數(shù)h、零均值噪聲以及重疊的其他星的干擾能量，去除重疊區(qū)域干擾的處理如下：

將Ωm(x,y)寫為極坐標(biāo)形式Ωm(φ,r)，在剔除不可用的區(qū)域之前，設(shè)保留區(qū)間為φn=U=[0,2π)。

1)剔除Ωm中與邊界重疊的區(qū)域，如圖8(a)所示，其中虛線為擴展前的邊界：如果(xm,ym)與某邊界b的間距Db(sm)

θmb=arccos(2Db(sm)/L-1)

(13)

更新保留區(qū)間為：

φm∶=φm(φmb-θmb,φmb+θmb)

(14)

式中：“”表示集合作差。剔除后的剩余區(qū)域都不與邊界重疊。

圖8 半圓心角示意圖Fig.8 Semi central angle

(15)

3)更新保留區(qū)間為：

φm∶=φm(φmm′-θmm′,φmm′+θmm′)

(16)

將與邊界重疊和模糊重疊的圓心角區(qū)域剔除后，Ωm可更新為：

(17)

假設(shè)2. 利用式(17)成功剔除了m星與其它星的模糊重疊區(qū)域。

當(dāng)假設(shè)2成立時，m星的能量近似為：

(18)

并以其作為歸一化系數(shù)，歸一化Ωm：

(19)

(20)

加權(quán)融合后即得到模糊核：

(21)

(22)

2.4 模糊核特性分析

式(12)在滿足假設(shè)2的條件下，可改寫為：

(23)

期望和方差為：

(24)

(25)

將式(20)、(23)代入式(21)，可得：

(26)

則期望和方差為：

E(hp(r,φ))=h(r,φ)

(27)

(28)

在滿足假設(shè)1和假設(shè)2的條件下，式(27)表明模糊核先驗估計hp是退化真實模糊核h的無偏估計。式(28)表明融合后的方差一定小于每顆星先驗區(qū)域內(nèi)的方差：

(29)

式(29)表明融合過程能降低其中的噪聲方差，起到降噪作用。

由式(28)可得：

(30)

3 高超星圖正則化非盲復(fù)原

基于提取并融合的模糊核，并將其用于非盲復(fù)原這一過程稱為半盲復(fù)原。一般正則化技術(shù)具有一定的快速性及模糊核估計誤差魯棒性。

分析目標(biāo)的先驗知識構(gòu)造正則化項，用以約束最優(yōu)化求解過程。將清晰的星圖的灰度分布統(tǒng)計經(jīng)對稱化處理后如圖9所示。

圖9 星圖灰度分布示意圖Fig.9 Gray distribution of astronomical image

圖10 星圖梯度分布示意圖Fig.10 Gradient distribution of astronomical image

(31)

(32)

交替計算f與輔助變量，其中f可由下式計算：

(33)

式中：F-1(·)為快速傅里葉逆變換，δx和δy分別為水平和垂直一階差分算子。由式(33)更新f后，對其進行半正定約束：

(34)

同時，輔助變量的子問題表示為：

(35)

式中：w∈{a,b1,b2},v∈{f,δxf,δyf}。采用文獻[10]給出的LUT算法快速求解。圖11為總體復(fù)原算法流程圖，記為算法1。

圖11 復(fù)原算法流程圖Fig.11 Restoration algorithm

4 仿真校驗

4.1 仿真條件

4.2 仿真結(jié)果

在模糊核提取過程中，剔除每顆星與邊界、其他星重疊的區(qū)域后，剩余部分可由圖12所示，白色框中的部分即為應(yīng)用于融合模糊核的圖像。

圖12 模糊核融合區(qū)域示意圖Fig.12 Fusion area of kernel

利用算法1和對比算法分別對湍流模糊星圖進行復(fù)原，估計模糊核和復(fù)原圖像分別如圖13和圖14所示。對比可以看出，CLSF算法對模糊核估計準(zhǔn)確度的要求較高，很難復(fù)原出質(zhì)量好的圖像。Krishnan等[11]的算法使用L1先驗，平滑度較高，可以抑制絕大多數(shù)噪聲，但同時復(fù)原的星點較粗，能量不夠集中，可能會使相鄰星點融為一體，難以區(qū)分相鄰的兩顆星。Pan等[12]的算法使用L0先驗，雖然可以得到更細致的星點，但過度保留了圖像梯度，使復(fù)原圖像中依然存在大量高頻噪聲，易使星點誤識別率高。算法1中使用的Lp先驗?zāi)Ｐ?，效果介于L1與L0之間，其既可抑制噪聲，同時也集中能量，提高了星點質(zhì)心提取精度。

圖13 模糊核對比Fig.13 Comparison of kernel

圖14 復(fù)原星圖對比Fig.14 Comparison of deblurred image

表1 不同復(fù)原算法的平均PSNR對比Table 1 Comparison of mean PSNR between different algorithm

如表2所示是對復(fù)原后的星圖進行星點質(zhì)心提取的精度比較。由復(fù)原圖像計算識別出的星點質(zhì)心坐標(biāo)，并與原始星點坐標(biāo)數(shù)據(jù)進行對比，將坐標(biāo)誤差小于2個像素記為識別正確，其個數(shù)占總星數(shù)的比例為識別正確率。同時將正確識別的星點坐標(biāo)與真實坐標(biāo)進行對比，統(tǒng)計二者的平均誤差角距。從表2可以看出，算法1的星點識別最精確，從而印證了此方法對模糊核的估計更準(zhǔn)確，并且復(fù)原精度更高，使得復(fù)原圖像更接近實際圖像，在提取星點時所受到的噪聲、模糊核多峰結(jié)構(gòu)、模糊重疊等影響因素更小，有助于提高識別精度。

表2 星點識別精度對比Table 2 Comparison of star recognition accuracy

采用文獻[13]給出的誤差率作為評價模糊核的估計準(zhǔn)確度，對不同算法統(tǒng)計積累誤差率繪制曲線如圖15所示。圖中曲線越高說明模糊核估計得越準(zhǔn)確，因此可以看出該算法1的模糊核估計準(zhǔn)確度最高。例如在誤差率1.2的指標(biāo)下，對比算法分別有50%和25%的復(fù)原圖像達到要求，而算法1有超過80%達到要求，由此可以明顯看出該算法的優(yōu)越性。

圖16為所有復(fù)原圖像的累積誤差率與星點質(zhì)心識別正確率的關(guān)系示意圖。從圖16可以看出，隨著誤差率的積累，對較差復(fù)原效果的容忍度越高，星點識別正確率自然隨之降低。在實際使用時可以根據(jù)不同識別正確率的任務(wù)需求，結(jié)合圖16和圖15得到不同算法的占比成功率，具體示例如表3所示。

圖15 積累誤差率示意圖Fig.15 Accumulated error rate

圖16 積累誤差率與識別精度的關(guān)系示意圖Fig.16 Relationship between accumulated error rate and recognition accuracy

識別正確率/%占比成功率/%文獻[11]文獻[12]算法19026253803446770502484

由表3的對比分析可知，在不同識別正確率需求下，算法1滿足任務(wù)的復(fù)原圖像占比均明顯高于對比算法，此對比進一步驗證了該算法對能有效提高星點識別能力。綜合上述對比可以看出，與傳統(tǒng)盲復(fù)原算法中僅利用圖像邊緣信息估計模糊核并復(fù)原的思想相比，針對高超湍流模糊復(fù)原問題，算法1能夠有效結(jié)合星圖特點及高超模糊核特點從而提取出更精確的模糊核，并且利用星圖灰度分布的特點使得復(fù)原效果達到最佳，因此星點質(zhì)心提取的正確率和精度更高。

5 結(jié) 論

針對高超聲速飛行器在平流層使用天文導(dǎo)航受氣動光學(xué)及運動模糊干擾的關(guān)鍵問題，本文結(jié)合氣動光學(xué)效應(yīng)對星圖復(fù)原算法進行深入研究，填補了氣動光學(xué)研究至復(fù)原算法之間的空白，得到以下結(jié)論：1)基于區(qū)域剔除與融合操作，模糊核提取方法對高超湍流模糊星圖模糊核的估計精度高于其他盲復(fù)原算法。2)基于星圖灰度與梯度的稀疏特性，算法1對高超湍流模糊星圖的復(fù)原精度高于其他復(fù)原算法。3)經(jīng)算法1復(fù)原后的星點識別正確率可達78%，誤差角距低至0.0249°。因此若將本文提出的復(fù)原算法應(yīng)用于高超聲速飛行器的天文導(dǎo)航領(lǐng)域，可大幅提升導(dǎo)航適應(yīng)能力及導(dǎo)航精度。同時可為高超聲速飛行器星圖復(fù)原的深入研究提供思路。