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

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

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

1 引言

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

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

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

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

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

2.1 基于分層異構模型的波形分解算法

機載LiDAR測深系統利用具有較強透水能力的藍綠色波段(532 nm)激光進行探測，其激光接收器能夠以數字化的形式記錄每個激光脈沖的全部回波波形信號（圖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

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

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

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

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

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

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

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

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

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

2.2 漫衰減系數Kd的提取

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

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

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

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

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

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

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

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

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

3 實驗與分析

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

圖4 研究區域位置Fig.4 Location of the study areas

表1 Aquarius與CZMIL的技術參數Table 1 Technical parameters of Aquarius and CZMIL

3.1 波形分解結果

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

圖5 波形處理結果Fig.5 Waveform processing results

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

3.2 漫衰減系數提取結果

為驗證文中漫衰減系數提取算法的有效性，運用Billard B傳統方法[9]與文中算法反演漫衰減系數，并作對比分析。文中以包含實驗區域深水和淺水數據為原則，隨機在西沙甘泉島、江蘇連云港各選取兩個區塊內的ALB測深數據，應用上述兩種算法進行處理，得到漫衰減系數的平均值如表3。

表2 分解算法性能指標Table 2 Performance indexes of waveform fitting algorithm

表3 漫衰減系數提取結果Table 3 Extraction results of diffuse attenuation coefficient

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

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

表4 漫衰減系數與水質之間的關系[25]Table 4 Relationship between diffuse attenuation coefficient and water quality[25]

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

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

綜上所述，文中利用ALB測深水體波形提取漫衰減系數的算法具有較強可行性和適應性，可有效提取每一個激光測深點的漫衰減系數。

4 結論

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