陳 希, 于 斌, , 何君姮, 孔凡龍, 李宗葉

(1.解放軍理工大學(xué) 氣象海洋學(xué)院, 江蘇 南京 211101; 2.91206部隊訓(xùn)練部, 山東 青島 266108; 3.東海艦隊37分隊, 浙江 寧波 315122; 4.71901部隊, 山東 聊城 252000)

海水透明度是指透明度盤鉛直沉入海水中的最大可見深度, 是軍事海洋環(huán)境的關(guān)鍵要素, 對探潛、反潛及潛艇隱蔽等活動具有重要的軍事意義。目前,獲取海水透明度包括, 海水透明度盤測量、經(jīng)驗公式計算、激光雷達探測、衛(wèi)星資料反演等多種方法。

利用透明度盤測量海水透明度是海洋調(diào)查工程中廣泛使用的一種常規(guī)方法, 但是受海況、人為差異等外界因素的影響較大, 觀測數(shù)據(jù)具有一定的主觀性和隨意性, 存在較大誤差; 海水透明度計算的另一種方法是基于經(jīng)驗公式, 如張緒琴等[1]基于海水均勻的假設(shè), 提出了海水透明度與光束衰減系數(shù)和漫射衰減系數(shù)之和成反比的經(jīng)驗公式, 適用于均勻海水透明度計算; 在激光測量透明度方面, 姜璐等[2]分析了機載激光雷達最大探測深度與海水透明度之間的比例關(guān)系, 該方法適用于水深 30 m以淺海域;在衛(wèi)星遙感反演透明度方面, 何賢強等[3]使用SeaWiFS衛(wèi)星遙感資料, 建立半經(jīng)驗?zāi)Ｊ椒囱莺Ｋ该鞫? 該方法適用于一類水體。

上述研究中, 透明度盤方法存在主觀性、隨意性缺點; 利用表觀和固有光學(xué)性質(zhì)數(shù)據(jù)計算透明度的方法一般是基于海水光學(xué)均勻性假設(shè)開展的相關(guān)研究。本文基于垂直非均勻海水分層思想, 引入蒙特卡羅方法, 提出并建立了一種基于海水固有光學(xué)性質(zhì)的海水透明度計算模型。

1 海水透明度計算模型

利用透明度盤測量海水透明度原理示意圖如圖1。當(dāng)入射光穿透海水照射到透明度盤上, 透明度盤反射太陽光, 產(chǎn)生透明度盤上涌光, 盤具有一定的亮度。同時海水由于對光有后向散射等作用, 產(chǎn)生海水上涌光, 海水具有一個背景亮度。于是在透明度盤處,透明度盤和海水存在亮度差異, 即固有對比度。透明度盤上涌光和海水上涌光傳遞到人眼處, 兩者仍然存在亮度差異, 即為相對對比度, 若此時人眼能感知到相對對比度, 則透明度盤可見。不斷下放透明度盤的深度, 當(dāng)人眼感知不到相對對比度時, 認(rèn)為透明度盤已經(jīng)下放到透明度深度, 則可測得海水透明度。本文在仿真晴空條件下透明度盤測量原理基礎(chǔ)上, 建立非均勻海水透明度計算模型的基本思路如下: 假設(shè)透明度盤位于深度D處, 其上涌光相對量為uD, 通過uD計算透明度盤的固有對比度, 代入非均勻水體對比度傳輸方程求解到達人眼處的相對對比度, 與人眼閾值進行比較, 若相對對比度大于人眼閾值, 則D達到了海水透明度所在深度,得到海水透明度值, 反之, 增大D, 進行新一輪計算, 計算流程如圖2所示。其中, 海水上涌光相對量、透明度盤固有對比度及相對對比度的計算方法是該模型的關(guān)鍵, 下文將對這3個問題進行詳細(xì)論述。

圖1 透明度盤測量海水透明度原理圖Fig.1 Schematic of transparency disc measuring water transparency

圖2 基于固有光學(xué)量計算海水透明度算法流程圖Fig.2 Calculating water transparency based on inherent optical quantization

1.1 海水上涌光相對量的計算方法

假設(shè)透明度盤在D處的上涌光相對量為uD, 向上輻照度為Eu, 向下輻照度為Ed, 則有:

Eu和Ed分別描述了水下光場上涌光和下涌光能量的分布情況, 而水下光場的能量分布, 可以利用蒙特卡羅方法, 跟蹤光子運動過程, 達到穩(wěn)定狀態(tài)時,對光子能量進行統(tǒng)計和累加得到,Eu和Ed的計算步驟如下(圖3)。

圖3 蒙特卡羅方法跟蹤光子水下運動過程流程圖Fig.3 Flow diagram of tracking photon underwater movement based on Monte Carlo method

第一步: 初始化光子能量、入射方向和位置;

第二步: 計算光子下一個散射點的位置及能量;

第三步: 根據(jù)散射前后位置關(guān)系, 若光子自下而上穿過深度D的海水薄層, 則將光子能量累加到向上輻照度, 若向下穿過, 則累加到向下輻照度, 若不能穿越, 光子繼續(xù)運動;

第四步: 根據(jù)能量閾值判斷光子能量是否耗盡,若耗盡則認(rèn)定光子死亡, 轉(zhuǎn)向第一步, 進行新一輪模擬;

第五步: 若光子并未死亡, 則判斷光子是否到達海表或海底, 若到達海表則判斷是否發(fā)生全反射,若發(fā)生全反射, 計算反射后光子的能量、坐標(biāo)和運動方向余弦, 若出水, 則進行新一輪模擬; 若到達海底,則計算反射或者散射后的方向和能量; 若未到達海表和海底, 則計算光子散射后的方向余弦, 轉(zhuǎn)向第二步。

在上涌光相對量計算過程中, 涉及光子散射點位置、運動方向余弦和能量的計算, 光子運動到海底、海面時運動方向余弦和能量的計算, 在透明度盤深度處的海水上涌光和下涌光輻照度的統(tǒng)計計算,相關(guān)計算方案如下。

1.1.1 光子散射位置、運動方向和能量的計算

1.1.1.1 光子散射位置

以水平的海表面為x-y軸, 垂直水平面向下為z軸建立直角坐標(biāo)系, 如圖4所示。假設(shè)光子在p= (x,y,z)處發(fā)生散射后, 其方向向量為μ={μx,μy,μz}, 能量為E, 下一次散射點的坐標(biāo)p′ = (x′,y′,z′), 方向向量為能量為E′, 由于光子運動的幾何路徑遠小于水深, 可認(rèn)為光子在兩次散射之間所經(jīng)歷的幾何路徑為直線, 則下一個散射點的坐標(biāo)為

圖4 光子水下運動參數(shù)幾何關(guān)系示意圖Fig.4 Schematic of relation among geometric parameters of photon underwater movement

式中l(wèi)表示光子走過的幾何路徑, 在均勻水體中,l可由光學(xué)衰減長度s和衰減系數(shù)c確定[6]:

在非均勻水體中, 可將水體分為若干層, 每一層中的水體光學(xué)性質(zhì)近似均勻, 假設(shè)每一層光子經(jīng)過的光學(xué)衰減長度為si, 衰減系數(shù)為ci, 則光子穿過的幾何路徑長度為:

根據(jù)蒙特卡羅方法的基本原理[7], 光學(xué)衰減長度s= -lnq, 式中q為區(qū)間[0,1]上服從均勻分布的隨機數(shù)。

1.1.1.2 光子散射后運動方向

由圖3中的散射后各參量的幾何關(guān)系, 通過余弦定理可計算光子散射后運動方向向量μ′:

式中,φ為水平散射角, 滿足φ=2πq,q為區(qū)間[0,1]上服從均勻分布的隨機數(shù)。sμ表示散射角sθ余弦, 根據(jù)Robert 等[9]的研究有

其中q為區(qū)間[0,1]上服從均勻分布的隨機數(shù),g是不對稱因子。

1.1.1.3 光子散射后能量

光子散射后的能量:

式中,sp分別為光子散射位置水體的散射系數(shù)。

1.1.2 海底對光子運動方向和能量的影響

假設(shè)光子碰撞海底前的坐標(biāo)為p= (x,y,z), 方向余弦為μ={μx,μy,μz}, 能量為E, 碰撞后的坐標(biāo)為p′ = (x′,y′,z′), 方向余弦為μ′= {μx′,μy′,μz′}, 能量為E′, 海底深度為H, 海底對光的反射可以看作是郎伯反射面, 則根據(jù)朗伯反射面的性質(zhì), 可以得到如下關(guān)系[10]:

其中,q和x分別為區(qū)間[0,1]上服從均勻分布的隨機數(shù)。若海底反射率為bρ, 則光子撞擊海底后, 其能量變化為

由方程組(7)可以計算光子碰撞海底后的坐標(biāo)以及運動的方向余弦, 根據(jù)(8)式, 可以更新光子的能量。

1.1.3 海面對光子運動方向和能量的影響

光子運動到海面后的方向和能量會受到波浪起伏的影響, 波浪起伏可以看作是若干小的傾斜平面組成的, 根據(jù)斯涅爾定律, 利用傾斜平面的斜率和入射方向可計算光束在海表處反射或折射的運動方向。由于出水后蒙特卡羅方法對單個光子模擬過程結(jié)束, 故在此只討論發(fā)生全反射時光子反射方向和能量變化情況。

設(shè)光子在坐標(biāo)p=(x0,y0,z0)處發(fā)生全反射,p所在平面的方程可由傾斜平面的仰角和方位角進行計算, 全反射前光子的運動軌跡直線方程可由全反射點、運動方向余弦唯一確定, 則可以確定直線上z坐標(biāo)高于p的任意一點X, 根據(jù)直線與平面位置關(guān)系可以計算X關(guān)于傾斜平面的對稱點X′, 向量up=X-p的單位向量即為反射方向向量, 計算得到:

θn、φn分別為波浪面的仰角和方位角。φn可認(rèn)為是[0,2π]上均勻分布,θn的概率分布函數(shù)[8]為:

其中,σ2= 0.003 + 0.00512v,v為海表風(fēng)速。根據(jù)概率分布函數(shù)的歸一性, 可以計算θn:

(11)式確定了傾斜平面的法線方向。

光束從海水中射向空氣中時, 若光線相對于波浪面法線的入射角iθ滿足光的全反射條件:

則在波浪面上, 光線發(fā)生全反射, 此時, 光子的運動方向與原方向關(guān)于波浪面法線對稱, 能量等于原能量與海表反射率的乘積。

反之, 若

則光子穿透水面折射入大氣, 光子運動符合折射定律, 通過折射定律可以計算光子的出射方向, 其能量等于原能量與海水折射率的乘積, 即則光子能量的變化為:

其中E為穿越前能量, 穿越后為E′,naw為海水相對于空氣的折射率, 一般情況下, 空氣的折射率為 1,而海水相對空氣的折射率可取為naw=1.34。

1.1.4 上涌光輻照度和下涌光輻照度的計算

深度D處的上涌光輻照度Eu和下涌光輻照度Ed即可通過對光子能量統(tǒng)計和累加的方式計算出來。若光子向上穿越深度D處的均勻薄水層, 即μz＜0, 則將光子在垂直方向上的能量累加到上涌光輻照度中, 即

若光子向下穿越深度D處的均勻薄水層, 即μz＞0,則將光子在垂直方向上的能量累加到下涌光輻照度中, 即

根據(jù)(1)式, 即可計算深度D處的上涌光相對量。

1.2 固有對比度和相對對比度的計算

使用上涌光相對量計算透明度盤在深度D處的固有對比度CD[1]:

式中Dρ為透明度盤漫反射率。

在垂直方向上, 因為人眼距離水面較近, 故人眼到水面之間的衰減可以忽略不計, 將固有對比度作為輸入, 通過對比度傳輸方程式(16), 即可計算透明度盤到達眼睛的相對對比度C為[9]:

式中l(wèi)i為透明度盤向上直到水面距離上各光學(xué)均勻?qū)拥暮穸?Ki為漫射系數(shù), 可在蒙特卡羅方法模擬光子水下運動時根據(jù)其定義進行計算[14]。

2 數(shù)值實驗

利用上述建立的透明度計算模型, 進行透明度計算實驗, 采用 Australian Antarctic Division的BROKE-WEST_ACS數(shù)據(jù)集中 2006年 1、2月份59.8°～68.1°S、30°～80°E 海區(qū)的實測固有光學(xué)性質(zhì)垂直剖面數(shù)據(jù)、海表面風(fēng)速和透明度盤觀測數(shù)據(jù)進行驗證, 具體步驟為:

第一步: 設(shè)置透明度盤的初始深度d為 0.1 m,固有光學(xué)性質(zhì)數(shù)據(jù)的垂直分辨率為2 m, 實驗海區(qū)水深取100 m(據(jù)海水透明度觀測數(shù)據(jù)統(tǒng)計表明, 一般不超過70 m), 海水分層為50層。

第二步: 運用蒙特卡羅方法模擬水下光場, 其光子總數(shù)設(shè)為 1 000萬個, 根據(jù) Duntely等[4]和Gilbert[5]的研究, 設(shè)光子的初始能量值為 1個單位,光子死亡能量閾值取為0.000 1, 初始位置隨機設(shè)定,初始方向為太陽天頂角。對模擬過程中穿過了d的薄水層的光子, 采用式(1)、式(13)、式(14)對計算d上的上涌光相對量uD, 并通過漫射衰減系數(shù)的定義,計算d上的漫射衰減系數(shù)k;

第三步: 將uD代入式(15), 計算透明度盤在d上的固有對比度Cd, 并通過式(16)計算到達人眼的相對對比度Cri;

第四步: 比較相對對比度Cri與人眼對比度閾值0.002, 若Cri＞ 0.002, 則將d增加0.1 m并轉(zhuǎn)向步驟1), 否則, 認(rèn)為d為海水透明度D并輸出結(jié)果。計算流程如圖5所示。

圖5 數(shù)值計算流程圖Fig.5 Circuit of numerical methods

將模型計算的海水透明度值與37個站點的透明度盤實測數(shù)據(jù)進行了對比, 如圖6所示, 計算值與實測值接近, 平均絕對誤差為1.1 m, 平均相對誤差為9%, 取得了較好的效果, 說明本文的模型是可行的。

圖6 計算透明度盤模型計算結(jié)果與實測值對比圖Fig.6 Comparison of transparency disk model calculation results and measurements

3 總結(jié)和討論

1) 本文建立了基于非均勻海水固有光學(xué)性質(zhì),仿真晴空條件下透明度盤測量海水透明度的基本過程, 將蒙特卡羅方法與對比度傳輸方程相結(jié)合, 計算海水透明度的數(shù)學(xué)模型。

2) 該模型克服了透明度盤測量方式易受天氣、海況和人工影響等缺點, 計算模型僅為依賴于海水固有光學(xué)性質(zhì), 對不同海區(qū)具有通用性, 為夜間探測海水透明度增加了一種新途徑。

3) 數(shù)值實驗結(jié)果與海水透明度觀測數(shù)據(jù)比較表明, 平均絕對誤差1.1 m, 證明該方法可以客觀、定量、準(zhǔn)確的計算海水透明度。

4) 受到海水固有光學(xué)性質(zhì)實測數(shù)據(jù)的限制, 數(shù)值實驗驗證樣本相對較少, 下一步研究中, 將利用不同海區(qū)、不同季節(jié)的海水透明盤觀測數(shù)據(jù)對模型計算精度進行進一步驗證和改進。

[1] 張緒琴.海水透明度[J].海洋湖沼通報, 1982, 4:14-17.

[2] 姜璐, 朱海, 李松, 等.機載激光雷達最大探測深度同海水透明度的關(guān)系[J].激光與紅外, 2005, 35(6):397-399.

[3] 何賢強.利用 SeaWiFS反演海水透明度的模式研究[J].海洋學(xué)報, 2004, 26(5): 55-63.

[4] Mobley C D.Light and Water [M].New York: Academic Press, 1994.

[5] Capone A, Digaetano T, Grimaldi A, et al.Measurements of light transmission in deep sea with the AC9 trasmissometer[J].Nuclear Instruments & Methods in Physics Research A, 2002, 487(3): 423-434.

[6] Robert A.Leathers, Monte Carlo Radiative Transfer Simulations for Ocean Optics[M].Washington: NRL,2004: 9.

[7] Lerner R M, Summers J D.Monte Carlo description of time and space resolved multiple forward scatter in natural water [J].Applied Optics, 1982, 21(5):861-869.

[8] Xu Z.A DNS capability for obtaining underwater light field and retrieving upper ocean conditions via in-water light measurements [D].Massachusetts: Massachusetts Institute of Technology Department of Mechanical Engineering, 2011.

[9] Jerlov N G.Marine Optics [M].Amsterdam: Elsevier Scientific Public, 1976.

[10] 陳文革, 黃鐵俠.機載海洋激光雷達的試驗研究[J].電子學(xué)報, 1998, 26(9): 21-24.