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

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

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

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

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

1 海水透明度計(jì)算模型

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

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

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

1.1 海水上涌光相對(duì)量的計(jì)算方法

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

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

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

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

第二步: 計(jì)算光子下一個(gè)散射點(diǎn)的位置及能量;

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

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

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

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

1.1.1 光子散射位置、運(yùn)動(dòng)方向和能量的計(jì)算

1.1.1.1 光子散射位置

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

圖4 光子水下運(yùn)動(dòng)參數(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]上服從均勻分布的隨機(jī)數(shù)。

1.1.1.2 光子散射后運(yùn)動(dòng)方向

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

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

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

1.1.1.3 光子散射后能量

光子散射后的能量:

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

1.1.2 海底對(duì)光子運(yùn)動(dòng)方向和能量的影響

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

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

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

1.1.3 海面對(duì)光子運(yùn)動(dòng)方向和能量的影響

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

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

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

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

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

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

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

反之, 若

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

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

1.1.4 上涌光輻照度和下涌光輻照度的計(jì)算

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

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

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

1.2 固有對(duì)比度和相對(duì)對(duì)比度的計(jì)算

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

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

在垂直方向上, 因?yàn)槿搜劬嚯x水面較近, 故人眼到水面之間的衰減可以忽略不計(jì), 將固有對(duì)比度作為輸入, 通過對(duì)比度傳輸方程式(16), 即可計(jì)算透明度盤到達(dá)眼睛的相對(duì)對(duì)比度C為[9]:

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

2 數(shù)值實(shí)驗(yàn)

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

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

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

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

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

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

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

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

3 總結(jié)和討論

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

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

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

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

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

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

[3] 何賢強(qiáng).利用 SeaWiFS反演海水透明度的模式研究[J].海洋學(xué)報(bào), 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ī)載海洋激光雷達(dá)的試驗(yàn)研究[J].電子學(xué)報(bào), 1998, 26(9): 21-24.