基于機(jī)載LiDAR測(cè)深水體波形的漫衰減系數(shù)提取方法

亓超，周豐年，吳敬文，宿殿鵬，王賢昆，陽(yáng)凡林*

( 1.山東科技大學(xué) 測(cè)繪與空間信息學(xué)院，山東 青島 266590；2.自然資源部海洋測(cè)繪重點(diǎn)實(shí)驗(yàn)室，山東 青島 266590；3.長(zhǎng)江水利委員會(huì)水文局長(zhǎng)江口水文水資源勘測(cè)局，上海 200136)

1 引言

漫衰減系數(shù)(Diffuse Attenuation Coefficient,Kd)是水體的固有光學(xué)屬性，反映了激光在水體中受散射、吸收所引起的能量衰減變化，是分析海水水質(zhì)的重要參數(shù)之一[1]。漫衰減系數(shù)的大小主要取決于水體中的浮游生物含量（葉綠素濃度）、懸浮泥沙含量（渾濁度）、營(yíng)養(yǎng)鹽含量（黃色物質(zhì)、溶解有機(jī)物質(zhì)、鹽度指標(biāo)）以及其他污染物等因素，研究漫衰減系數(shù)可對(duì)水體富營(yíng)養(yǎng)化、水體渾濁度驟增以及水體污染等環(huán)境變化實(shí)現(xiàn)實(shí)時(shí)、有效監(jiān)測(cè)，在海洋科學(xué)研究、海洋工程應(yīng)用、海洋漁業(yè)研究和海岸帶環(huán)境監(jiān)測(cè)等方面發(fā)揮重要的作用[2–3]。

由于水體的流動(dòng)性，導(dǎo)致了水中物質(zhì)的非均勻分布，增加了水體中漫衰減系數(shù)的獲取難度。目前，國(guó)內(nèi)外針對(duì)漫衰減系數(shù)的提取主要通過(guò)實(shí)地測(cè)量、被動(dòng)遙感以及主動(dòng)遙感3種方式進(jìn)行。實(shí)地測(cè)量利用賽克盤(pán)和透射計(jì)等工具測(cè)量水體透明度，進(jìn)而反演得到漫衰減系數(shù)，這些方式需以測(cè)量船為載體，效率較低且覆蓋范圍小，難以快速獲取大面積水域的漫衰減系數(shù)空間分布[2,4–5]；通過(guò)衛(wèi)星遙感影像數(shù)據(jù)可間接反演漫衰減系數(shù)，這種被動(dòng)遙感方式可對(duì)大范圍海域的漫衰減系數(shù)進(jìn)行求取，但是無(wú)法在云霧天氣、夜晚實(shí)施測(cè)量，且由于衛(wèi)星設(shè)備本身的限制，測(cè)量精度難以保證，成果分辨率較低[6–8]；全波形機(jī)載LiDAR測(cè)深（Airborne LiDAR Bathymetry, ALB）技術(shù)的出現(xiàn)彌補(bǔ)了上述方法的缺陷，運(yùn)用水體散射回波的強(qiáng)度變化可實(shí)現(xiàn)漫衰減系數(shù)的快速反演[4,9–13]。Billard等[9]利用ALB技術(shù)在不同深度下的輻照度進(jìn)行指數(shù)擬合，進(jìn)而反演得到漫衰減系數(shù)的方法，該方法需要水底底部回波強(qiáng)度和相應(yīng)深度信息；Ding等[10]提出了一種基于單波長(zhǎng)ALB全波形數(shù)據(jù)的漫衰減系數(shù)計(jì)算方法，克服了傳統(tǒng)方法的局限性，但未對(duì)水體散射回波進(jìn)行擬合分解，易受噪聲影響。

因此，針對(duì)目前船載實(shí)地測(cè)量效率與分辨率低、衛(wèi)星遙感反演精度與分辨率較低的局限性，本文結(jié)合ALB技術(shù)高精度、高分辨率、靈活機(jī)動(dòng)、快速高效的特點(diǎn)，提出了一種基于機(jī)載LiDAR測(cè)深水體波形的漫衰減系數(shù)提取方法。ALB原始波形數(shù)據(jù)經(jīng)小波去噪之后，利用基于分層異構(gòu)模型的波形分解算法[14–15]（水面?高斯函數(shù)、水體?雙指數(shù)函數(shù)及水底?高斯函數(shù)）分解得到水體散射回波信號(hào)，進(jìn)而獲取大面積水域漫衰減系數(shù)的空間分布。

2 基于機(jī)載LiDAR測(cè)深水體波形的漫衰減系數(shù)提取方法

機(jī)載LiDAR測(cè)深水體波形的漫衰減系數(shù)提取方法充分利用了ALB技術(shù)高精度、高分辨率、靈活機(jī)動(dòng)、快速高效的獨(dú)特優(yōu)勢(shì)，考慮到激光在水體中的衰減特性以及計(jì)算效率，通過(guò)分層異構(gòu)模型的機(jī)載LiDAR波形分解算法分解得到水體散射回波信號(hào)，并用其與水體散射回波強(qiáng)度方程構(gòu)建相應(yīng)模型，用以計(jì)算每個(gè)激光測(cè)深點(diǎn)的漫衰減系數(shù)。該算法的步驟為：首先獲取機(jī)載LiDAR測(cè)深原始波形數(shù)據(jù)，對(duì)其進(jìn)行去噪；然后利用分層異構(gòu)模型對(duì)去噪后的波形數(shù)據(jù)進(jìn)行分解，分解得到水體散射回波信號(hào)；再通過(guò)水體散射回波信號(hào)分別計(jì)算得到雙段指數(shù)所對(duì)應(yīng)的初始漫衰減系數(shù)；最后根據(jù)雙段指數(shù)的采樣時(shí)間確定相應(yīng)權(quán)重，對(duì)初始漫衰減系數(shù)取加權(quán)平均，得到該激光測(cè)深點(diǎn)最終的漫衰減系數(shù)值。本文算法具體流程如圖1所示。

2.1 基于分層異構(gòu)模型的波形分解算法

機(jī)載LiDAR測(cè)深系統(tǒng)利用具有較強(qiáng)透水能力的藍(lán)綠色波段(532 nm)激光進(jìn)行探測(cè)，其激光接收器能夠以數(shù)字化的形式記錄每個(gè)激光脈沖的全部回波波形信號(hào)（圖2）[15–16]，包括水面反射回波、水體散射回波、水底反射回波以及噪聲4部分[17]，可描述為

式中，t表示時(shí)間刻度；fs(t)表示水面反射回波；fc(t)表示水體散射回波；fb(t)表示水底反射回波；fN(t)表示原始波形的噪聲；f(t)為t時(shí)刻回波的振幅。

圖1 本文方法流程圖Fig.1 Flow diagram of the new method in this study

圖2 機(jī)載LiDAR測(cè)深波形示意圖Fig.2 Airborne LiDAR bathymetric waveform

其中，噪聲包括背景噪聲和傳感器內(nèi)部噪聲，其存在會(huì)導(dǎo)致波形分解函數(shù)初始參數(shù)估計(jì)不準(zhǔn)等問(wèn)題，容易獲得含有粗差的波形，影響波形分解精度。因此，在波形分解之前，需對(duì)原始回波信號(hào)中的噪聲進(jìn)行相應(yīng)的濾除，文中使用小波算法對(duì)噪聲進(jìn)行處理[15,18]。

波形分解是機(jī)載LiDAR測(cè)深數(shù)據(jù)處理的重要環(huán)節(jié)，是漫衰減系數(shù)提取的基礎(chǔ)。根據(jù)機(jī)載LiDAR測(cè)深回波信號(hào)各部分的波形特征，選擇合適的函數(shù)構(gòu)建分層異構(gòu)模型（水面?高斯函數(shù)、水體?雙指數(shù)函數(shù)及水底?高斯函數(shù)）進(jìn)行分解，從而將疊加在水面反射回波和水底反射回波中的水體散射回波精確地分解出來(lái)，得到水體散射回波所對(duì)應(yīng)的函數(shù)。

考慮到ALB系統(tǒng)發(fā)射和接收的激光脈沖信號(hào)近似服從高斯分布[19–20]，因此，根據(jù)水面和水底反射回波相應(yīng)特性分析，文中采用高斯函數(shù)（公式(2)）分解得到水面、水底反射回波信號(hào)，

式中，i表示水面（水底）反射回波信號(hào)的分解結(jié)果；Ai、μi和σi分別為水面（水底）高斯函數(shù)的波峰值、波峰位置和半幅波寬。

從去噪后的ALB波形數(shù)據(jù)中剔除水面、水底反射回波信號(hào)，即可得到初始水體散射回波信號(hào)。根據(jù)光在水體中輻射傳輸定理可知，激光在水體中傳播時(shí)，其脈沖能量隨水深的增加而呈現(xiàn)指數(shù)衰減。針對(duì)水體散射回波的分解，Ding等[21]在三角形函數(shù)[22]和四邊形函數(shù)[23]分解方法的基礎(chǔ)之上，提出了一種改進(jìn)的四邊形函數(shù)，將水體散射回波下降沿的一段采用指數(shù)函數(shù)描述。由于水體渾濁度在垂直剖面上非均勻分布導(dǎo)致激光脈沖的衰減程度不同[24]，綜合激光在水體中衰減特性和分解效率的問(wèn)題，如圖3所示，文中將水體沿垂直剖面分為兩層，利用雙指數(shù)函數(shù)ABCD（公式(3)）分解得到水體散射回波信號(hào)

式中，ax、bx、cx、dx分別表示雙指數(shù)函數(shù) 4 個(gè)頂點(diǎn) (A、B、C、D)的位置時(shí)刻；by、cy、dy為 B、C、D的振幅值。

非線性擬合需要預(yù)設(shè)初值，以降低迭代次數(shù)，增加擬合速度。文中遵循非線性最小二乘(Non-linear Least Squares, NLS)方法（公式 (4)），采用 Levenberg-Marquardt算法對(duì)13個(gè)初始參數(shù)（水面–高斯函數(shù)3個(gè)、水底?高斯函數(shù)3個(gè)、水體?雙指數(shù)函數(shù)7個(gè)）進(jìn)行迭代優(yōu)化，直到擬合偏差達(dá)到最小為止，以求得各參數(shù)的精確值。將優(yōu)化后的水體反向散射回波參數(shù)代入公式(3)，得到最優(yōu)的水體散射回波分解結(jié)果，進(jìn)而用于漫衰減系數(shù)的提取。

圖3 雙指數(shù)函數(shù)分解水體散射回波示意圖Fig.3 Decomposition of water column contribution by doubleexponential function

式中，λ表示所求參數(shù)；S(λ)表示擬合偏差的平方和；Fj表示波形分解結(jié)果在j時(shí)刻的振幅。

2.2 漫衰減系數(shù)Kd的提取

機(jī)載LiDAR測(cè)深系統(tǒng)在深度為Di處的水體散射回波強(qiáng)度方程可由下式來(lái)描述[4,17]

式中，Di表示水層深度；Pc(Di)表示深度Di時(shí)的水體散射回波強(qiáng)度；P表示激光發(fā)射強(qiáng)度；表示大氣雙

e程損失因子；AR表示激光接收器的視場(chǎng)面積；ηe和ηR分別表示激光發(fā)射器的光學(xué)效率和激光接收器的光學(xué)效率；F表示接收視場(chǎng)角損失因子；LS是通過(guò)表面?zhèn)鬏數(shù)膿p失，即表面反照率；β（?）表示體散射函數(shù)；θ表示激光傳入水體時(shí)的入射角，即天底角；θw表示射入水體后的折射角；nw表示水體折射率；H表示飛機(jī)航高；Kd表示532 nm波段處的漫衰減系數(shù)；c表示激光在真空的傳播速度；ti表示水層深度Di處的時(shí)刻；ts表示水面反射回波的峰值位置。

對(duì)該方程分析可知，激光脈沖回波強(qiáng)度隨水深深度和漫衰減系數(shù)增加呈指數(shù)衰減，因此，可以通過(guò)水體散射回波所對(duì)應(yīng)的雙指數(shù)函數(shù)反演水體的漫衰減系數(shù)。文中假設(shè)外業(yè)測(cè)量過(guò)程中飛機(jī)航高保持穩(wěn)定，激光發(fā)射器的天底角以及接收視場(chǎng)角保持恒定，且所有損耗因素均得到較好的控制，對(duì)水體散射回波強(qiáng)度方程進(jìn)行簡(jiǎn)化，并將公式(6)代入公式(5)，并簡(jiǎn)化公式(5)，得到簡(jiǎn)化后的水體散射回波強(qiáng)度方程為

由于近岸水深通常小于30 m，與飛機(jī)航高H（通常為400～500 m）相比，參數(shù)Di對(duì)W的影響可以忽略不計(jì)，因此，W可認(rèn)為是一個(gè)常數(shù)。但是W涉及未知參數(shù)太多，不能直接確定。通過(guò)分析發(fā)現(xiàn)，公式(3)和公式(7)均可描述為如下指數(shù)型函數(shù)的形式，

式中，y(t)為指數(shù)型函數(shù)；m、n分別為指數(shù)型函數(shù)的系數(shù)。

考慮到傳統(tǒng)漫衰減系數(shù)提取算法所存在的局限性[9–10]，結(jié)合2.1節(jié)中水體散射回波雙指數(shù)函數(shù)（公式(3)）與簡(jiǎn)化后的水體散射回波強(qiáng)度方程(公式(7))構(gòu)建模型，從而計(jì)算漫衰減系數(shù)Kd值。以雙指數(shù)函數(shù)（圖3）BC段為例，當(dāng)bx≤t

對(duì)公式(9)做方程變換，得到BC段指數(shù)所對(duì)應(yīng)的初始漫衰減系數(shù)Kd1，同理，計(jì)算得到CD段指數(shù)所對(duì)應(yīng)的初始漫衰減系數(shù)Kd2，如下式所示：

為了保證所得漫衰減系數(shù)的準(zhǔn)確性，以BC段、CD 段的采樣時(shí)間（Δt1=cx?bx與 Δt2=dx?cx）確定相應(yīng)權(quán)重，根據(jù)公式(12)，對(duì)Kd1與Kd2取加權(quán)平均，得到Kd1與Kd2的加權(quán)平均值，將此加權(quán)平均值作為該激光測(cè)深點(diǎn)最終的漫衰減系數(shù)值

3 實(shí)驗(yàn)與分析

為分析所提算法的性能，文中分別應(yīng)用2013年1月西沙甘泉島(Optech Aquarius)與2014年12月江蘇連云港(Optech CZMIL)兩個(gè)航次的ALB實(shí)驗(yàn)數(shù)據(jù)，兩個(gè)航次的測(cè)區(qū)位置如圖4所示。Aquarius和CZMIL均為全波形ALB系統(tǒng)，如表1所示，二者相應(yīng)的技術(shù)參數(shù)是不同的。數(shù)據(jù)采集期間，由于西沙甘泉島遠(yuǎn)離大陸，人類活動(dòng)較少，相比于連云港沿岸海域，甘泉島海域水質(zhì)好，海水能見(jiàn)度高。

圖4 研究區(qū)域位置Fig.4 Location of the study areas

表1 Aquarius與CZMIL的技術(shù)參數(shù)Table 1 Technical parameters of Aquarius and CZMIL

3.1 波形分解結(jié)果

通過(guò)2.1節(jié)分解算法分別對(duì)西沙甘泉島和江蘇連云港兩個(gè)航次的ALB波形數(shù)據(jù)（西沙甘泉島8 034個(gè)，連云港6 011個(gè)）進(jìn)行處理，分別從兩個(gè)航次的波形處理結(jié)果中隨機(jī)選取一個(gè)進(jìn)行展示，如圖5所示。此外，利用分解結(jié)果與去噪后的波形數(shù)據(jù)對(duì)比計(jì)算得到均方根誤差(Root Mean Square Error, RMSE)和確定系數(shù) (Coefficient of Determination,R2)[14–15]，從而對(duì)文中分解算法性能加以驗(yàn)證，如表2所示。

圖5 波形處理結(jié)果Fig.5 Waveform processing results

通過(guò)對(duì)波形去噪結(jié)果分析可知，采用小波去噪算法可將西沙甘泉島和江蘇連云港兩個(gè)航次的原始回波信號(hào)中的噪聲信號(hào)明顯地濾除，去噪效果理想。連云港屬于高渾濁海域，其ALB波形的水底反射回波較弱，且由于該測(cè)區(qū)水深較淺，水體散射回波與水面、水底反射回波疊加在一起，導(dǎo)致波形分解的難度變大，由波形分解結(jié)果可知，使用文中所提波形分解算法能夠?qū)⑺w散射回波精確地分解出來(lái)，其RMSE、R2的平均值分別為10.019 7、0.994 5。西沙甘泉島周圍海域水質(zhì)較好，ALB最大探測(cè)深度相對(duì)較深，因而運(yùn)用文中分解算法處理甘泉島波形數(shù)據(jù)的精度相對(duì)更高，其RMSE、R2的平均值分別為7.471 0、0.994 7。綜上所述，無(wú)論是水質(zhì)清澈的西沙甘泉島海域，還是較渾濁的江蘇連云港沿岸，采用文中所提基于分層異構(gòu)模型的波形分解算法均能夠?qū)B加在水面反射回波和水底反射回波中的水體散射回波精確地提取出來(lái)，得到水體散射回波所對(duì)應(yīng)的函數(shù)（雙指數(shù)函數(shù)），為進(jìn)一步提取漫衰減系數(shù)提供了可靠保障。

3.2 漫衰減系數(shù)提取結(jié)果

為驗(yàn)證文中漫衰減系數(shù)提取算法的有效性，運(yùn)用Billard B傳統(tǒng)方法[9]與文中算法反演漫衰減系數(shù)，并作對(duì)比分析。文中以包含實(shí)驗(yàn)區(qū)域深水和淺水?dāng)?shù)據(jù)為原則，隨機(jī)在西沙甘泉島、江蘇連云港各選取兩個(gè)區(qū)塊內(nèi)的ALB測(cè)深數(shù)據(jù)，應(yīng)用上述兩種算法進(jìn)行處理，得到漫衰減系數(shù)的平均值如表3。

表2 分解算法性能指標(biāo)Table 2 Performance indexes of waveform fitting algorithm

表3 漫衰減系數(shù)提取結(jié)果Table 3 Extraction results of diffuse attenuation coefficient

圖6 甘泉島西南角樣本區(qū)域空間分布Fig.6 Distribution of the sample area in the Ganquan Island

通過(guò)表3分析可見(jiàn)，運(yùn)用文中算法與傳統(tǒng)方法在西沙甘泉島和江蘇連云港海域所得的漫衰減系數(shù)平均值相近，證明了利用機(jī)載LiDAR測(cè)深波形進(jìn)行漫衰減系數(shù)提取算法的可行性。此外，通過(guò)表4可知[25]，基于文中算法所得漫衰減系數(shù)分別反演連云港沿岸和甘泉島海域的水質(zhì)情況，顯示出了較好的區(qū)分度。

表4 漫衰減系數(shù)與水質(zhì)之間的關(guān)系[25]Table 4 Relationship between diffuse attenuation coefficient and water quality[25]

為進(jìn)一步驗(yàn)證文中算法的可行性，文中選擇西沙甘泉島西南角沿岸作為樣本區(qū)域（圖6a），利用其ALB數(shù)據(jù)進(jìn)行處理，并生成該區(qū)域漫衰減系數(shù)空間分布圖（圖6c）與頻率分布直方圖（圖7）、相應(yīng)水深的空間分布圖（圖6b）。

通過(guò)圖6和圖7可知，運(yùn)用文中算法所得漫衰減系數(shù)的分布規(guī)律與該海域相應(yīng)水深的分布規(guī)律基本保持一致，即當(dāng)甘泉島沿岸水越來(lái)越來(lái)深時(shí)，漫衰減系數(shù)會(huì)越來(lái)越?。辉搮^(qū)塊漫衰減系數(shù)主要集中在0～0.5 m?1，其平均值為 0.215 6 m?1。淺水處由于潮水的作用，海水會(huì)對(duì)海底底質(zhì)產(chǎn)生沖刷，大量泥沙等底質(zhì)在水中懸浮，因此，淺水處漫衰減系數(shù)較高；隨著水深的增加，海底底質(zhì)受潮水的沖刷作用減弱，漫衰減系數(shù)逐漸變小。此外，文中還繪制了連云港區(qū)塊A的漫衰減系數(shù)頻率分布直方圖（圖8），結(jié)果表明，該區(qū)域漫衰減系數(shù)主要集中在0.3～0.7 m?1，其平均值為 0.378 8 m?1。

綜上所述，文中利用ALB測(cè)深水體波形提取漫衰減系數(shù)的算法具有較強(qiáng)可行性和適應(yīng)性，可有效提取每一個(gè)激光測(cè)深點(diǎn)的漫衰減系數(shù)。

4 結(jié)論

結(jié)合ALB技術(shù)高精度、高分辨率、靈活機(jī)動(dòng)、快速高效的特點(diǎn)，提出了一種基于機(jī)載LiDAR測(cè)深水體波形的漫衰減系數(shù)提取方法。通過(guò)分析ALB波形數(shù)據(jù)，針對(duì)水面、水體及水底3部分回波的相應(yīng)特性，構(gòu)建分層異構(gòu)模型，分離得到混疊在水面回波和水底回波之間的水體后向散射回波，進(jìn)而利用所得水體后向散射回波實(shí)現(xiàn)大面積水域漫衰減系數(shù)空間分布的精確獲取，為快速、精細(xì)化分析大面積水域的水體環(huán)境變化提供了一種新的解決方案。