何大華,李陽(yáng)陽(yáng),周少杰
(華中光電技術(shù)研究所—武漢光電國(guó)家研究中心, 湖北 武漢 430223)
水下光電成像是利用可見(jiàn)光對(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é)?、LED燈或者藍(lán)綠激光器等,水下被動(dòng)光電成像利用自然光進(jìn)行水下照明,包括太陽(yáng)或者天空背景光等??梢?jiàn)光在水下傳輸時(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)鍵信息。
一般認(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)之和。
暫不考慮水氣界面、水下物體等邊界條件以及瞬態(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)。
Ⅴ.以此類(lèi)推,將形成水下初始光場(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)的求解方程式,它屬于第二類(lèi)Fredholm積分方程,且未知函數(shù)有3個(gè)變量。在滿足一定條件時(shí),可利用行列式求解得到近似解析解[30]。本文從工程應(yīng)用出發(fā),給出該類(lèi)方程的數(shù)值迭代解法,具體計(jì)算方法在3.1節(jié)中給出,僅考慮水下光場(chǎng)含一個(gè)變量的情形,對(duì)多個(gè)變量的情形可依此類(lèi)推。
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)畫(huà)半徑為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)為第二類(lèi)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.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)為
編寫(xiě)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.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(mén)α。
圖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(mén)α,立體角為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小節(jié)分別給出了平行光、天球光、水下和水上點(diǎn)光源形成的水下光場(chǎng)的計(jì)算方法,并考慮到水面下表面反射的影響。假定光線的傳輸及散射現(xiàn)象具有獨(dú)立性和疊加性,則在以上各類(lèi)光源的任意組合下求解水下光場(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)。
水下光電成像技術(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)行分析,而將本文的方法作為必要的參考。