何大華，李陽(yáng)陽(yáng)，周少杰

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

1 引 言

水下光電成像是利用可見光對(duì)水下目標(biāo)進(jìn)行成像的技術(shù)[1-12]，在國(guó)民經(jīng)濟(jì)乃至軍事上有廣泛的應(yīng)用前景。由于水體對(duì)光線存在散射，無(wú)論是人工光源還是自然光源照射到水體時(shí)，都會(huì)在整個(gè)水域中形成水下光場(chǎng)[13-14]，其中處于目標(biāo)和水下光電成像系統(tǒng)之間的水體散射將在成像系統(tǒng)靶面上形成背景干擾，降低目標(biāo)對(duì)比度和成像質(zhì)量。水下光電成像可分為主動(dòng)和被動(dòng)兩種成像方式，水下主動(dòng)光電成像利用系統(tǒng)自帶的人工光源對(duì)水下目標(biāo)進(jìn)行照明，可以是鹵素?zé)簟ED燈或者藍(lán)綠激光器等，水下被動(dòng)光電成像利用自然光進(jìn)行水下照明，包括太陽(yáng)或者天空背景光等。可見光在水下傳輸時(shí)由于受到水體的強(qiáng)烈吸收和散射，其能量隨距離增大按指數(shù)規(guī)律衰減，使得無(wú)論是被動(dòng)還是主動(dòng)方式，甚至近年來(lái)取得重要進(jìn)展的水下激光成像技術(shù)[2-5]，其水下作用距離最多不過(guò)幾十米。在另一個(gè)重要的水下應(yīng)用領(lǐng)域，水下激光通信[15-16]中，水體的吸收和散射現(xiàn)象對(duì)通信距離和誤碼率同樣會(huì)產(chǎn)生嚴(yán)重影響。

對(duì)水體散射等光學(xué)現(xiàn)象的研究可追溯到20世紀(jì)60年代，Duntley深入研究了光線被海水吸收和散射的基本性質(zhì)[17]；GROSSO在理論建模和實(shí)驗(yàn)對(duì)比的基礎(chǔ)上給出了體散射函數(shù)（VSF）與調(diào)制傳遞函數(shù)（MTF）的關(guān)系[18]；MERTENS給出水體點(diǎn)擴(kuò)展函數(shù)（PSF）和光束擴(kuò)展函數(shù)（BSF）的定義，并研究了二者與水體的固有光學(xué)參數(shù)（IOP）之間的關(guān)系[19]；JAFFE建立了水下成像的最優(yōu)模型，發(fā)現(xiàn)調(diào)整光源和成像系統(tǒng)的相對(duì)位置能提升圖像對(duì)比度[20]；HOU建立了簡(jiǎn)單的水下成像模型，并研究了PSF在自然水體中成像的有效性，利用PSF對(duì)圖像進(jìn)行去卷積復(fù)原，獲得了較好效果[21-22]；VOSS給出了一個(gè)估算海水PSF的經(jīng)驗(yàn)公式，可以用來(lái)預(yù)測(cè)任意長(zhǎng)度水體的PSF[23]。國(guó)內(nèi)方面，陳養(yǎng)渭對(duì)水下激光散射場(chǎng)進(jìn)行了實(shí)驗(yàn)研究，并建立了散射橢球模型[24]；閆旭光利用Monte Carlo方法分析了激光前向散射的時(shí)間和空間特性[25]；宋慶君等研究了黃海、東海海區(qū)水體的后向散射系數(shù)與總懸浮物濃度的關(guān)系，為水色遙感提供了支撐[26]；李彩等對(duì)水體體散射函數(shù)的測(cè)量技術(shù)進(jìn)展進(jìn)行了綜述[27]；許廷發(fā)等利用連續(xù)幀圖像信息估算PSF，結(jié)合凸集投影法進(jìn)行超分辨率成像重建[28]；褚金奎等引入偏振技術(shù)，提出了一種抑制水下圖像后向散射的方法[29]。

本文以水體吸收和散射的基本過(guò)程為基礎(chǔ)，在穩(wěn)定光源照射以及體散射函數(shù)為球形對(duì)稱的情況下，定量求解光線經(jīng)水體多次散射后形成的水下光場(chǎng)分布。求解水下光場(chǎng)的過(guò)程闡明了水體對(duì)光線吸收和散射的微觀機(jī)制，精確的水下光場(chǎng)數(shù)據(jù)為嚴(yán)格推導(dǎo)水體PSF和MTF提供了關(guān)鍵信息。

2 水下光場(chǎng)計(jì)算模型

2.1 基本概念

一般認(rèn)為，準(zhǔn)直窄光束在水下傳輸時(shí)能量衰減遵循朗伯-比爾定律[13]，即

其中P0為準(zhǔn)直窄光束的初始功率，P為準(zhǔn)直窄光束在水下傳輸距離r之后的功率，c為水體衰減系數(shù)，其倒數(shù)1/c稱為衰減長(zhǎng)度（Attenuation Length，AL），a為吸收系數(shù)，b為散射系數(shù)。

對(duì)于水下點(diǎn)光源，設(shè)光強(qiáng)為I0，則在水下傳輸距離r后未被吸收和散射的能量形成的徑向照度為

多次散射：光束或光子在水下傳輸過(guò)程中被水分子或水中雜質(zhì)散射從而多次改變傳輸方向的現(xiàn)象。

球面照度：空間中某點(diǎn)的球面照度為全空間范圍內(nèi)入射到該點(diǎn)處一個(gè)小球內(nèi)的光通量總和與該小球的表面積之比，可表示為4π立體角范圍內(nèi)對(duì)亮度的積分。球面照度也稱為標(biāo)量照度。

水下初始光場(chǎng)：在指定的水域中，光源發(fā)射的光能量進(jìn)入水體傳輸過(guò)程中未被水體吸收和散射的部分在水下形成的球面照度分布。

水下散射光場(chǎng)：在指定的水域中，水體散射光形成的球面照度分布。

水下光場(chǎng)：在指定的水域中，水下初始光場(chǎng)與水下散射光場(chǎng)之和。

2.2 水下光場(chǎng)的形成過(guò)程

暫不考慮水氣界面、水下物體等邊界條件以及瞬態(tài)光源或光源功率波動(dòng)的影響。假定在三維直角坐標(biāo)系oxyz中充滿了水體，衰減系數(shù)、吸收系數(shù)與散射系數(shù)分別為c、a和b，其中b=4πβ，式中 β為體散射函數(shù)[13]。下面以單點(diǎn)光源為例說(shuō)明水下光場(chǎng)的形成過(guò)程，設(shè)在(x0,y0,z0)處有光強(qiáng)為I0的穩(wěn)恒點(diǎn)光源S，則水下光場(chǎng)按如下步驟形成：

Ⅰ.點(diǎn)光源發(fā)出的光線經(jīng)擴(kuò)散和水體衰減后形成水下初始光場(chǎng)Et0。

Ⅱ.水下初始光場(chǎng)Et0使整個(gè)水域產(chǎn)生散射，形成一次水下散射光場(chǎng)Es1，疊加在Et0上形成一次水下光場(chǎng)Et1=Et0+Es1。

Ⅲ.一次水下光場(chǎng)Et1使整個(gè)水域產(chǎn)生散射，形成二次水下散射光場(chǎng)Es2，疊加在Et0上形成二次水下光場(chǎng)Et2=Et0+Es2。

Ⅳ.i次水下光場(chǎng)Eti使整個(gè)水域產(chǎn)生散射，形成i+1次水下散射光場(chǎng)Es(i+1)，疊加在Et0上形成i+1次水下光場(chǎng)Et(i+1)=Et0+Es(i+1)。

Ⅴ.以此類推，將形成水下初始光場(chǎng)Et0及水下光場(chǎng)序列Eti,i=1,2,3,······。

Ⅵ.Eti將收斂于水下光場(chǎng)Et。

根據(jù)以上思路，下面給出一般意義上的水下光場(chǎng)求解方程式。

如圖1，點(diǎn)光源S發(fā)出的光在水下傳輸過(guò)程中會(huì)進(jìn)行球面擴(kuò)散，還會(huì)被水體吸收、散射，一部分光線繼續(xù)沿原來(lái)的路徑傳播，在P點(diǎn)處形成的水下初始光場(chǎng)為

由于受到光源照射，水體也成為散射光源，M點(diǎn)處體積元dudvdw的散射會(huì)在P點(diǎn)處產(chǎn)生球面照度，對(duì)所有體積元散射產(chǎn)生的球面照度進(jìn)行疊加可得到P點(diǎn)處的水下散射光場(chǎng)Es(x,y,z)。

圖1 水下光場(chǎng)示意圖Fig.1 Diagram of the underwater light field

設(shè)水下光場(chǎng)為Et(x,y,z)，則有

P點(diǎn)處的水下光場(chǎng)等于水下初始光場(chǎng)與水下散射光場(chǎng)之和，因而有

式（5）即為水下光場(chǎng)Et(x,y,z)的求解方程式，它屬于第二類Fredholm積分方程，且未知函數(shù)有3個(gè)變量。在滿足一定條件時(shí)，可利用行列式求解得到近似解析解[30]。本文從工程應(yīng)用出發(fā)，給出該類方程的數(shù)值迭代解法，具體計(jì)算方法在3.1節(jié)中給出，僅考慮水下光場(chǎng)含一個(gè)變量的情形，對(duì)多個(gè)變量的情形可依此類推。

3 典型水下光場(chǎng)的計(jì)算方法

3.1 太陽(yáng)照射下水下光場(chǎng)的計(jì)算方法

2.2節(jié)給出了理想情況下水下光場(chǎng)的計(jì)算方法，實(shí)際中總會(huì)存在一些邊界條件，下面給出太陽(yáng)照射下足夠大平靜水面下的水下光場(chǎng)。

僅考慮太陽(yáng)照度，忽略天空背景亮度的影響，太陽(yáng)光以一定角度平行入射進(jìn)入水中，并發(fā)生反射和折射。折射光線經(jīng)吸收和散射后，剩下部分繼續(xù)向前傳播，隨深度增加能量越來(lái)越小。散射光線同樣會(huì)被吸收、散射或者繼續(xù)傳播，在水下形成多次散射，進(jìn)而形成穩(wěn)定的水下光場(chǎng)。

如圖2所示，從對(duì)稱性考慮，太陽(yáng)照射下的水下光場(chǎng)的數(shù)值僅僅與深度有關(guān)，而與水平坐標(biāo)無(wú)關(guān)，因此，水下光場(chǎng)函數(shù)可用深度h為自變量的一元函數(shù)Et(h)表示。

圖2 太陽(yáng)照射下的水下光場(chǎng)Fig.2 Underwater light field from the sun

設(shè)太陽(yáng)光的入射角為α，水體的折射率為nw，由折射定律可求出折射角γ，再由Fresnel公式[31]可求出入水時(shí)的反射率Rα和透射率Tα。

設(shè)與入射太陽(yáng)光垂直方向上的照度為Esun，則水下深度h處的水下初始光場(chǎng)為

以圖3來(lái)說(shuō)明水下深度h處P點(diǎn)的水下光場(chǎng)Et(h)的形成過(guò)程，Et(h)=Ed(h)+Es(h)，Ed(h)為水下初始光場(chǎng)，表達(dá)式為式（6），Es(h)為水下散射光場(chǎng)，在數(shù)值上等于水下所有微體積元散射在P點(diǎn)形成的球面照度之和。為此，建立直角坐標(biāo)系，以P點(diǎn)正上方與水面相交點(diǎn)O為坐標(biāo)原點(diǎn)，水平向右為x軸正向，垂直于紙面向外為y軸正向，豎直向下為z軸正向，則P點(diǎn)的坐標(biāo)為P(0,0,h)。

圖3 水下球面照度形成原理Fig.3 Formation of underwater spherical illuminance

PO與深度為z的水平面相交于C，以C為圓心在水平面內(nèi)畫半徑為r的圓C，MN為直徑，其中M與M'、N與N'分別關(guān)于水面對(duì)稱，PM'交水面于A點(diǎn)，PN'交水面于B點(diǎn)，以M'N'為直徑的水平圓的圓心為D，由圖3可知，在xOy平面內(nèi)，以MN為直徑的圓周上所有點(diǎn)到P點(diǎn)距離相等，記為。

圓C上的微體積元dV=2πr·drdz，體積元dV在P點(diǎn)形成的球面照度dEs1為

還應(yīng)考慮圓C上的微體積元散射光經(jīng)水面下表面反射后對(duì)P點(diǎn)形成的球面照度。M點(diǎn)發(fā)出的光線經(jīng)A點(diǎn)反射可到達(dá)P，可等效為從M'發(fā)出的光線直接到達(dá)P。同理，N點(diǎn)發(fā)出的光線經(jīng)B點(diǎn)反射也可到達(dá)P，可等效為從N'發(fā)出的光線直接到達(dá)P，故這組反射光線可等效從圓D發(fā)出，并經(jīng)過(guò)了長(zhǎng)度為PM'的水程到達(dá)P點(diǎn)，記為。

注意到光線MA在A點(diǎn)反射時(shí)有能量損失，由入射角α根據(jù)Fresnel公式可計(jì)算出反射率Rα。由此可知圓C上的微體積元經(jīng)水面下表面反射在P點(diǎn)形成的球面照度dEs2為

綜合以上結(jié)果，圓C上的微體積元在P點(diǎn)形成的球面照度dEs為

其中Rα、l、m均為h、r和z的函數(shù)。

由此可知P點(diǎn)的水下散射光場(chǎng)Es(h)為

故水下光場(chǎng)Et(h)為

式（11）為第二類Fredholm積分方程，對(duì)式(11)進(jìn)行代換簡(jiǎn)化得到

其中令深度h及積分變量r和z的細(xì)分間隔均為Δ，并對(duì)式（12）進(jìn)行數(shù)值化得

其中，k=1,2,3,······。

由于水體對(duì)光線的衰減嚴(yán)重，水體散射對(duì)距離超過(guò)10AL的區(qū)域影響可以忽略，因此在計(jì)算時(shí)，可適當(dāng)限定計(jì)算范圍，基本不影響計(jì)算精度，例如計(jì)算10AL深度內(nèi)的水下光場(chǎng)時(shí)，可令k=1,2,3,···,N，并滿足NΔ＞20AL。據(jù)此，將式（14）改為以下迭代形式

式（15）的迭代求解流程圖如圖4所示。

圖4 水下光場(chǎng)迭代求解流程圖Fig.4 Iterative solution flow chart of the underwater light field

迭代6次后Et(kΔ)已基本收斂，可輸出結(jié)果。下面利用Matlab仿真程序迭代計(jì)算太陽(yáng)照射下的水下光場(chǎng)，選取晴朗天氣、潔凈海水條件下的水下光場(chǎng)進(jìn)行計(jì)算，可以從仿真曲線直觀看出這一典型條件下的水下光場(chǎng)隨深度的變化情況。

仿真程序的相關(guān)參數(shù)取值如下：太陽(yáng)光照度：Esun=105lx；太陽(yáng)光入射角：α=45°；水體折射率：nw=1.33；衰減系數(shù)：c=0.1/m；吸收系數(shù)：a=0.05/m；散射系數(shù)：b=0.05/m；體散射函數(shù)：β=1/80π/(m·sr)；Et(h)深度精度：Δ=0.2m；

Δ取值的依據(jù)是，當(dāng)Δ=0.02 AL時(shí)，水下光場(chǎng)精度足夠，此例中AL=10m，因而Δ=0.2m。

圖5是水深100 m（10AL ）以內(nèi)仿真程序輸出的水下光場(chǎng)迭代曲線。

圖5 太陽(yáng)照射下水下光場(chǎng)迭代曲線Fig.5 Iteration curves of the underwater light field from the sun

由于Et(h)=Ed(h)+Es(h)，為了直觀表達(dá)，不直接顯示Et(h)，而顯示分量Ed(h)和Es(h)，其中Ed(h)是定值，在i從0到5這6次迭代過(guò)程中保持不變，而Es(h)在迭代過(guò)程中逐漸收斂。圖6最上面的曲線為Ed(h)，下面6條曲線為Es(h)的6次迭代變化（包括Es(h)=0），順序?yàn)閺南轮辽?，從圖6中可以看出6次迭代后已基本收斂。

3.2 均勻亮度天空背景下的水下初始光場(chǎng)

由3.1節(jié)對(duì)水下光場(chǎng)的求解過(guò)程可知，當(dāng)光源選定為均勻亮度天空背景時(shí)，計(jì)算過(guò)程中唯一差別是求解水下初始光場(chǎng)Ed(h)，在不同的光源照射下，只要求出了Ed(h)，則可利用上述方法得到迭代方程式對(duì)水下光場(chǎng)Et(h)進(jìn)行求解。下面推導(dǎo)均勻亮度天空背景下的Ed(h)表達(dá)式，并給出Matlab程序的仿真結(jié)果。

由于天空背景為大型面光源，各部分光線入水時(shí)的入射角 α各不相同，因此，求解Ed(h)時(shí)需對(duì)上半球天空背景進(jìn)行積分。

如圖6所示，天頂角為α的微球帶對(duì)應(yīng)的立體角為dΩ=2πsinαdα。

圖6 天空背景在水面形成照度Fig.6 Illuminance on the surface with the sky as the background

設(shè)均勻天空的亮度為L(zhǎng)sky，對(duì)水下初始光場(chǎng)求解時(shí)，可將微球帶發(fā)出的光線等效為入射角為α的一組平行光，而且垂直于入射方向上的照度為dE=2Lskyπsinαdα，則可得均勻亮度天空背景下水下初始光場(chǎng)為

編寫Matlab仿真程序計(jì)算Et(h)，假設(shè)天空亮度為L(zhǎng)sky=5000cd/m2，水體光學(xué)參數(shù)及Et(h)的深度精度同3.1節(jié)，則Et(h)在深度10AL之內(nèi)的迭代曲線如圖7所示。其中各曲線的意義與圖5相同。

圖7 天空背景照射下水下光場(chǎng)迭代曲線Fig.7 Iteration curves of the underwater light field with the sky as the background

3.3 點(diǎn)光源形成的水下初始光場(chǎng)

3.3.1 點(diǎn)光源在水面以下

設(shè)水深h0處有均勻發(fā)光點(diǎn)光源S，光強(qiáng)為I0，水體衰減系數(shù)、吸收系數(shù)及散射系數(shù)同上，計(jì)算水面平靜時(shí)光源S形成的水下初始光場(chǎng)。

如圖8所示，由于對(duì)稱性，僅分析通過(guò)S點(diǎn)的一個(gè)豎直平面內(nèi)的情形，對(duì)于水下任意點(diǎn)P，設(shè)其深度為h，離S點(diǎn)的水平距離為x，則P點(diǎn)處的水下初始光場(chǎng)記為Ed(h,x)。

圖8 水下點(diǎn)光源形成的水下初始光場(chǎng)Fig.8 Underwater initial light field by an underwater point light source

光源S直射到P點(diǎn)處的球面照度為Ed1(h,x)，光源S點(diǎn)發(fā)出的光線經(jīng)水面下表面反射后照射到P點(diǎn)處的球面照度為Ed2(h,x)，則

S點(diǎn)發(fā)出的光線經(jīng)水面下表面反射的光線可以等效為S點(diǎn)關(guān)于水面的對(duì)稱點(diǎn)S′點(diǎn)發(fā)出，經(jīng)過(guò)了水程S′P的衰減以及反射率Rα的修正，則

3.3.2 點(diǎn)光源在水面以上

設(shè)水面以上高h(yuǎn)0處有均勻發(fā)光點(diǎn)光源S，光強(qiáng)為I0，水體衰減系數(shù)、吸收系數(shù)及散射系數(shù)同上，下面計(jì)算水面平靜時(shí)光源S形成的水下初始光場(chǎng)。

如圖9所示，通過(guò)S作豎直線與水面交于O點(diǎn)，顯然，水下初始光場(chǎng)關(guān)于直線SO旋轉(zhuǎn)對(duì)稱，因此只分析通過(guò)S點(diǎn)的豎直平面內(nèi)的情形。對(duì)于水下任意點(diǎn)P，設(shè)其深度為h，離O點(diǎn)的水平距離為x，則P點(diǎn)處的水下初始光場(chǎng)值記為Ed(h,x)。SO與深度為h的水平面交于M，則以M為圓心，x為半徑的水平圓上的水下初始光場(chǎng)值均為Ed(h,x)。

光源S發(fā)出的光線中入射角為 α的光線形成一個(gè)錐角為2α的圓錐面，其水面入射點(diǎn)集合是圓，圓心為O，直徑為AB，設(shè)折射角為γ，透射率為Tα。

圖9 水上點(diǎn)光源形成的水下初始光場(chǎng)Fig.9 Underwater initial light field by an overwater point light source

錐角為2α的圓錐立體角為2π(1-cosα)，考慮圓錐2(α+dα)與2α之間的光線折射入水能量傳輸情況，其入射角均為α，透射率均為Tα，立體角為2πsinα·dα，對(duì)應(yīng)能量為Φ0=2πI0sinα·dα，該能量經(jīng)水面反射和水體衰減后到達(dá)深度h處時(shí)變?yōu)?/p>

由圖9可知x=h0·tanα+h·tanγ，在深度h處，折射光線是一個(gè)圓環(huán)，內(nèi)半徑為x，外半徑為x+dx，其面積ds=2πx·dx。由于P處的球面照度為Ed(h,x)，則Φh又可表示為

比較式（20）和式（21）并化簡(jiǎn)得

3.4 小結(jié)

以上3小節(jié)分別給出了平行光、天球光、水下和水上點(diǎn)光源形成的水下光場(chǎng)的計(jì)算方法，并考慮到水面下表面反射的影響。假定光線的傳輸及散射現(xiàn)象具有獨(dú)立性和疊加性，則在以上各類光源的任意組合下求解水下光場(chǎng)也是可行的。在復(fù)色光入射的情況下，也可以對(duì)各個(gè)波段分別計(jì)算然后進(jìn)行積分得到結(jié)果。這些光源的選取具有一定的代表性，因此，利用該計(jì)算模型可解決一般意義上的水下光場(chǎng)分布問(wèn)題。

在自帶光源的主動(dòng)水下光電成像應(yīng)用中，為了提高功率密度和減小后向散射，光源的照明區(qū)域往往局限在有限的角度范圍內(nèi)，這時(shí)水下光場(chǎng)也可以通過(guò)上述方法求解，因此該計(jì)算模型可擴(kuò)展應(yīng)用到邊界條件更為復(fù)雜一些的情形。

已知水下光場(chǎng)分布，則易求得水體的表觀光學(xué)參數(shù)[13]，因此，該計(jì)算模型還可直接用于求解水體的表觀光學(xué)參數(shù)（AOP）。

4 結(jié) 論

水下光電成像技術(shù)目前得到了廣泛的關(guān)注，在未來(lái)具有可觀的經(jīng)濟(jì)和軍事價(jià)值。近年來(lái)軍事方面的一個(gè)緊迫需求是水下對(duì)空跨介質(zhì)成像[32]，由此水下航行器得以在安全深度下對(duì)空中威脅目標(biāo)進(jìn)行預(yù)警。其主要難點(diǎn)在于：一是海面波浪引起的圖像碎化問(wèn)題[33-34]；二是水體對(duì)光線的吸收和散射效應(yīng)引起的圖像模糊。這需要研究水下光場(chǎng)分布，改善圖像質(zhì)量，以提升水下航行器的預(yù)警深度。

從基本的吸收和散射模型出發(fā)，考慮到有無(wú)邊界條件，從理論上嚴(yán)格推導(dǎo)了幾種典型光源照射下的水下光場(chǎng)分布，這對(duì)于分析水體PSF以及MTF具有重要意義，也為研究水下光電圖像的退化模型奠定了基礎(chǔ)。雖然只給出了水體VSF為球形對(duì)稱且水面平靜情形下水下光場(chǎng)的求解方法，但對(duì)于VSF具有任意形式、水面波形具有一階連續(xù)可導(dǎo)的情形也能求解，這為一般意義上水下光場(chǎng)問(wèn)題的定量求解提供了借鑒。

然而，在實(shí)際問(wèn)題中，光源特性并不穩(wěn)定，水體散射特性在空間上也不均勻，且存在湍流時(shí)變擾動(dòng)[35-37]，故水下光場(chǎng)是動(dòng)態(tài)變化的，且水下各種噪聲也對(duì)計(jì)算結(jié)果有較大影響。因此，在具體的水下光場(chǎng)求解過(guò)程中，需要結(jié)合估計(jì)、現(xiàn)場(chǎng)測(cè)量、大數(shù)據(jù)處理等經(jīng)驗(yàn)方法進(jìn)行分析，而將本文的方法作為必要的參考。