基于光譜分層的淺海水深遙感反演方法

楚森森 ，程亮 , ，程儉，張雪東，劉晉銘

( 1. 南京大學 地理與海洋科學學院，江蘇 南京 210023；2. 江蘇省地理信息技術重點實驗室，江蘇 南京 210023；3. 中國南海研究協同創新中心，江蘇 南京 210023；4. 自然資源部國土衛星遙感應用重點實驗室，江蘇 南京 210023；5. 軍事科學院國防工程研究院，北京 100036)

1 引言

淺海水深數據是航行安全、工程建設、資源開發、海洋救護、生態保護等應用的基礎保障[1-3]。傳統的船載聲吶和機載激光測深方法可獲取高精度淺海水深數據，但是成本昂貴，不適用于大規模、周期性測量，并且無法應用于存在爭議或被他國非法侵占的島礁海域[4-6]。

隨著衛星遙感技術的發展，利用衛星影像反演水深已成為水深測量的重要手段。衛星遙感有不受地理空間約束、費用低等優勢，很好地克服了傳統方法的局限性[7]。衛星遙感水深反演的物理基礎是可見光具有一定穿透水體的能力，根據比爾-朗伯特定律可知，可見光通過水體時因水柱吸收和散射產生衰減，利用光衰減特性建立光在水中的傳輸路程（水深）與衛星影像輻亮度關系[8]。衛星遙感水深反演誤差普遍為1～3 m[9]，雖然精度不能滿足我國海圖測量標準，但遙感水深反演方法可解決特殊海域水深信息有無的難題[10]，因此備受關注。目前已有機構實驗性地出版基于衛星遙感反演的海圖，如英國海道測量局已發布一版基于衛星遙感水深反演的安提瓜島海圖[11]，未來，如何進一步提高衛星遙感水深反演精度仍然是研究重點。

衛星遙感水深反演模型包括解析模型和統計模型。解析模型不需要實測水深數據參與，主要利用水體光傳輸機理建立反演模型[12]，但模型復雜且未知參數較多，導致反演結果不確定性增加，這些問題阻礙了解析模型的推廣應用[13]，簡化模型并提高普適性是當前解析模型的研究熱點[14-16]。不同于解析模型，統計模型基于少量實測水深點標定模型，建立水深與遙感反射率之間統計關系，因具有結構簡潔易操作的優點而被廣泛應用[17-18]。統計模型中比較經典的算法有Stumpf比值算法[19]和Lyzenga多項式算法[20]，雖然兩種算法參數由統計模型的實測水深樣本訓練確立，但其數學模型可通過解析模型結合經驗公式推導出來[19-21]，因此相對于神經網絡、隨機森林等算法，Stumpf比值算法和Lyzenga多項式算法，對實測水深樣本數量和分布要求較低，且可以一定程度上抑制底質異構、水體異質帶來的干擾，算法可靠性與遷移性較強[22]。特別是最近兩年，最新研究表明ICESat-2激光衛星數據可以提供少量的實測水深數據用于標定反演模型[23-24]，解決了統計模型算法需要實測水深數據的局限，因此Stumpf比值算法和Lyzenga多項式算法備受關注，并誕生了大量改進算法[25-30]。這些改進算法可歸納為兩類，一是非線性擬合改進，即從機器學習角度考慮，如將Stumpf波段對數調節因子由1個增加到2個[31]，將一次多項式模型提升為多次多項式模型[32-36]，將線性擬合改進為支持向量機擬合等[37-38]；二是地理自適應改進，即從地理學角度考慮，如顧及底質類型水深反演[39]、水深分段反演[40]、地理分區反演[41]等。

統計模型的特點是無需理解水體光傳輸理論，因此，當前統計模型中經典算法的改進主要從機器學習和地理學角度開展，沒有重視光譜測深性能，光譜的測深極限與適用區間會影響水深反演精度，如果進一步從光譜測深性能角度考慮，統計模型反演算法仍有較大精度提升空間。不同光譜測深極限不同，對于清澈水體，藍光（波長：440～540 nm）最大穿透深度接近30 m、綠光（波長：500～600 nm）最大穿透深度接近15 m、紅光（波長：600～700 nm）最大穿透深度接近5 m[42]，但現有水深反演算法沒有顧及不同光譜測深特性差異，在不同水深處采用一套反演系數，如在紅光可達的淺水區與紅光不可達深水區，均賦予紅光波段相同的反演系數或權重，影響水深反演精度。對此，本文從光譜測深性能角度改進Stumpf比值算法和Lyzenga多項式算法，提出一種基于多光譜遙感影像的光譜分層策略和一種基于光譜分層的淺海水深遙感反演方法，以我國南沙群島長線礁和美屬維京群島巴克島礁為實驗區，與經典Stumpf比值算法和Lyzenga多項式算法進行對比分析，證明了本文方法的有效性。

2 實驗區與數據

實驗區為我國南沙海域長線礁和美屬維爾京群島巴克島（圖1），遙感影像數據為Sentinel-2 L2A級多光譜衛星影像（圖1a，圖1b），長線礁、巴克島的Sentinel-2影像獲取時間分別為2018年5月9日2時35分（UTC）和2019年10月17日14時57分（UTC），影像坐標系統分別為WGS84 UTM 49N和WGS84 UTM 20N，采用的Sentinel-2波段為近紅外、紅、綠、藍波段，空間分辨率為10 m。

實測水深數據包括船載聲吶測深數據和機載激光測深數據，其中，長線礁實測水深數據采集儀器為Odom Hydrotrac II單波束聲吶設備和SONIC 2024多波束聲吶設備，采集時間為2018年1月18-20日，水深基準面為理論深度基準面，水深范圍70 m以淺。巴克島實測水深數據來自美國國家海洋和大氣管理局海岸海洋科學中心公開數據集，采集儀器為Laser Airborne Depth Sounder Mark II 激光雷達設備，采集時間為2011年2月21-22日，激光點間距為3 m×3 m，水深基準面為平均低低潮面，水深范圍50 m以淺。測深數據垂直精度和水平精度均符合第五版國際海道測量標準1等測量精度，即垂直定位不確定度小于[0.52+（0.013×水深）2]1/2m，水平定位不確定度小于（5+0.05×水深）m。經查潮汐表，長線礁和巴克島影像獲取時刻潮高分別為1.48 m和0.38 m，實驗區實測水深已校正為影像獲取時刻瞬時海面水深，并提取0～30 m水深作為本文實驗水深數據（圖1）。水深訓練樣本和測試樣本選用Sentinel-2影像像元均值水深，即對單個像元內部的實測水深數據求均值作為該像元的水深值，長線礁水深像元共計1 127個，巴克島水深像元共計357 931個。

圖1 實驗區與水深數據Fig. 1 Study areas and depth data

3 方法介紹

3.1 經典水深反演算法

3.1.1 Stumpf比值算法

Stumpf比值算法的數學表達式為

式 中，z為 水 深；m1、m0為 模 型 的 回 歸 系 數；R(λ1)、R(λ2)分 別為藍、綠波段的反射率；n為藍、綠波段的固定系數，為保證任何情況下對數值為正數且保證對數比值與深度呈線性關系而設定，本文取n=1 000。

3.1.2 Lyzenga多項式算法

Lyzenga多項式算法的數學表達式為

式中，z為水深；a0、ai為 模型的回歸系數；N為參加反演的波段數量；R(λi)為 波段i的反射率；R∞(λi)為波段i對應的深水區反射率均值，可采用直方圖統計方法，選擇直方圖中像元最小的部分作為深水區像元。

3.2 基于光譜分層的水深反演算法

本文提出基于光譜分層的淺海水深遙感反演方法，首先根據光譜對水體的穿透能力，構建光譜分層策略，提取紅光波段可穿透的水深層為紅光層，綠光波段可穿透但紅光不可達的水深層為綠光層，藍光波段可穿透但紅光與綠光不可達的水深層為藍光層，然后基于光譜分層構建水深反演優化模型，獲取水深反演結果。

3.2.1 光譜水體穿透特性

不同光譜在水中衰減系數存在差異，在清澈水體中，藍光、綠光、紅光、近紅外波段對水體穿透能力依次減弱。圖2以研究區長線礁為例展示了遙感影像反射率與水深的關系，圖2a中近紅外波段反射率較小，且隨著水深增加沒有變化，說明近紅外波段基本無法穿透水體；受水體吸收的影響，圖2b至圖2d反射率均隨著水深增加逐漸降低，其中紅光波段和綠光波段分別在水深大于5 m和15 m上下時，反射率值達到最小且趨于穩定，說明紅光波段和綠光波段最大穿透水深在5 m和15 m上下；而藍光波段在0～30 m水深區間反射率值仍未達到最小值，說明藍光波段最大穿透水深可達30 m上下。由于水體懸浮物含量、溫度、鹽度、密度等因素都會影響光譜水體穿透能力，因此不同區域、不同時間的光譜水體穿透能力可能存在不同，但是藍光、綠光、紅光、近紅外波段之間的水深穿透能力差異必然存在，這種不同光譜對水體穿透特性差異是本文構建光譜分層和提高反演模型精度的關鍵。

圖2 反射率與水深散點圖Fig. 2 Scatter plot of reflection values versus depth values

3.2.2 光譜分層策略

提出一種基于影像本身的無參數（無需輸入任何參數）光譜分層策略。將多光譜遙感影像作為輸入，利用Ostu二值化分割算法[43]，分別對近紅外、紅光、綠光波段進行二值化處理，提取近紅外層、紅光層、綠光層、藍光層（圖3）。具體步驟如下：

圖3 光譜分層示例圖Fig. 3 Spectral stratification example diagram

（1）近紅外層Lnr提取。利用近紅外波段無法穿透水體的特性，通過Ostu算法二值化近紅外波段影像Inr，分割出的前景色區域為近紅外層Lnr，Lnr層記錄了水面高頻噪聲（如耀斑、船只、厚云等）與陸地區域，用于水面噪聲和陸地掩膜。

（2）紅光層Lr提取。利用紅光最大可穿透5 m上下水深的特性，基于多光譜影像中的紅色光譜影像Ir，首先，掩膜掉近紅外層Lnr區域；然后，利用Ostu算法對掩膜后的紅色光譜影像Ir進行二值化處理，其前景色為水下紅光信號可達區域，即紅光層Lr。

（3）綠光層Lg提取。利用綠光最大可穿透15 m上下水深的特性，基于多光譜影像中的綠光波段影像Ig， 首先，掩膜掉近紅外層Lnr、 紅光層Lr區域；然后，利用Ostu算法對掩膜后的綠光波段影像Ig進行二值化處理，其前景色為水下綠光信號可達但紅光不可達區域，即綠光層Lg。

（4）藍光層Lb提取。利用藍光最大可穿透30 m上下水深的特性，基于多光譜影像中的藍光波段影像Ib， 首先，掩膜掉近紅外層Lnr、 紅光層Lr和 綠光層Lg區域；然后，將剩余區域作為藍光層Lb。

3.2.3 分層水深反演

基于光譜分層提取的紅光層Lr、 綠光層Lg、藍光層Lb，建立基于光譜分層的Stumpf比值反演改進算法，和基于光譜分層的Lyzenga多項式反演改進算法。其中基于光譜分層的Stumpf比值反演算法數學公式如下：

式中，z為水深；mr1、mg1、mb1分別是紅光層、綠光層、藍光層中比值模型斜率常數；mr0、mg0、mb0分別是紅光層、綠光層、藍光層中水深0 m時，對數比值的偏移量；n為固定常數，本文取n=1 000；R(λr)、R(λg)、R(λb)分別為紅、綠、藍波段的反射率。此處需注意的是，不同于經典Stumpf比值算法（公式（1））中僅綠、藍波段參與反演，在本文的改進算法中，紅光層水深反演采用了紅、綠波段。

基于光譜分層的Lyzenga多項式反演算法數學公式如下：

式 中，z為 水 深；ar0和ari是 紅 光 層 線 性回 歸 系 數；和是綠光層線性回歸系數；ab0和abi是藍光層線性回歸系數；N為參加反演的波段數量；R(λi)為 波段i的反射率；R∞(λi)為 波段i對應的深水區反射率均值。

3.3 精度評價

為定量評價水深反演結果，采用以下3個定量指標：均方根誤差（Root Mean Square Error, RMSE）、平均絕對誤差（Mean Absolute Error, MAE）、平均相對誤差（Mean Relative Error, MRE）。RMSE、MAE、MRE越小，水深反演效果越好。

式中，ei為反演水深值與對應的實測水深值zi之間的差值；n為實測水深測試樣本個數。

4 研究結果

4.1 水深反演結果對比

本節將提出的基于光譜分層的反演算法與經典Stumpf比值算法和Lyzenga多項式算法對比分析?；诠庾V分層的Stumpf比值算法和Lyzenga多項式算法實驗設置如下：分別在紅光層、綠光層、藍光層隨機選取50個水深像元作為訓練樣本，其余水深像元作為測試樣本。經典Stumpf比值算法和Lyzenga多項式算法設置如下：上述紅光層、綠光層、藍光層3個光譜層選取的總計150個水深像元作為訓練樣本，同樣其余水深像元作為測試樣本。即所有算法共用一套訓練樣本和測試樣本。

在長線礁實驗區，對比分析經典Stumpf比值算法與基于光譜分層的Stumpf比值算法，以及經典Lyzenga多項式算法與基于光譜分層的Lyzenga多項式算法的水深反演結果（圖4）。從圖4來看，基于光譜分層的反演算法優于經典算法，相對于經典算法反演結果（圖4a，圖4c），基于光譜分層算法的反演結果在藍色深水區斑點噪聲更少，且紅色淺水區域礁盤輪廓更加清晰（圖4b，圖4d）。圖5展示了以上反演結果分別在紅光層、綠光層、藍光層和全局的誤差對比，從圖5整體來看，相對于經典算法反演誤差圖（圖5a1至圖5a4，圖5c1至圖5c4），基于光譜分層算法的反演水深與實測水深誤差更聚集于0值（圖5b1至圖5b4，圖5d1至圖5d4），說明誤差更小。圖5中每個子圖右上角誤差指標可以定量地看出，基于光譜分層的反演算法明顯優于經典算法，相對于經典Stumpf比值算法，基于光譜分層的Stumpf比值算法RMSE、MAE、MRE分別降低0.89 m、0.65 m、19%（圖5a4，圖5b4），從紅光層、綠光層、藍光層來看，誤差均有降低，其中紅光層MRE降低了149%（圖5a1，圖5b1）；相對于經典Lyzenga多項式算法，基于光譜分層的Lyzenga多項 式 算 法RMSE、MAE、MRE分 別 降 低0.64 m、0.35 m、8%（圖5c4，圖5d4），從紅光層、綠光層、藍光層來看，誤差同樣均有降低，其中紅光層MRE降低了90%（圖5c1，圖5d1），改進效果顯著。

圖4 長線礁水深反演結果圖Fig. 4 Bathymetric maps of the Changxian Reef produced by inversion

圖5 長線礁反演水深與實測水深值誤差分布Fig. 5 Distribution of differences obtained by subtracting the measured values from the inversion water depth for the Changxian Reef

在巴克島實驗區，圖6對比分析了經典Stumpf比值算法與基于光譜分層的Stumpf比值算法，以及經典Lyzenga多項式算法與基于光譜分層的Lyzenga多項式算法的水深反演結果。對比巴克島激光測深地形圖（圖1c），從巴克島右側深水區域來看，基于光譜分層的反演算法斑點噪聲較少（圖6b，圖6d），優于對應的經典算法反演結果（圖6a，圖6c）。圖7展示了以上反演結果分別在紅光層、綠光層、藍光層和全局的誤差對比，從每個子圖右上角誤差指標可以定量地看出，基于光譜分層的反演算法同樣明顯優于經典算法，相對于經典Stumpf比值算法，基于光譜分層的Stumpf比值算法RMSE、MAE、MRE分別降低0.41 m、0.36 m、4%（圖7a4，圖7b4），從紅光層、綠光層、藍光層來看，誤差均有降低，其中紅光層MRE降低了58%（圖7a1，圖7b1）；相對于經典Lyzenga多項式算法，基于光譜分層的Lyzenga多項式算法RMSE、MAE、MRE分別降低0.45 m、0.52 m、8%（圖7c4，圖7d4），從紅光層、綠光層、藍光層來看，誤差均有降低，其中紅光層MRE降低了141%（圖7c1，圖7d1），改進效果顯著。

圖6 巴克島水深反演結果圖Fig. 6 Bathymetric maps of the Buck Island Reef produced by inversion

圖7 巴克島反演水深與實測水深值誤差分布圖Fig. 7 Distribution of differences obtained by subtracting the measured values from the inversion water depth for the Buck Island Reef

4.2 訓練樣本數量對反演精度影響

為對比分析訓練樣本數量N對反演算法的影響，設置不同數量的訓練樣本，即在長線礁實驗區為15～180個，巴克島實驗區為15～600個，以15為間隔依次取值，并對于每個訓練樣本數量N，做100次水深反演重復實驗，對100次重復實驗誤差值求平均作為訓練樣本數量為N時的反演精度，其中每次樣本取值規則為：分別在紅光層、綠光層、藍光層隨機等量取值，如訓練樣本數量為15時，則分別在紅光層、綠光層、藍光層隨機取5個訓練樣本。

長線礁和巴克島實驗區訓練樣本數量對反演精度的影響如圖8和圖9所示，圖中誤差棒代表100次重復實驗誤差值的標準差，當標準差較大時，以數字表示，數字所在橫軸位置表示對應的樣本數量，數字顏色表示對應的算法。從圖8和圖9來看，隨著訓練樣本數量增加，RMSE、MAE、MRE整體上呈現先降低后趨于穩定的趨勢，訓練樣本數量N＜60時，誤差呈現降低趨勢，其中經典Stumpf比值算法和Lyzenga多項式算法比基于光譜分層的反演算法可以更快地實現誤差穩定；訓練樣本數量N≥60時，誤差基本穩定且誤差棒較小，該階段中基于光譜分層的Stumpf比值算法和Lyzenga多項式算法誤差低于相應的經典算法，表明訓練樣本數量N≥60時，可充分發揮基于光譜分層反演算法優勢。

圖8 長線礁實驗區訓練樣本數量對反演精度的影響Fig. 8 Effects of the number of training samples on the accuracy of bathymetric inversion for the Changxian Reef experimental area

5 討論

相對于經典Stumpf比值算法和Lyzenga多項式算法，基于光譜分層的水深反演算法擁有更高的反演精度，主要原因是基于光譜分層的水深反演算法顧及到不同光譜層的波段特性差異，紅、綠、藍波段在不同光譜層中采用了最適宜的反演系數（權重），從而取得較高的反演精度。即水深反演常用紅、綠、藍波段的測深性能存在差異，不同光譜層各波段貢獻程度不一樣，在紅光層，因紅光在水中衰減強于藍、綠光，所以紅光對水深區分度高于藍、綠光，因此，紅光層中紅光對水深反演貢獻高于藍光和綠光；在綠光層，因紅光不可達，所以紅光對水深反演貢獻極小，因此，綠光層中藍光和綠光對水深反演貢獻最高；在藍光層，因紅光和綠光不可達，所以紅光和綠光對水深反演貢獻極小，因此，藍光層中藍光對水深反演貢獻最高。經典Stumpf比值算法和Lyzenga多項式算法全局采用一套反演系數，無法兼顧不同光譜層波段對水深反演的貢獻差異，而基于光譜分層的水深反演算法根據紅、綠、藍波段極限測深能力劃分光譜層，分別建立最適宜反演系數（權重），解決了上述問題，提高了水深反演性能。

從紅光層MRE來看，相對于經典Stumpf比值算法和Lyzenga多項式算法，基于光譜分層的水深反演算法精度有巨大提升，MRE降低了58%～149%。主要原因如下：紅光層水深較淺，極淺水域MRE反演誤差的控制是當前水深反演難點，微小的水深反演偏差可引起較大的MRE，如實際水深為0.5 m的區域被反演為水深2 m，其MRE為300%。而經典Stumpf比值算法和Lyzenga多項式算法在水深較淺的紅光層仍然為米級誤差，造成較大的MRE?；诠庾V分層的水深反演算法在水深較淺的紅光層建立獨立的反演系數，充分發揮紅光波段優異的淺水測深性能，將紅光層誤差控制在分米級，大幅降低了淺水區MRE。

圖8和圖9顯示隨著訓練樣本數量的增加，水深反演誤差整體先降低后趨于穩定，這是因為隨著訓練樣本數量的增加，訓練樣本規模上可以映射出更多光譜到水深的特征規律，因此誤差降低，當訓練樣本呈現的特征被掌握后，訓練樣本繼續增加不會呈現新的特征規律，因此趨于穩定。其中經典Stumpf比值算法和Lyzenga多項式算法比基于光譜分層反演算法可以更快地實現誤差穩定，這是因為基于光譜分層的水深反演算法對樣本數量有更高的需求，即相同數量樣本情況下，經典Stumpf比值算法和Lyzenga多項式算法可以使用全部樣本訓練出全局反演系數，而基于光譜分層的水深反演算法需要拆分為3份，分別訓練對應紅光層、綠光層、藍光層的反演系數，需要更多的訓練樣本數據。但是隨著訓練樣本數量增加，基于光譜分層的反演算法的反演誤差迅速降低，當訓練樣本數量N≥60時，其誤差可低于相應的經典算法，并實現基本穩定且誤差棒較小，證明基于光譜分層的反演算法對訓練樣本數量需求并未增加太多，從水深反演應用案例來看[22-25,30-31]，N≥60的要求不會成為本文算法推廣應用的障礙，訓練樣本數量N≥60可作為其他實驗參考值。

對比基于光譜分層的Stumpf比值算法和Lyzenga多項式算法，基于光譜分層的Stumpf比值算法繼承了經典Stumpf比值算法的優點，即一定程度上抑制水體不均與底質異構引起的干擾，同時也削弱太陽高度角、水面波動及衛星姿態、掃描角變化等差異對水深反演的影響[9]?；诠庾V分層的Lyzenga多項式算法繼承了經典Lyzenga多項式算法的優點，同樣可以一定程度上抑制水體不均與底質異構問題，但相對于僅有兩個波段參與運算的Stumpf比值算法，Lyzenga多項式算法可以綜合多個波段的水深信息，這也是本文實驗Lyzenga多項式算法精度優于Stumpf比值算法的一個重要因素（圖5, 圖7至圖9）?？偟膩碚f，兩種算法各有優勢，結合實驗結果與算法本身特點來看，當研究區較小且影像成像質量較佳時，即太陽高度角、水面波動及衛星姿態、掃描角變化等影響較弱時，優先推薦基于光譜分層的Lyzenga多項式算法；當研究區較大或需要模型遷移時，有文獻表明Stumpf比值算法可以取得更優的反演精度[22]。但是由于水深反演干擾因素多且海洋環境復雜，兩種算法對比與適用場景值得進一步研究。

圖9 巴克島實驗區訓練樣本數量對反演精度的影響Fig. 9 Effects of the number of training samples on the accuracy of bathymetric inversion the Buck Island Reef experimental area

6 結論

本文在經典Stumpf比值算法和Lyzenga多項式算法基礎上，從光譜測深性能角度，提出一種基于光譜分層的水深遙感反演優化方法。首先根據光譜對水體的穿透能力，提出一種基于多光譜本身的無參數光譜分層策略，提取紅光層、綠光層、藍光層；然后分光譜層構建水深反演模型，建立各光譜層最適宜波段和其反演系數，解決現有水深反演算法全局一套反演系數難以適應不同光譜層的局限，最終達到提高水深反演精度的目的。以長線礁和巴克島作為實驗區開展驗證分析，得出以下結論：

（1）提出的基于光譜分層的水深反演算法精度高于經典Stumpf比值算法和Lyzenga多項式算法：相對于經典算法，基于光譜分層的水深反演算法RMSE、MAE、MRE可 降 低0.41～0.89 m、0.35～0.65 m、4%～19%，證明了本文算法的有效性。

（2）在紅光層，即水深較淺區域，基于光譜分層的水深反演算法MRE大幅降低：相對于經典算法，基于光譜分層的水深反演算法在紅光層MRE降低58%～149%，精度提升明顯。

（3）基于光譜分層的水深反演算法對樣本數量與分布要求較高，但構不成算法推廣應用的障礙：訓練樣本數量N≥60，即紅光層、綠光層、藍光層訓練樣本數量均大于等于20，可作為其他實驗參考值。

本文從光譜測深性能角度改進了經典Stumpf比值算法和Lyzenga多項式算法并取得較好反演效果，但本文算法仍有較大改進空間。如所提光譜分層反演策略理論上可以與現有機器學習角度、地理學角度改進算法疊加使用；如針對更多光譜的高光譜衛星數據時，通過光譜優選等方法建立更細的光譜分層策略，這些策略可進一步提升反演精度，值得未來研究。