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

亓超，周豐年，吳敬文，宿殿鵬，王賢昆，陽凡林*

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

1 引言

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

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

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

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

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

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

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

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

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

圖2 機載LiDAR測深波形示意圖Fig.2 Airborne LiDAR bathymetric waveform

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由于近岸水深通常小于30 m，與飛機航高H（通常為400～500 m）相比，參數(shù)Di對W的影響可以忽略不計，因此，W可認(rèn)為是一個常數(shù)。但是W涉及未知參數(shù)太多，不能直接確定。通過分析發(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)）與簡化后的水體散射回波強度方程(公式(7))構(gòu)建模型，從而計算漫衰減系數(shù)Kd值。以雙指數(shù)函數(shù)（圖3）BC段為例，當(dāng)bx≤t

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

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

3 實驗與分析

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

圖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é)果

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

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

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

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

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

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

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

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

通過圖6和圖7可知，運用文中算法所得漫衰減系數(shù)的分布規(guī)律與該海域相應(yīng)水深的分布規(guī)律基本保持一致，即當(dāng)甘泉島沿岸水越來越來深時，漫衰減系數(shù)會越來越小；該區(qū)塊漫衰減系數(shù)主要集中在0～0.5 m?1，其平均值為 0.215 6 m?1。淺水處由于潮水的作用，海水會對海底底質(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測深水體波形提取漫衰減系數(shù)的算法具有較強可行性和適應(yīng)性，可有效提取每一個激光測深點的漫衰減系數(shù)。

4 結(jié)論

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