反向追跡法獲取水下對空全景圖像

（華中光電技術(shù)研究所 武漢光電國家研究中心，湖北 武漢 430223）

引言

從水下透過波浪觀察海面及空中目標(biāo)的方法由美國Areté公司提出[1]，為水下對空探測提供了一種有效的途徑，被稱為虛擬潛望鏡技術(shù)。該技術(shù)可用于水下平臺對?？漳繕?biāo)的隱蔽探測，在軍事上有重要應(yīng)用，潛艇借此可以在不浮出海面的情況下實現(xiàn)對海面船只及空中反潛機的探測，實現(xiàn)潛艇上浮避碰及在水下安全深度對空警戒，極大地增加了潛艇的安全性。

水下對空成像屬跨介質(zhì)成像，光線經(jīng)過水-空氣界面時會發(fā)生折射，且海面常處于動態(tài)變化過程中，給水下對空成像帶來畸變和像差等不利影響。當(dāng)光線從空氣中折射入水后，會產(chǎn)生色散效應(yīng)，海水吸收和散射會使能量衰減并帶來圖像模糊，在這些因素的綜合作用下，獲取高質(zhì)量的水下對空圖像十分困難，因此，虛擬潛望鏡技術(shù)至今離實際應(yīng)用還有一定距離。然而，國內(nèi)外學(xué)者根據(jù)水下對空成像的特點，重點從某一方面進(jìn)行研究，仍然取得了重要進(jìn)展。例如，文獻(xiàn)[2]分析了接近水天線目標(biāo)的色散情況，并給出了利用楔形鏡校正色差的具體方法；文獻(xiàn)[3]提出了一種水下對空成像的光學(xué)計算模型，利用光線追跡分析波浪對光線的擾動，得到了水下像面照度的分布；文獻(xiàn)[4]中趙文強等研究了水下對空成像Snell 窗口的光學(xué)問題，解決了水下對空等效入瞳、波面曲率的等效光學(xué)特性以及天空亮度分布的計算問題；文獻(xiàn)[5]通過計算Snell圓域的擴展程度對海面波浪大小進(jìn)行估算；文獻(xiàn)[6]以一維海浪為例，通過定量計算太陽像的變化來估算海面海情，文獻(xiàn)[7]研究了從空中觀察水底目標(biāo)時產(chǎn)生的畸變，并給出了一種圖像畸變校正算法，文獻(xiàn)[8]則分析了波浪斜率對太陽反射光斑的影響，這些工作主要從定量計算的角度分析波浪對水下對空成像的影響，是本文研究的基礎(chǔ)。文獻(xiàn)[9-11]分析了動態(tài)波面對水下圖像重建的影響，提出了水下光學(xué)系統(tǒng)的設(shè)計思路，美國專利[12]提出了透過波動海面成像的方法及系統(tǒng)設(shè)計，文獻(xiàn)[13-14]主要基于統(tǒng)計原理，闡述了利用中心法對水下對空畸變圖像進(jìn)行復(fù)原處理，文獻(xiàn)[15-16]則重點分析了畸變向量場的幾何變換方法，運用波動視頻中相鄰兩幅圖像之間的相關(guān)關(guān)系對圖像進(jìn)行復(fù)原，取得了較好的效果，這些工作也為本文研究提供了有益的借鑒。

本文以海面波形對水下對空圖像的畸變影響為重點，通過MATLAB 仿真，利用主光線反向追跡獲取水下對空全景圖像，在假定海浪具有4種典型的三維表達(dá)式（海面平靜、徑向波浪、柱面波浪、二維正弦波浪）的情況下給出了仿真圖像的計算結(jié)果，為分析水下對空圖像的畸變特性提供了一種有效手段，該仿真方法可適用于具有連續(xù)一階偏導(dǎo)波形的波浪，為研究更復(fù)雜波浪下的水下對空圖像畸變特性奠定了基礎(chǔ)。

1 水下對空全景成像計算模型

1.1 水下對空全景成像系統(tǒng)組成

為簡化水下對空成像的計算過程，可建立如下理想化的成像模型。如圖1為水下對空成像系統(tǒng)示意圖，系統(tǒng)主要包括水下成像鏡頭和成像靶面，工作時假定光軸豎直向上。目標(biāo)光線從空中入射到海面上，經(jīng)水-氣界面折射后在水中向下傳輸，直到被水下鏡頭收集并折射到靶面上形成像點。從成像過程來看，應(yīng)將海面及系統(tǒng)上方海水看作水下光學(xué)系統(tǒng)的一部分，且波動的海面是一塊動態(tài)變化的透鏡，會使圖像產(chǎn)生動態(tài)的畸變。

圖1 水下對空成像系統(tǒng)示意圖Fig.1 Schematic of water-to-air imaging system

由于光線在水下會發(fā)生全反射現(xiàn)象，因此在海面平靜時，要實現(xiàn)對半球全空域的觀察，水下對空成像系統(tǒng)的視場至少為97.2°[4]。當(dāng)海面有波浪時，為了觀察到全空域，則需要適當(dāng)擴大觀察視場，在本文中，設(shè)定水下對空成像視場橫縱兩個方向的大小均為120°。

1.2 基于反向追跡法的水下對空成像計算模型

為使水下對空成像的計算簡潔，假定目標(biāo)處于無窮遠(yuǎn)處的天球之上，而水下對空成像系統(tǒng)所在區(qū)域的海平面為天球赤道，海面平靜時光軸與海平面的交點為天球中心，下面分條詳細(xì)描述該計算模型。

1） 建立坐標(biāo)系

天球坐標(biāo)系：以水下對空成像系統(tǒng)光軸與海平面交點O為天球中心，當(dāng)?shù)睾Ｆ矫鏋樘烨虺嗟?，豎直向上方向為天球北極方向。不失一般性，設(shè)正東方為赤經(jīng)0°方向（若當(dāng)?shù)貫槟媳眱蓸O，則取0°經(jīng)線方向為赤經(jīng)0°方向），逆時針為赤經(jīng)正向。設(shè)天球赤道的赤緯為0°，往天球北極方向為赤緯正向。

水下對空成像系統(tǒng)坐標(biāo)系：靶面水平放置，光軸豎直向上為z軸正向，正東方為x軸正向，正北方為y軸正向，成像系統(tǒng)的光心為坐標(biāo)原點，由此建立了天球與成像系統(tǒng)之間的坐標(biāo)聯(lián)系。

2） 主光線反向追跡

假定水下對空成像系統(tǒng)滿足小孔成像模型，且對成像靶面上的每個像素點只考慮其主光線的折射路徑，設(shè)入射光線由水下對空成像系統(tǒng)發(fā)出，其方向矢量由成像系統(tǒng)光心C指向某一像素點p，如圖2。入射光線經(jīng)海面折射后射入空中，折射光線矢量具有赤經(jīng)和赤緯兩個參量，這對應(yīng)于天球上的一點，這一點與像素點p正好滿足物像對應(yīng)關(guān)系。

圖2 靶面像素點的主光線反向追跡示意圖Fig.2 Schematic of main-ray inverse tracing of target surface pixel points

設(shè)成像靶面為正方形，水平放置，E點為靶面中心，像素點p的二維坐標(biāo)為（cx，cy），C點為成像系統(tǒng)光心，CE為焦距f，Cp則為像素點p的主光線的反向追跡，本文把它作為入射光線，Cp的方向矢量為

已知C點空間坐標(biāo)及Cp的空間方向矢量，易得Cp的空間直線方程，設(shè)海面瞬時波面方程為

求出入射光線Cp與海面的交點(x0,y0)，得到點(x0,y0)處海面的法線方向矢量

根據(jù)Snell定律

在三維空間中計算出折射光線的方向矢量，得到p點在天球的對應(yīng)點P。其中i為入射角；r為折射角；n為海水折射率。

3） 建立物像對應(yīng)關(guān)系

由上可知天球上P點與靶面像素點p為物像共軛關(guān)系，作過E、C、p3點的豎直截面進(jìn)行說明。如圖3，入射光線Cp延長與海面交于F點，出水后得到折射光線FQ，過O點作FQ的平行線，延長后與天球交于P，則P與矢量FQ對應(yīng)，即天球上點P與像素點p具有物像共軛關(guān)系。

圖3 靶面像素點與天球目標(biāo)點的對應(yīng)Fig.3 Correspondence of pixel point on target surface and target point on celestial sphere

4） 遍歷全部像素點得到對空全景圖像

若對成像靶面上全部像素點進(jìn)行遍歷，則得到一幅完整的水下對空圖像。若成像視場足夠大，則該圖像為水下對空全景圖像。

2 典型波面及場景下的水下對空全景圖像仿真

2.1 水下對空全景圖像計算流程

根據(jù)以上建立的水下對空成像模型，可以利用反向追跡法逐點重構(gòu)出整幅圖像，具體如下：

1） 設(shè)定天球目標(biāo)場景；

2） 給出海面的波面方程；

3） 設(shè)定成像系統(tǒng)探測視場、分辨率、焦距、工作深度等成像參數(shù)；

4） 對成像靶面所有像素點進(jìn)行遍歷，每個像素點需完成5）～9）之間的操作；

5） 求出像素點對應(yīng)的入射光線的直線方程；

6） 求出入射光線與波面的交點(x0,y0)；

7） 求出過點(x0,y0)的波面法線方程；

8） 根據(jù)Snell定律求出折射光線方向矢量；

9） 求出與折射光線對應(yīng)的天球目標(biāo)點，該點的灰度值即為對應(yīng)像素點的灰度值；

10） 輸出水下對空全景圖像。

入射光線與波面的交點可用二分法進(jìn)行求解，如圖4所示，具體過程：設(shè)入射光線的兩個端點C和D處于波面異側(cè)，其坐標(biāo)分別為和則與異號，不妨設(shè)V1＞0，V2＜0。設(shè)并設(shè)F點的坐標(biāo)為計算若V0＞0，則將C點坐標(biāo)更新為若V0＜0，則將D點坐標(biāo)更新為繼續(xù)對端點C和D進(jìn)行以上迭代計算，直至V0的絕對值小于給定的小正數(shù)為止，此時即為直線與波面的交點近似解。

圖4 二分法求解直線與曲面交點Fig.4 Solving intersection of straight line and curved surface by bisection method

根據(jù)Snell定律計算折射光線矢量需要在三維空間中求解，如圖5所示，過程如下：設(shè)入射光線、折射光線及法線的方向矢量分別為R=（rx，ry，rz）和N=（nx，ny，nz），由于事先已求得F點的坐標(biāo)故I和N均為已知量，設(shè)入射角為i，折射角為r，則有i=arccos(I·N)，r=arccos(R·N)，由于R為單位矢量，故|R|=1，由Snell定律，有n·sin(i)=sin(r)，由于I、R和N3個矢量必須共面，故I×N·R=0，由于入射光線與折射光線必須位于在法線的兩側(cè)，故(I×N)·(R×N)≥0。根據(jù)以上條件，可以聯(lián)立解出折射光線的方向矢量R，即（rx，ry，rz）。

圖5 三維空間中求解折射光線矢量Fig.5 Solving direction vector of refracted ray in 3D space

2.2 邊界條件設(shè)定

對天球目標(biāo)場景、成像參數(shù)及波面方程進(jìn)行如下設(shè)定。

1） 天球目標(biāo)場景

設(shè)定4個天球目標(biāo)及背景如表1，即灰度值為255、角直徑為10°的4個圓形亮目標(biāo)懸掛在灰度值為128的均勻亮度天空中，位于正東方，高度角分別為10°、30°、60°、90°。

表1 天球目標(biāo)場景Table1 Celestial sphere target scenes

2） 成像系統(tǒng)參數(shù)

海面平靜時Snell圓錐的張角為97.2°，海面有波浪時水下視場會增大，故設(shè)定成像系統(tǒng)的探測視場大小為120°×120°，靶面分辨率為1024×1024像素，鏡頭焦距為3 mm，工作水深為50 m。在以下仿真計算中，假定成像系統(tǒng)沒有像差和畸變，滿足小孔成像原理。

3） 波面方程

設(shè)定波面方程為4種形式：海面平靜；徑向波浪；柱面波浪；二維正弦波浪。具體在下文詳細(xì)描述。

2.3 水下對空全景圖像計算結(jié)果

2.3.1 海面平靜下的計算結(jié)果

海面平靜時的波面方程為

此時半球天空在靶面上的像正好是一個圓，且其對應(yīng)的水下視場為97.2°，由于水下對空全景成像系統(tǒng)的視場為120°，因此在此Snell圓區(qū)域之外由于發(fā)生全反射現(xiàn)象而為純黑色。如圖6，灰色大圓為Snell圓，即半球天空所成像的區(qū)域，天球上的4個圓形目標(biāo)在靶面上所成的像從右往左依次對應(yīng)1、2、3和4 號目標(biāo)，即圖中的白色區(qū)域，赤緯角小的1 號目標(biāo)像靠近Snell圓邊緣，而正當(dāng)空的4 號目標(biāo)像在Snell圓中心。

圖6 海面平靜時水下對空全景圖像Fig.6 Water-to-air panoramic image under calm sea surface

2.3.2 徑向波浪下的計算結(jié)果

海面為徑向波浪時，其波面方程為

圖7是該波面方程的直觀外形示意圖。

圖7 徑向波浪波形示意圖Fig.7 Schematic of radial waveform

由圖可見徑向波浪是以一點為中心的旋轉(zhuǎn)對稱的同心圓，它由平面余弦曲線繞z軸旋轉(zhuǎn)一周生成，如（7）式所示：

水下對空全景成像系統(tǒng)處于徑向波浪中心正下方一定深度處，通過以下MATLAB 代碼來確定徑向波浪的波形參數(shù)：

其中：A為波浪振幅；w為波浪角頻率；T為波長；steep為波陡，即波高與波長的比值；depth表示成像系統(tǒng)光心C的深度；DIA為觀察視場在海面的寬度；N表示視場寬度范圍內(nèi)波浪的周期數(shù)；波高為振幅A的兩倍。在以上成像參數(shù)條件下，通過MATLAB 仿真得到圖8所示全景圖像。

圖8 徑向波浪下水下對空全景圖像Fig.8 Water-to-air panoramic image under radial wave

如圖8，由于波浪的旋轉(zhuǎn)對稱特征，天球區(qū)域的圖像仍然是灰色大圓，直觀看來圖8與圖6相似，其中目標(biāo)1到目標(biāo)4的白色像次序也一致，但仔細(xì)對照可以發(fā)現(xiàn)圖8的Snell圓直徑比圖6大，4個目標(biāo)像的形態(tài)也不相同，其位置也存在一定的偏移，特別地，兩圖中目標(biāo)4的圖像直徑不同。

2.3.3 柱面波浪下的計算結(jié)果

海面為柱面波浪時，其波面方程為

圖9是該波面方程的直觀外形示意圖。

圖9 柱面波浪波形示意圖Fig.9 Schematic of cylindrical waveform

波形參數(shù)取值與2.3.2節(jié)一致，通過MATLAB仿真得到圖10所示水下對空全景圖像，可以明顯看出Snell圓的邊緣已經(jīng)發(fā)生明顯的變形，4個天球目標(biāo)像的位置次序雖一致，但其形態(tài)和位置都與圖6有一定變化，特別地，最左邊目標(biāo)4的白色像不再是圓。

圖10 柱面波浪下水下對空全景圖像Fig.10 Water-to-air panoramic image under cylindrical wave

2.3.4 二維正弦波浪下的計算結(jié)果

海面為二維正弦波浪時，其波面方程為

圖11是該波面方程的直觀外形示意圖。

圖11 二維正弦波浪波形示意圖Fig.11 Schematic of 2D sinusoidal waveform

波形參數(shù)取值與2.3.2節(jié)一致，且取

由此通過MATLAB 仿真得到圖12所示水下對空全景圖像?？梢钥吹?，Snell圓的邊緣發(fā)生了更復(fù)雜的變形，4個天球目標(biāo)的位置次序雖一致，但其形態(tài)和位置都與圖6有一定變化，其中目標(biāo)1的圖像的形變特別大，目標(biāo)4的圖像也不是圓。

圖12 二維正弦波浪下水下對空全景圖像Fig.12 Water-to-air panoramic image under 2D sinusoidal wave

2.4 波面參數(shù)對水下對空全景圖像的影響

上文分別給出了海面平靜、徑向波浪、柱面波浪及二維正弦波浪下的水下對空全景仿真圖像，仿真中波陡steep取了較小值，視場內(nèi)波浪的周期數(shù)N也較小，這種情況下4個天球目標(biāo)的像畸變不是特別嚴(yán)重。事實上，在較大波浪下，每個圓形目標(biāo)的圖像極易分裂成多個，為避免碎化后的目標(biāo)像混疊難以分辨，下文僅對天球目標(biāo)2 進(jìn)行仿真。由于當(dāng)波陡steep和視場內(nèi)波浪周期數(shù)N這兩個參數(shù)確定后，水下對空全景圖像也就確定，與成像系統(tǒng)的深度depth無關(guān)，因此，下面只考察steep和N的變化對全景圖像的影響。

2.4.1 徑向波浪波面參數(shù)變化的影響

嚴(yán)格來說，應(yīng)對steep和N分別取一系列值來仿真水下對空圖像，為了直觀和簡潔起見，下面各節(jié)均在兩種情況下給出仿真結(jié)果：一是steep取固定值，N取4個不同值；二是N取固定值，steep取4個不同值。具體取值以能直觀展示圖像的變化趨勢來確定。

圖13為steep=0.05時目標(biāo)2的仿真圖像。由圖可見，隨著N值增大，天球成像區(qū)域大小變化較小，但目標(biāo)2的圖像逐漸碎化。

圖13 steep=0.05時目標(biāo)2的仿真圖像Fig.13 Simulated images of target No.2 when steepis 0.05

圖14為N=10時目標(biāo)2的仿真圖像。由圖可見，隨著steep值增大，目標(biāo)2的圖像逐漸碎化，同時Snell圓域也向外擴展。

圖14 N=10時目標(biāo)2的仿真圖像Fig.14 Simulated images of target No.2 when Nis 10

2.4.2 柱面波浪波面參數(shù)變化的影響

圖15為steep=0.05時目標(biāo)2的仿真圖像。由圖可見，隨著N值增大，天球成像區(qū)域大小變化較小，但目標(biāo)2的圖像逐漸碎化。

圖15 steep=0.05時目標(biāo)2的仿真圖像Fig.15 Simulated images of target No.2 when steep=0.05

圖16為N=10時目標(biāo)2的仿真圖像。由圖可見，隨著steep值增大，目標(biāo)2的圖像逐漸碎化，同時Snell圓域也向外擴展。

圖16 N=10時目標(biāo)2的仿真圖像Fig.16 Simulated images of target No.2 when N=10

2.4.3 二維正弦波浪波面參數(shù)變化的影響

圖17為steep=0.05時目標(biāo)2的仿真圖像。由圖可見，隨著N值增大，天球成像區(qū)域大小變化較小，但目標(biāo)2的圖像逐漸碎化。

圖17 steep=0.05時目標(biāo)2的仿真圖像Fig.17 Simulated images of target No.2 when steep=0.05

圖18為N=10時目標(biāo)2的仿真圖像。由圖可見，隨著steep值增大，目標(biāo)2的圖像逐漸碎化，同時Snell圓域也向外擴展。

圖18 N=10時目標(biāo)2的仿真圖像Fig.18 Simulated images of target No.2 when N=10

2.4.4 波浪參數(shù)對水下對空成像的定量影響

由以上3小節(jié)給出的仿真結(jié)果可以看出，波浪周期數(shù)N顯著影響圖像的碎化程度，但對整個成像區(qū)域的擴展程度影響很小，而波陡steep既能影響圖像的碎化程度，同時又顯著影響圖像的擴展程度，這是波浪參數(shù)對水下對空成像的定性影響。然而，根據(jù)本文提出的計算方法，只要給定波浪參數(shù)，仿真程序?qū)嶋H上已經(jīng)輸出了對空圖像定量結(jié)果，可以從仿真結(jié)果確切計算出成像區(qū)域的擴展程度與波浪參數(shù)的對應(yīng)關(guān)系，而且當(dāng)給出波浪的運動方程時，也完全可以仿真出對空圖像隨波浪變化的視頻。

3 結(jié)論

本文利用主光線反向追跡法仿真水下對空全景圖像，給出了根據(jù)Snell定律求解水下對空全景圖像的一般過程。在海面平靜、徑向波浪、柱面波浪和二維正弦波浪下得到了典型天球目標(biāo)及背景的圖像，同時，通過改變波陡和視場內(nèi)波浪周期數(shù)兩個波浪參數(shù)，直觀地表現(xiàn)了波浪的變化對水下對空全景圖像的影響。根據(jù)建立的成像模型，在水下對空全景圖像、波浪方程和天球目標(biāo)場景三者之間，若已知其中兩個，則可求出另一個，這為分析水下對空成像的畸變特性提供了定量的方法。該方法適用于分析透過具有一階可導(dǎo)波浪的水下對空全景圖像的畸變特征，為研究更復(fù)雜波浪下的水下對空圖像畸變特征奠定了基礎(chǔ)。

然而，波浪對水下對空成像的影響十分復(fù)雜，除了利用主光線反向追跡所描述的圖像畸變之外，在實際成像過程中還存在離焦、色散、偏振、散射、吸收、成像噪聲等多方面的不利影響，太陽不可避免地進(jìn)入視場后也會對成像產(chǎn)生干擾，這些因素將嚴(yán)重降低水下對空成像質(zhì)量，增加目標(biāo)探測和識別的難度。因此，還需對這些問題進(jìn)行深入研究。