水下光場的迭代求解

何大華，李陽陽，周少杰

（華中光電技術研究所—武漢光電國家研究中心, 湖北 武漢 430223）

1 引 言

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

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

本文以水體吸收和散射的基本過程為基礎，在穩定光源照射以及體散射函數為球形對稱的情況下，定量求解光線經水體多次散射后形成的水下光場分布。求解水下光場的過程闡明了水體對光線吸收和散射的微觀機制，精確的水下光場數據為嚴格推導水體PSF和MTF提供了關鍵信息。

2 水下光場計算模型

2.1 基本概念

一般認為，準直窄光束在水下傳輸時能量衰減遵循朗伯-比爾定律[13]，即

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

對于水下點光源，設光強為I0，則在水下傳輸距離r后未被吸收和散射的能量形成的徑向照度為

多次散射：光束或光子在水下傳輸過程中被水分子或水中雜質散射從而多次改變傳輸方向的現象。

球面照度：空間中某點的球面照度為全空間范圍內入射到該點處一個小球內的光通量總和與該小球的表面積之比，可表示為4π立體角范圍內對亮度的積分。球面照度也稱為標量照度。

水下初始光場：在指定的水域中，光源發射的光能量進入水體傳輸過程中未被水體吸收和散射的部分在水下形成的球面照度分布。

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

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

2.2 水下光場的形成過程

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

Ⅰ.點光源發出的光線經擴散和水體衰減后形成水下初始光場Et0。

Ⅱ.水下初始光場Et0使整個水域產生散射，形成一次水下散射光場Es1，疊加在Et0上形成一次水下光場Et1=Et0+Es1。

Ⅲ.一次水下光場Et1使整個水域產生散射，形成二次水下散射光場Es2，疊加在Et0上形成二次水下光場Et2=Et0+Es2。

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

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

Ⅵ.Eti將收斂于水下光場Et。

根據以上思路，下面給出一般意義上的水下光場求解方程式。

如圖1，點光源S發出的光在水下傳輸過程中會進行球面擴散，還會被水體吸收、散射，一部分光線繼續沿原來的路徑傳播，在P點處形成的水下初始光場為

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

圖1 水下光場示意圖Fig.1 Diagram of the underwater light field

設水下光場為Et(x,y,z)，則有

P點處的水下光場等于水下初始光場與水下散射光場之和，因而有

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

3 典型水下光場的計算方法

3.1 太陽照射下水下光場的計算方法

2.2節給出了理想情況下水下光場的計算方法，實際中總會存在一些邊界條件，下面給出太陽照射下足夠大平靜水面下的水下光場。

僅考慮太陽照度，忽略天空背景亮度的影響，太陽光以一定角度平行入射進入水中，并發生反射和折射。折射光線經吸收和散射后，剩下部分繼續向前傳播，隨深度增加能量越來越小。散射光線同樣會被吸收、散射或者繼續傳播，在水下形成多次散射，進而形成穩定的水下光場。

如圖2所示，從對稱性考慮，太陽照射下的水下光場的數值僅僅與深度有關，而與水平坐標無關，因此，水下光場函數可用深度h為自變量的一元函數Et(h)表示。

圖2 太陽照射下的水下光場Fig.2 Underwater light field from the sun

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

設與入射太陽光垂直方向上的照度為Esun，則水下深度h處的水下初始光場為

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

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

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

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

還應考慮圓C上的微體積元散射光經水面下表面反射后對P點形成的球面照度。M點發出的光線經A點反射可到達P，可等效為從M'發出的光線直接到達P。同理，N點發出的光線經B點反射也可到達P，可等效為從N'發出的光線直接到達P，故這組反射光線可等效從圓D發出，并經過了長度為PM'的水程到達P點，記為。

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

綜合以上結果，圓C上的微體積元在P點形成的球面照度dEs為

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

由此可知P點的水下散射光場Es(h)為

故水下光場Et(h)為

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

其中令深度h及積分變量r和z的細分間隔均為Δ，并對式（12）進行數值化得

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

由于水體對光線的衰減嚴重，水體散射對距離超過10AL的區域影響可以忽略，因此在計算時，可適當限定計算范圍，基本不影響計算精度，例如計算10AL深度內的水下光場時，可令k=1,2,3,···,N，并滿足NΔ＞20AL。據此，將式（14）改為以下迭代形式

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

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

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

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

Δ取值的依據是，當Δ=0.02 AL時，水下光場精度足夠，此例中AL=10m，因而Δ=0.2m。

圖5是水深100 m（10AL ）以內仿真程序輸出的水下光場迭代曲線。

圖5 太陽照射下水下光場迭代曲線Fig.5 Iteration curves of the underwater light field from the sun

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

3.2 均勻亮度天空背景下的水下初始光場

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

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

如圖6所示，天頂角為α的微球帶對應的立體角為dΩ=2πsinαdα。

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

設均勻天空的亮度為Lsky，對水下初始光場求解時，可將微球帶發出的光線等效為入射角為α的一組平行光，而且垂直于入射方向上的照度為dE=2Lskyπsinαdα，則可得均勻亮度天空背景下水下初始光場為

編寫Matlab仿真程序計算Et(h)，假設天空亮度為Lsky=5000cd/m2，水體光學參數及Et(h)的深度精度同3.1節，則Et(h)在深度10AL之內的迭代曲線如圖7所示。其中各曲線的意義與圖5相同。

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

3.3 點光源形成的水下初始光場

3.3.1 點光源在水面以下

設水深h0處有均勻發光點光源S，光強為I0，水體衰減系數、吸收系數及散射系數同上，計算水面平靜時光源S形成的水下初始光場。

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

圖8 水下點光源形成的水下初始光場Fig.8 Underwater initial light field by an underwater point light source

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

S點發出的光線經水面下表面反射的光線可以等效為S點關于水面的對稱點S′點發出，經過了水程S′P的衰減以及反射率Rα的修正，則

3.3.2 點光源在水面以上

設水面以上高h0處有均勻發光點光源S，光強為I0，水體衰減系數、吸收系數及散射系數同上，下面計算水面平靜時光源S形成的水下初始光場。

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

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

圖9 水上點光源形成的水下初始光場Fig.9 Underwater initial light field by an overwater point light source

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

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

比較式（20）和式（21）并化簡得

3.4 小結

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

在自帶光源的主動水下光電成像應用中，為了提高功率密度和減小后向散射，光源的照明區域往往局限在有限的角度范圍內，這時水下光場也可以通過上述方法求解，因此該計算模型可擴展應用到邊界條件更為復雜一些的情形。

已知水下光場分布，則易求得水體的表觀光學參數[13]，因此，該計算模型還可直接用于求解水體的表觀光學參數（AOP）。

4 結 論

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

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

然而，在實際問題中，光源特性并不穩定，水體散射特性在空間上也不均勻，且存在湍流時變擾動[35-37]，故水下光場是動態變化的，且水下各種噪聲也對計算結果有較大影響。因此，在具體的水下光場求解過程中，需要結合估計、現場測量、大數據處理等經驗方法進行分析，而將本文的方法作為必要的參考。