李嘉興,王大軼,鄂薇,葛東明

北京空間飛行器總體設(shè)計部,北京 100094

高超聲速飛行器近些年來的發(fā)展方向是射程更遠、精度更高[1],這無疑需要更高精度的導航能力?？紤]到這種軍事意義十分明顯的飛行器要做到自主導航的迫切需求,探索適合于高超聲速飛行器的自主導航技術(shù)成為了目前中國對該領(lǐng)域攻關(guān)研究的主要方向之一。天文導航是一種可靠的自主導航方法,在航天領(lǐng)域有很成熟的應(yīng)用經(jīng)驗,將其與慣性敏感器進行組合導航可以達到長期穩(wěn)定且高精度的導航效果。要將天文導航應(yīng)用于高超聲速飛行器仍處于研究階段,主要受到了氣動光學效應(yīng)的影響。因為在飛行器高速飛行過程中,星敏感器窗口的大氣被加熱,流場發(fā)生巨大改變,產(chǎn)生激波、附面層等復(fù)雜流場,光線通過湍流后發(fā)生偏折,致使星圖產(chǎn)生偏移、抖動和模糊等負面效果[1],嚴重降低導航能力。為了抑制氣動光學效應(yīng),需要首先模擬在高超湍流成像條件下的天文導航圖像,在此基礎(chǔ)上研究圖像復(fù)原算法。

目前針對氣動光學效應(yīng)的研究方法主要包括數(shù)值模擬和實驗測量2種方法。數(shù)值模擬方法是利用幾何光學、物理光學、統(tǒng)計光學等物理模型,在計算機中模擬光線穿過高速流場后的畸變波前,進而模擬出成像時的點擴散函數(shù)?，F(xiàn)有的湍流數(shù)值模擬方法包括直接模擬法、大渦模擬法和雷諾平均法。直接模擬法能獲得瞬態(tài)湍流但計算量巨大；雷諾平均法計算量小但無法模擬出瞬態(tài)湍流對光線的影響效果；大渦模擬法計算量居中且能夠表達瞬時流場中的大渦結(jié)構(gòu),因此本文選用該方法進行氣動光學效應(yīng)仿真。近些年有不少數(shù)值仿真研究,例如Pond和Sutton[2]采用標準k-ε湍流模型對三維凸臺周圍的流場建立了數(shù)值仿真,利用相位差、斯特列爾比等參數(shù)對氣動光學效應(yīng)進行了評價。利用大渦模擬法研究氣動光學效應(yīng)已為主流,White和Visbal[3]對熱壁面和冷壁面可壓縮湍流邊界層引起的航空光學像差進行了大渦模擬,還表明當壁面受熱時光程差會增加。Kamel等[4]驗證了一種帶有壁面模型的大渦模擬方法的有效性。實驗測量方法是直接測量光線穿過真實流場后的光程差數(shù)據(jù),研究氣流參數(shù)對圖像模糊的影響機理。Gordeyev等[5]通過統(tǒng)計實驗數(shù)據(jù)獲得了非隔熱壁面上邊界層處波動密度場的預(yù)測模型。但這些針對氣動光學效應(yīng)的研究主要集中于數(shù)值仿真模型的建立、流場參數(shù)對光學成像指標的影響等方面,很少有面對氣動光學影響下的圖像復(fù)原研究。

近些年利用圖像復(fù)原技術(shù)解決高超聲速飛行器的氣動光學效應(yīng)問題受到了很高的重視。洪漢玉和張?zhí)煨騕6]在極大似然估計準則的基礎(chǔ)上采用正則化方法復(fù)原紅外圖像,文中基于圖像和湍流點擴展函數(shù)的梯度變化先驗知識,提出了針對圖像與點擴散函數(shù)特點的非線性各向異性正則化函數(shù),使其估計模糊核與清晰圖像時能自適應(yīng)地平滑梯度。正則化復(fù)原方法近些年在國際上也備受關(guān)注,學者們利用不同的先驗?zāi)Ｐ驮O(shè)計正則函數(shù)。大多數(shù)文獻針對圖像和梯度的分布特點設(shè)計模型,Fergus等[7]提出基于變分貝葉斯的圖像盲復(fù)原方法,Krishnan和Fergus[8]發(fā)現(xiàn)了圖像梯度的重尾分布并用超拉普拉斯分布來擬合。Li等[9]引入了分裂布雷格曼迭代法用于解決非凸優(yōu)化模型。Ohkoshi等[10]提出結(jié)合全變差模型和Shock濾波器的盲復(fù)原算法。Zhao等[11]利用Lp范數(shù)擬合車牌圖像的分布特征,模糊圖像中同時考慮了離焦與運動模糊問題,并用交替方向乘子(ADMM)算法求解復(fù)原模型的最優(yōu)值。Sun等[12]利用迭代支持檢測方法精細化模糊核,以降低圖像噪聲對模糊核估計結(jié)果的精度影響。Pan等[13]提出基于暗通道稀疏先驗知識的盲復(fù)原算法。而根據(jù)低秩先驗原理,Ma等[14]設(shè)計了基于加權(quán)核范數(shù)正則化和全變差正則化的復(fù)原模型,前者能夠降低由相似局部圖像堆疊而成的矩陣的秩,進而約束去噪與模糊復(fù)原過程。以上這些圖像模糊復(fù)原方法大多是為自然圖像設(shè)計,而為了解決氣動光學大動態(tài)對天文導航圖像影響,需要針對這種圖像的特點進行專門的復(fù)原模型設(shè)計,與氣動光學的研究成果相結(jié)合。

模糊核的先驗信息有助于獲得更接近真實值的模糊核估計結(jié)果[15]。文獻[16]為解決氣動光學效應(yīng)專門設(shè)計了星圖復(fù)原算法,利用從圖像中提取的先驗?zāi)：诉M行非盲復(fù)原。但該方法無法適應(yīng)先驗?zāi)：颂崛【鹊突驘o先驗信息等情況。為了解決這一問題,本文提出一種新的改進算法,在對氣動光學影響下星圖特點研究的基礎(chǔ)上,把先驗?zāi)：伺c盲復(fù)原方法相結(jié)合,設(shè)計了新的正則化圖像復(fù)原模型,既具備了盲復(fù)原算法可以在未知模糊核情形下的復(fù)原能力,又具備了非盲復(fù)原算法對先驗?zāi)：酥袌D像信息的利用能力,提高算法應(yīng)用時的魯棒性和復(fù)原精度。

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

基于高超湍流的大渦模擬(LES)流場數(shù)據(jù),利用幾何光學和物理光學,湍流場影響下的星圖模糊核,引入隨機運動模糊,建立具有氣動光學效應(yīng)和運動模糊的高超湍流的運動模糊核。

1.1 高超湍流的大渦模擬

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

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

1.2 湍流影響的星圖模擬

基于上述大渦流場模擬結(jié)果,利用物理光學原理,采用光線追擊法,計算平行光穿過流場后的光程差,如圖2所示。根據(jù)文獻[1]中的方法,由光程差進一步計算氣動光學效應(yīng)造成的高超湍流模糊核hH。

圖2 光程差

當高超聲速飛行器飛行時受到氣流擾動和機體振動導致光軸發(fā)生抖動,造成的運動模糊核記為hM,將其與高超湍流模糊核卷積即可得到高超湍流場運動模糊核h,如式(1)所示,卷積結(jié)果如圖3所示。

圖3 高超湍流場運動模糊核

h=hH?hM

(1)

式中：?為空間卷積運算,根據(jù)湍流凍結(jié)理論,在短曝光前提下可將模糊退化過程表示為目標圖像與模糊核的卷積過程,同時考慮噪聲:

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

(2)

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

2 先驗?zāi)：颂崛∨c分析

圖像復(fù)原是模擬模糊圖像過程的逆過程,在實際導航過程中,需要從模糊圖像當中提取有用信息,利用復(fù)原算法獲得清晰圖像。區(qū)別于自然圖像當中利用盲復(fù)原算法估計模糊核,基于天文導航圖像的獨有特點,可以先從星圖當中提取模糊核的先驗信息來指導復(fù)原過程,稱之為先驗?zāi)：恕?/p>

2.1 先驗?zāi)：颂崛?/h3>

天文圖像的模糊核提取與光學系統(tǒng)點擴散函數(shù)測量問題中的點脈沖法十分相近,星圖中每顆星的圖像與模糊核非常相近,借助這一特性從模糊星圖中提取模糊核的先驗信息,為復(fù)原提供依據(jù)。但在氣動光學作用下,需解決模糊核重疊、圖像邊界處不完整及成像噪聲等問題。

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

(3)

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

xm,y-ym)+G(μn,σn)

(4)

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

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

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

(5)

式中：L為估計模糊核的支持域邊長,單位為像素個數(shù)；Δxmm′=xm′-xm；Δymm′=ym′-ym；D(sm,sm′)為m星到m′星之間的距離。

由式(5)知,Ωm(x,y)中包含m星對應(yīng)的點擴散函數(shù)h、零均值噪聲以及重疊的其他星的干擾能量,將Ωm(x,y)寫為極坐標形式Ωm(r,φ),按照文獻[16]中的算法提取出先驗?zāi)：恕?/p>

假設(shè) 2利用文獻[16]的算法成功剔除了m星與其他星的模糊重疊區(qū)域。

加權(quán)融合后的模糊核表示為

(6)

(7)

(8)

利用文獻[16]中的先驗?zāi)：颂崛∷惴?能夠剔除相鄰星之間圖像的重疊區(qū)域,從而避免了大擾動條件下模糊核重疊的問題。對圖像進行邊界延拓可以使圖像邊界處模糊核不完整時也能正常提取有用信息。

2.2 先驗?zāi)：颂匦苑治?/h3>

式(5)在滿足假設(shè)2的條件下,可改寫為

(9)

期望和方差為

(10)

(11)

將式(7)、式(9)代入式(6)可得

(12)

則期望和方差為

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

(13)

(14)

在滿足假設(shè)1和假設(shè)2的條件下,式(13)表明模糊核先驗估計hp是退化真實模糊核h的無偏估計。且由式(14)可以看出融合后的hp方差一定不會大于融合之前,而且所利用的圖像信息越多、能量越豐富,方差越小,去噪效果越強。下面定性描述hp與h的差異,為用于表征先驗?zāi)：说目尚湃味取?/p>

由式(14)得到hp每點處的方差為

(15)

為描述hp與h的誤差,采用二者差值的均方誤差

(16)

將式(15)代入式(16)可得

(17)

同時,均方誤差亦可表示為

(18)

(19)

將式(19)代入式(18)可得

(20)

(21)

hp-h為融合后的噪聲,與h相獨立,則

(22)

式中:

(23)

(24)

(25)

將式(17)、式(25)代入式(20)可得hp與h相似度的估計為

(26)

在估計相似度時,只需要先驗?zāi)：薶p和噪聲方差估計值這些已知量,不需要實際模糊核h。因此式(26)的估計方法適用于在未知h的前提下盲復(fù)原模糊星圖時,衡量所提取的先驗?zāi)：说目尚湃纬潭?為復(fù)原算法在利用hp時提供參數(shù)設(shè)計依據(jù)。

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

3.1 半盲復(fù)原模型設(shè)計

由式(2)所示的圖像退化模型可知,圖像盲復(fù)原的過程用數(shù)學可以表示成為求解未知量f與h使得后驗概率P(f,h|g)取得最大值,即可表示為如下最大后驗概率問題:

(27)

式中:f*表示清晰圖像的最優(yōu)估計;h*表示模糊核的最優(yōu)估計。根據(jù)貝葉斯定理,后驗概率P(f,h|g)可以轉(zhuǎn)化為

(28)

由于在圖像復(fù)原過程中模糊圖像g為已知量,因此式(28)中的分母是一個固定的歸一化常數(shù),可以忽略,于是可以得到

P(f|g,h)∝P(g|f,h)P(f)P(h)

(29)

將式(29)代入式(27)可得

(30)

(31)

綜上,本章提出的高超星圖正則化半盲復(fù)原模型設(shè)計為

(32)

3.2 最優(yōu)化求解

為了使模糊核估計更精確,總體求解思路采用文獻[7]中提出的多尺度框架策略,從分辨率由低到高進行逐層復(fù)原,特別針對嚴重模糊的圖像,能夠有效避免陷入局部最優(yōu)。式(32)是一個非凸優(yōu)化問題,需要交替迭代g與f進行更新。下面為該算法的具體過程,記為算法 1。

3.2.1 目標圖像估計

在更新目標圖像f的階段,固定上一次迭代得到的h,則f的子問題可轉(zhuǎn)化為如下最優(yōu)化問題:

(33)

式中：Hf=h?f。

式(33)為Lp范數(shù)的優(yōu)化問題,采用半二次懲罰技術(shù),引入輔助變量a=f,b1=δx?f,b2=δy?f,得到新的目標函數(shù)為

(34)

交替計算f與輔助變量。式中：δx和δy分別表示水平和垂直一階差分算子,δx?f和δy?f簡寫為δxf和δyf。懲罰參數(shù)α越大時式(34)的解越接近式(33)的解,最終f可由式(35)計算:

(35)

(36)

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

(37)

式中：w∈{a,b1,b2}；v∈{f,δxf,δyf}。采用文獻[8]給出的LUT(Look up Table)算法快速求解。

3.2.2 模糊核估計

在更新模糊核h的階段,固定上一次迭代得到的f,則h的子問題可轉(zhuǎn)化為如下最優(yōu)化問題:

(38)

式中：Fh=f?h。

式(38)為L1范數(shù)的優(yōu)化問題,利用分裂Bregman迭代來求解。采用半二次懲罰技術(shù),引入輔助變量c=h,式(38)可改寫為迭代框架

(39)

t:=t+h-c

(40)

固定c、t,更新h的方法為

(41)

利用式(41)更新模糊核后對其進行歸一化

(42)

并為了提高抗噪魯棒性,施加動態(tài)閾值約束

(43)

式中：δh表示閾值參數(shù)。固定h、t,利用快速收縮算法求解c的最優(yōu)解:

(44)

綜上,利用算法 1復(fù)原圖像的過程如圖4所示。先利用經(jīng)預(yù)處理后的模糊圖像進行目標圖像估計,再利用先驗?zāi)：撕拖嗨贫冗M行模糊核估計,這一過程作為一個循環(huán)。每次循環(huán)結(jié)束前將參數(shù)α擴大αInc倍,但不超過αMax,一共進行N次循環(huán)。再將算法 1代入文獻[17]中的金字塔算法,來提高大尺度模糊時算法的魯棒性。

圖4 算法1流程圖

4 算法收斂性驗證

改進的半盲復(fù)原算法中,目標圖像估計階段,hp給定,f的子問題為凸優(yōu)化,其他參數(shù)均由查表法獲得最優(yōu)值,收斂性不必證明。因此只有計算h時的L1約束項通過Bregman分裂法進行近似,因此下文將證明其收斂性。式(38)的遞推式可寫為

(45)

式(45)為L1范數(shù)的優(yōu)化問題,利用分裂Bregman迭代來求解。采用半二次懲罰技術(shù),引入輔助變量c=h,式(45)轉(zhuǎn)化為無約束分裂形式

(46)

式(46)可改寫為迭代框架

(47)

(48)

tn+1=tn+(hn+1-cn+1)

(49)

式(47)～式(49)子問題是凸性可微分的,可得一階最優(yōu)條件為

0=FT(Fhn+1-g)+θη1hn+1-θη1hp+

(50)

0=rn+1+2β(cn+1-hn+1-tn)

(51)

tn+1=tn+(hn+1-cn+1)

(52)

(53)

(54)

并且此解是唯一解。

證明：因為式(46)存在至少一個解h*,其一階最優(yōu)條件滿足:

(55)

0=FT(Fh*-g)+θη1h*-θη1hp+β(1-θ)×

(56)

0=r*+2β(c*-h*-t*)

(57)

t*=t*+(h*-c*)

(58)

(59)

(60)

(61)

(62)

將式(61)與式(62)相加,得到

(63)

用式(52)減去式(58),可得到

(64)

對式(64)的等號兩邊求平方,得到

(65)

(66)

將式(66)代入式(63)中,得到

(67)

(68)

(69)

(70)

對式(70)的等號兩邊分別從0到N進行求和運算,可以得到

(71)

式(71)中的各個項都是非負性的,于是有

(72)

假設(shè)β,η1,η2,η3>0, 0<θ<1,并且使得N趨近于無窮,對于不等式(72)得到結(jié)論

因此可以得到

(73)

(74)

(75)

(76)

(77)

(78)

(79)

(80)

hn-h*=0

(81)

(82)

hn-h*=0

(83)

(84)

(85)

綜合式(80)～式(85),可得

(86)

由式(55),可得

(87)

再結(jié)合式(82),即可證明收斂至h*,下面證明解得唯一性：定義

(88)

E(hni)>ζE(h*)+(1-ζ)E(hni)≥E(ζh*+

(1-ζ)hni)=E(v)≥min{E(h):

(89)

(90)

存在矛盾,因此h*是唯一最優(yōu)解。

5 仿真驗證

5.1 仿真條件

5.2 仿真結(jié)果

表1 不同算法的平均PSNR對比

從圖5、圖6中可以直觀地得到以下對比結(jié)果。CLSF是一種非盲復(fù)原算法,因此直接使用先驗?zāi)：俗鰪?fù)原,但這種方法對模糊核估計準確度的要求較高,當模糊核不夠準確時復(fù)原圖像精度較低,沒有先驗?zāi)：藭r更無法進行復(fù)原。Krishna等[18]的算法使用L1先驗,對模糊核和圖像的平滑度略高,可以抑制絕大多數(shù)噪聲,但模糊核的形狀比較窄,同時復(fù)原的星點較粗,能量不夠集中,出現(xiàn)了相鄰星互相重疊的現(xiàn)象,有時難以區(qū)分。Pan和Su[19]的算法使用L0先驗,對圖像表征的稀疏程度更高,雖然復(fù)原的圖像能量更集中,但過度保留了圖像和模糊核中梯度的分布,使復(fù)原圖像和模糊核中依然存在大量高頻噪聲,可能會對星圖識別造成干擾。文獻[16]是采用了Lp先驗?zāi)Ｐ偷姆敲?fù)原方法,平滑程度介于L1和L0之間,在抑制噪聲的同時也使得目標星點能量更集中。相比文獻[16]直接使用先驗?zāi)：诉M行復(fù)原,算法 2加入了模糊核估計項,是一個盲復(fù)原方法,由于其僅依靠傳統(tǒng)的模糊核估計項,得到的模糊核精度依然不高,但同樣采用了Lp先驗?zāi)Ｐ褪沟闷鋸?fù)原效果低于文獻[16]又強于文獻[18-19]。在先驗?zāi)：说膸椭?算法 1估計出的模糊核與實際模糊核相差不多,并且由于其具備盲復(fù)原的功能,不完全依賴先驗?zāi)：?能夠在復(fù)原過程中對模糊核進行調(diào)整,因此其復(fù)原效果更佳。

圖5 不同算法的模糊核對比

圖6 不同算法的復(fù)原星圖對比

星圖復(fù)原后提取質(zhì)心坐標,與真實坐標進行對比,若二者相差小于2個像素就認為提取正確,并統(tǒng)計正確識別的個數(shù)占總個數(shù)的識別正確率,如表2所示。造成提取失敗的因素可能有復(fù)原后噪聲被放大,被錯誤的提取出來；模糊程度較大的星沒有正確復(fù)原,沒有被提取出來；復(fù)原后的質(zhì)心比原始質(zhì)心偏移較大；兩顆相鄰的星模糊后重疊在一起,復(fù)原時沒有將二者分離開等等。將正確識別的星點坐標與真實坐標進行對比,統(tǒng)計二者的平均誤差角距如表2所示。從表中的數(shù)據(jù)分析算法的復(fù)原能力,可以得到與前文相似的結(jié)論。同時發(fā)現(xiàn)文獻[18-19]和算法 2等盲復(fù)原方法質(zhì)心提取偏差較大的現(xiàn)象,是因為在缺乏先驗?zāi)：溯o助時估計的模糊核不夠準確,沒有擬合出模糊核中的多峰結(jié)構(gòu),致使復(fù)原出來的星點圖像依然存在多峰結(jié)構(gòu),除主峰以外的次峰會嚴重影響星點質(zhì)心的計算精度。

表2 不同算法的質(zhì)心提取精度對比

采用文獻[20]給出的誤差率作為評價模糊核的估計準確度,對不同算法統(tǒng)計積累誤差率繪制曲線如圖7所示,圖中曲線越高說明模糊核估計得越準確。從圖中可以看出算法 1的模糊核估計準確度更高。

圖7 不同方法的積累誤差率

在6馬赫的飛行速度下,改變飛行高度,利用算法 1復(fù)原后的質(zhì)心偏移曲線如圖8所示,可以看出隨著飛行高度的提升,大氣密度隨指數(shù)降低,高超湍流帶來的模糊強度也逐漸降低,因此對質(zhì)心偏移的影響也逐漸變小?？梢钥闯鋈粝虢档蜌鈩庸鈱W效應(yīng)對天文導航精度的影響,除了圖像復(fù)原之外,還可以提升飛行高度。

圖8 質(zhì)心偏移隨飛行高度變化

綜合上述對比分析,本文提出的半盲復(fù)原算法相較于其他對比算法,對模糊核的估計更加準確,針對高超湍流模糊星圖的復(fù)原效果更好,用于星點質(zhì)心提取的正確率和精度更高。本文的兩種算法相互對比,帶有模糊核先驗項的算法 1比算法 2效果更好,說明在先驗?zāi)：说淖饔弥?估計出的模糊核準確度比完全盲復(fù)原時要更高,從而使得復(fù)原精度也更高。

6 結(jié) 論

針對高超聲速飛行器在平流層中使用天文導航進行自主導航過程中,星圖受氣動光學效應(yīng)及運動模糊干擾的關(guān)鍵問題,深入研究了高超星圖的復(fù)原方法,抑制氣動光學效應(yīng)給天文導航帶來的負面影響。

1) 本文在盲復(fù)原算法的基礎(chǔ)上,利用先驗?zāi)：颂峁┑男畔⒁龑：斯烙嫷倪^程,形成了半盲復(fù)原模型。

2) 設(shè)計了針對復(fù)原模型的優(yōu)化求解算法,交替估計模糊核與清晰圖像,并驗證了該方法的可收斂性。

3) 通過實驗發(fā)現(xiàn)提高了模糊核估計的精度,以及復(fù)原圖像的精度,最終平均PSNR可達到36.97,質(zhì)心提取正確率達到99.23%,質(zhì)心偏差縮小至0.018 6°。

因此該算法可進一步提高高超聲速飛行器的天文導航精度,提升大動態(tài)下光學自主導航能力。