王鑫，貝祎軒，陳卓，張凱,2*

( 1. 山東科技大學 測繪與空間信息學院，山東 青島 266590；2. 自然資源部第二海洋研究所，浙江 杭州 310012；3. 中兵勘察設計研究院有限公司，北京 100053；4. 廣州三海海洋工程勘察設計有限公司，廣東 廣州 510220)

1 引言

隨著我國海洋強國戰(zhàn)略的逐步實施，近岸水深數據的需求顯著增加。因此，如何高效、準確地獲取近岸水深信息成為當前研究熱點。目前，水深測量主要依靠船載聲吶測量[1–2]與機載激光雷達測量[3–4]。其中，船載多波束測深系統(tǒng)橫向覆蓋寬度達深度的3～4 倍，可高效獲取海底的準確水深信息，但在淺水區(qū)效率較低且存在擱淺風險。機載激光雷達測量系統(tǒng)可高效測量淺水區(qū)域水深信息，但其廣泛應用受到測量成本高和飛行區(qū)域合法性的限制。相比之下，利用遙感衛(wèi)星的多光譜影像反演近岸水深信息[5–6]，因其高效、高時空分辨率、低成本的特點，表現出顯著的應用價值。

自20 世紀80 年代Lyzenga[7]提出雙層流理論水深反演模型以來，遙感衛(wèi)星技術的進步有力推動了多光譜遙感影像反演水深理論的發(fā)展[8]，所形成的反演方法大致可分為3 類：理論解釋模型、半理論半經驗模型、統(tǒng)計模型。其中，理論解釋模型基于水光場輻射傳輸方程[7,9–11]，通過建立光學傳感器接收到的輻亮度與水深和底質反射的解析表達式來反演水深。這類方法理論嚴密，但所需水體光學參數較多且獲取困難，模型構造復雜。而半理論半經驗模型通過結合理論模型和經驗參數來實現水深反演，其中以Stumpf對數比值[12]為代表的方法，將輻射亮度表示為深水區(qū)與海底反射輻亮度之和，建立其與水深之間的解析表達式來預測水深，得益于模型構造簡單與物理機制清晰的優(yōu)勢被廣泛應用[13–16]。然而對數比值法特征維度較少，導致反演精度受到模型表達能力的限制[17]。近年來，隨著計算機技術的迅速發(fā)展，基于機器學習算法的統(tǒng)計模型逐漸成為水深反演最熱門的前沿研究領域[18–21]。該方法不需要考慮水深遙感的物理機制，而是從水深與圖像輻亮度值之間的統(tǒng)計關系出發(fā)，自發(fā)學習數據之間的本質聯系來建立模型。憑借其在解決多變量、非線性復雜問題等方面的優(yōu)勢，被引入統(tǒng)計水深反演模型中，取得了優(yōu)于對數比值法的反演效果。然而受衛(wèi)星影像信噪比低和實測數據量少的制約，多光譜水深反演精度提升有限，成為限制該方法廣泛應用的瓶頸。因此，通過深入挖掘數據潛在信息提高水深反演精度，是多光譜水深反演方法研究的一個發(fā)展方向。

上述方法通過挖掘影像光譜值與水深之間的函數關系與統(tǒng)計互相關特征進行水深反演。除此之外，水深信息本身的空間統(tǒng)計自相關特征亦可用于估計未采樣點位上的水深信息，如基于稀疏水深數據進行空間插值（如反距離加權法[22]、克里金法[23–24]）估計未采樣位置的水深信息。因此，在上述傳統(tǒng)水深反演方法基礎上，通過挖掘水深信息的空間自相關特征，有望有效提高水深反演精度。對此，Su 等[25]提出回歸克里金法（Regression Kriging, RK），將空間相關性特征納入了Stumpf 對數比值模型中，有效提升了水深反演精度。然而，RK 和Stumpf 模型皆屬于線性模型，在描述復雜非線性映射關系時表達能力有限。同時，RK 基于平穩(wěn)假設的前提在實際觀測環(huán)境中也常常難以得到有效滿足[26]。因此，該方法的效果在環(huán)境各向異性顯著的區(qū)域明顯下降。

針對上述問題，本文以提高多光譜水深反演精度為目標，利用機器學習在解決多變量與非線性問題方面的優(yōu)勢，通過融合數據源的空間自相關特征與統(tǒng)計互相關性來研究高精度水深反演方法，以期在現有研究基礎上，進一步提高預測水深的準確性。

2 方法

2.1 對數比值模型（Stumpf 模型）

對數比值法是代表性的半理論半經驗水深反演模型[12]，該方法以藍、綠波段反射率的對數比值作為反演因子，其模型表達式為

式中，d為反演出的水深；β0和β1是通過回歸方程得到的系數；n為縮放比例因子；Rw(λb)和Rw(λg)分別是藍、綠波段的反射率。

該方法具有構造簡單，對環(huán)境因素干擾（如大氣、水體和海底反射率的變化）不敏感等優(yōu)勢。但作為線性擬合模型，該方法在實際應用中仍然受到復雜水質和海底底質變化的顯著影響，難以準確描述復雜海區(qū)反射率與水深呈現出的非線性關系[27]。同時，未能綜合利用多個波段的光譜測量信息亦是該方法的一個缺點。

2.2 隨機森林

機器學習方法可以對不同變量之間的復雜映射機制進行有效學習和表達。其中，隨機森林（Random Forest，RF）模型憑借其優(yōu)秀的非線性回歸和泛化能力、多變量特征融合能力以及運算速度快等優(yōu)勢[28–30]，得到廣泛應用。用于多光譜反演水深時，RF 模型通過并行處理策略將一組回歸決策樹組合生成學習器，并利用未抽取樣本來泛化模型誤差，綜合分析每個模型的預測結果。通過上述學習過程，RF 模型可以獲得優(yōu)于對數比值法的水深反演精度，其算法流程如圖1 所示。

圖1 隨機森林算法示意圖Fig. 1 Schematic diagram of random forest

光學水深遙感反演的物理基礎是光對水體的穿透能力，受水體的漫衰減影響，光在水體中的穿透深度與水體特性、不同譜段的漫衰減系數密切相關性，可通過機器學習發(fā)掘其中顯著而復雜的函數關系。RF 模型的數學表達式為

式中，h（ ）代表隨機森林算法的函數關系式；N為決策樹的總數；輸入向量X為波段光譜值；Θk為樣本通過模型訓練得到的第k棵樹的參數向量；d為綜合所有決策樹預測結果的平均值所得到的反演水深。

2.3 空間自相關隨機森林

2.3.1 空間自相關隨機森林模型

傳統(tǒng)RF 模型通過分析多光譜數據與實測水深數據之間的統(tǒng)計互相關特征預測水深，但忽視了數據自身的空間自相關特征，故未能充分利用有效信息進行學習和表達[31]。為此，本文在訓練模型過程中，以實測水深及對應的波段光譜值作為基礎因/自變量前提下，將待測點周圍具有高度空間自相關性的已知點所對應的屬性值作為自相關變量納入RF 的學習框架中，以此實現統(tǒng)計互相關特征與空間自相關特征的融合。其算法流程圖如圖2 所示。

圖2 空間自相關隨機森林算法示意圖Fig. 2 Schematic diagram of spatial autocorrelation random forest

2.3.2 空間自相關與最佳空間間隔

空間自相關特征在遙感圖像分析中的應用十分廣泛，常用于變化檢測與分類優(yōu)化[32–33]?？臻g自相關性指的是研究區(qū)域內要素的屬性值之間潛在的相互依賴性，可用來衡量地理現象空間聚集程度。研究表明，預測點周圍具有較高自相關性的已知點所包含的屬性信息可用來提高預測的精度[34]，理論上，數據之間的統(tǒng)計相關性越強，預測精度提升幅度越大。對此，全局莫蘭指數（Moran’s I）[35]常用于度量研究數據的聚集性于自相關程度。莫蘭指數越大，表明研究對象的聚集性越強，自相關程度越高。其計算公式為

式中，N為要素的總數目；xi和xj分別為第i和第j個要素的光譜值與水深值；x為上述屬性值的均值；wij為要素i和j的空間權重值，這里為基于距離的鄰接矩陣。

由于Moran’s I 基于空間隨機性分布的零假設，因此需要利用Z得分和P值來判斷假設是否成立，并檢驗自相關指數的顯著性[36]。其中，P值是由已知分布的曲線得出的面積近似值，而Z值按以下形式計算：

式中，E[I]和V[I]分別為Moran’s I 指數的期望值和方差。

受空間距離的影響，變量的自相關程度僅在一定范圍內較為顯著，若超過該距離，則認為發(fā)生了空間變異。為確定局部最佳自相關的空間間隔，采用半方差來體現空間變異程度[37]，半方差越大，自相關越弱，其計算公式為

式中，h為樣本點之間的距離；n為距離h范圍內成對樣本點的數量；Z為樣本點所對應的不同波段光譜值與水深值。

2.3.3 空間自相關變量

為了獲取訓練樣本的自相關特征信息，需要建立局域子窗進行提取，具體步驟如下：

（1）計算訓練樣本的水深與不同波段值的半方差函數和莫蘭指數，用于確定搜索窗口初始尺寸。

（2）根據計算窗口中待測點到m個已知點的二維歐式距離進行排序，得到衡量自相關程度的距離向量Dn，并將遍歷所有待測點得到的n個距離向量，合成為空間自相關距離矩陣（EDF）。

（3）由于訓練數據分布的不均勻性，Dn中元素數量k應取所有距離向量的最少元素數目，且需保證不等于0，以滿足向量組成矩陣的條件；此外，為在不影響預測精度條件下，提升模型訓練效率，還需對搜索窗口大小進行調整并重復步驟（2），將得到的EDF作為表征數據空間自相關特征的補充變量，EDF的具體計算公式計算為

2.3.4 空間自相關隨機森林模型訓練/驗證

基于上述步驟，空間自相關隨機森林（Spatial Autocorrelation Random Forest，SARF）模型訓練與驗證過程存在一定差異，下面列出關鍵步驟說明。

（1）對水深點個數為a的訓練數據而言，所用的水深點既充當待測點，也作為已知點參與模型訓練，即按給定的搜索半徑計算每個當前水深點到周圍鄰近點的距離，用以構成大小為a×k的自相關矩陣EDF。將當前點水深值作為監(jiān)督學習的因變量，與空間自相關變量、基礎自變量輸入RF 模型中訓練，得到訓練參數。

（2）對于驗證數據（數量為b）來說，每個水深點都為待測點，故需將訓練數據中的所有水深點作為已知點來構建自相關矩陣（大小為b×k），并將待測點對應的空間自相關變量與基礎自變量輸入訓練好的模型中，即可求得待測點水深。

3 實驗區(qū)域與數據源

3.1 實驗區(qū)域概況

本文研究區(qū)域為北島周邊淺海水域。北島位于南海西沙群島海域七連嶼的中部，面積約為3.6 km2，島礁呈長條形（圖3a）。作為典型的珊瑚島，其底質由珊瑚、砂和貝屑組成，水下地形特征復雜。由于遠離大陸且受人類活動影響小，島嶼周邊水質清澈，適合多光譜水深反演。式中，Dn為距離向量；EDF為n個距離向量的集合。

圖3 北島地理位置（a）及原位水深點測量分布（b）Fig. 3 Location of Beidao (a) and distribution of in situ depth measurements (b)

以每個距離向量中k個元素為前提，按其自相關程度的強弱順序，依次獲取每個已知點所對應的屬性值并將其與EDF共同組成空間自相關變量X2，與基礎變量X1構成反演因子。

式中，r1、r2、r3為預測點對應的藍、綠、紅波段光譜值比值，即Bblue/Bgreen，Bgreen/Bred，Bblue/Bred；r1k、r2k、r3k、zk分別為k個已知點的光譜值比值與水深值。

3.2 數據源

本研究使用的原始影像為高分辨率的World-View-2 衛(wèi)星影像（圖3b），拍攝于2017 年3 月11 日11 時。多光譜分辨率為1.84 m，包含藍、綠、紅以及近紅外4 個波段。使用的數據來自于2013 年4 月測量得到的北島1:2 000 陸域及水下地形測量圖，對其進行矢量化與篩選后得到1 700 個水深點（圖3b）。

3.3 數據預處理

數據的獲取過程中存在時間、空間以及儀器測量帶來的誤差，因此在利用遙感影像與水深數據進行定量水深研究前需要對這兩類數據源進行預處理。對于Worldview-2 影像，首先通過輻射定標和大氣校正得到真實地物反射率，并對影像進行幾何校正以及云霧和陸地部分的掩膜，使用Hedley 法[38]消除受到太陽光、海面波浪和光線入射角度等因素造成的太陽耀斑。對于水深數據而言，將其坐標轉換到與影像同一坐標系UTM WGS-84 下，并進行地理配準，以提取水深點對應的像元值。由于遙感影像采用的是衛(wèi)星過境時瞬時海面深度，而實測水深數據是以理論深度基準面為基準的穩(wěn)態(tài)水深，因此還需要根據潮位預報表提供的潮高進行潮汐改正。

本文中水深分布區(qū)如圖4 所示，主要集中于0～6 m淺水區(qū)域。為模擬小樣本淺海水深反演這一常見應用場景，同時確保數據具有代表性，以1 m 為間隔分層抽樣，共選取150 個水深數據作為訓練樣本，而剩余1 550 個水深點則作為驗證樣本。為定量比較4 種反演模型效果，本文選取決定系數（R2），均方根誤差（RMSE）和平均絕對誤差（MAE）進行精度評價。

圖4 水深點分布區(qū)間與數量Fig. 4 Distribution interval and number of water depth points

4 實驗結果分析

4.1 變量的空間自相關性

通過計算訓練樣本的Moran’s I 指數、Z值和P值對水深與波段光譜值的空間自相關性加以評估（表1）。以藍光波段為例，Moran’s I 指數為0.631，說明該波段光譜值呈現出較強的聚類型式，而其Z值（10.967）和P值（0.001）則表明可以拒絕空間隨機性分布的零假設，并通過了置信度為99%的顯著性檢驗。表1 中結果顯示，訓練數據的水深、藍、綠波段皆具有顯著的自相關性，而紅光波段自相關性較弱。

表1 研究區(qū)變量的全局莫蘭指數、Z 值和P 值Table 1 Global Moran’s I, normalized Z value and P value of variables in the study area

4.2 搜索窗口尺寸確定

研究區(qū)域變量的全局Moran’s I 及半方差隨空間間隔的變化如圖5 所示，通過計算波段值的半方差，并除以最大值歸一化后與全局Moran’s I 對比發(fā)現，隨著空間間隔變大，全局Moran’s I 減小，而半方差值呈上升趨勢，這表明變量的空間變異程度增加，而空間自相關性逐漸減小。值得注意的是，二者相交所產生的交點反映了空間自相關與空間變異的平衡。對于水深及藍、綠、紅波段而言，對應的最佳空間間隔分別為175 m、156 m、152 m 和76 m，這表明上述變量在該局域內具有顯著空間自相關性?；谏鲜鲇^察，本文將初始搜索窗口尺寸定為175 m，并根據距離向量的最小元素數原則，將搜索窗口尺寸擴大為220 m，用以獲取空間自相關變量并構建SARF 水深反演模型。

圖5 全局的全局莫蘭指數及半方差的變化Fig. 5 Global Moran’s I and semivariance with different lags of variable

4.3 精度評價

為驗證實驗方法效果，本文首先將反演結果與傳統(tǒng)方法進行了對比（表2）?？梢钥吹?，受限于模型表達能力，Stumpf 模型反演精度最低，RMSE 為2.067 m，R2僅為0.797，相關性也最低。普通克里金（Ordinary Kriging, OK）模型受限于空間插值方法對樣本點密度等要求，在小樣本實測水深數據集的條件下，其預測精度亦不甚理想（RMSE 為1.894 m）。RF 模型得益于綜合多個波段光譜信息，RMSE 為1.635 m，決定系數達到0.888，相比前兩種模型預測結果更準確。在這4 種模型中，SARF 模型表現效果最好，RMSE 僅為1.338 m，相對于RF 模型預測精度提升18%，同時相關性參數也提高到0.923。

表2 驗證數據測試精度對比Table 2 Accuracy comparison of different methods

4.3.1 誤差分析

為了比較4 種模型的反演精度和優(yōu)劣性，根據水深反演圖（圖6）對不同方法的預測結果與真實水深進行對比分析。圖6 中黑線是斜率為1 的指標線，黑線上方表示預測水深值大于實測水深值，下方則表示預測值偏小；紅線是通過最小二乘擬合的趨勢線，用以衡量與指標線的擬合程度；為體現預測值聚集區(qū)間，采用色棒來標識，聚集程度越高，則數值越大且顏色越淺。

Stumpf 模型的水深反演效果如圖6a 所示，在水深淺于10 m 的區(qū)域，反演結果與實測水深符合較好。但是受限于線性表達方式，出現反演水深值小于0 的不合理情況，且對于深于10 m 的數據，水深反演結果產生了明顯的偏差。很顯然，該偏差現象產生的根源在于線性模型無法準確擬合數據的非線性特征；類似地，作為非線性局部擬合方法，OK 模型較為準確地預測了0～8 m 區(qū)域水深（圖6b），而受制于8 m以深實測數據量的減少，無法真實還原深水區(qū)情況。雖然該方法顧及了水深的自相關特征信息，但未考慮水深與光譜值之間的相關性，因此只能通過增加實測數據量來彌補模型的缺陷。對比前兩種模型，圖6c和圖6d 中機器學習預測精度提升效果明顯，較為完整地還原了研究區(qū)域內真實水深分布情況。但從細節(jié)上比較發(fā)現，RF 模型受訓練過程中低信噪比與實測數據量不足影響，指標線上方出現大量離散點且與趨勢線相距較遠，并影響了預測結果可靠性；與之相比，SARF 模型具有明顯的抗離群值效果以及更高的反演精度。

圖6 對數比值（a）、普通克里金（b）、隨機森林（c）和空間自相關隨機森林（d）模型反演水深散點圖（訓練數量為150）Fig. 6 Scatter diagram of predicted depth values by Stumpf (a), ordinary Kriging (b), random forest (c), and spatial autocorrelation random forest (d) models (the number of training data points is 150)

殘差分布圖分別給出了不同方法在表達水深-波段值耦合機制過程中的差異，從而通過殘差的分布形式來對4 種模型的細節(jié)差異進行比較（圖7）。圖7a顯示，Stumpf 模型的殘差在10 m 以深區(qū)域呈現出明顯的遞增趨勢，預測值與真實值存在較大差異；由圖7b 可知，OK 模型整體殘差值較大，特別是在8 m 以深區(qū)間，大量的離散點明確顯示出模型自身的缺陷；而RF 模型在淺水區(qū)域的效果較好（圖7c），殘差值在0 兩側均勻分布，但在水深超過10 m 的區(qū)域，其預測結果受到離群值干擾；與前者相比，無論是從殘差整體分布形式，亦或是深水區(qū)數據受自身統(tǒng)計值波動影響所產生的離群點數量，SARF 模型都表現出更明顯的優(yōu)勢（圖7d）。

圖7 對數比值（a）、普通克里金（b）、隨機森林（c）和空間自相關隨機森林（d）模型的殘差散點圖（訓練數量為150）Fig. 7 Scatter diagram of residual error by Stumpf (a), ordinary Kriging (b), random forest (c), and spatial autocorrelation random forest(d) models (the number of training data points is 150)

4.3.2 不同訓練樣本數量下反演結果的精度對比

訓練樣本數量的變化可顯著影響機器學習方法的預測效果。通過計算訓練樣本數為500 時SARF 模型和RF 模型的殘差標準差（圖8a）可以看到，RF 模型殘差標準差（1.295）遠高于SARF 模型（1.012），因此后者殘差離散程度較小，分布更為集中；其次，通過觀察0.05 顯著性水平下殘差分布區(qū)間發(fā)現，RF 模型殘差主要分布于[1.25, 1.35]，而SARF 模型殘差分布區(qū)間集中在[0.97, 1.05]，后者的殘差分布幅度更小，說明預測值與實測值擬合效果更好。除此之外，通過觀察訓練樣本占全部樣本比例下的均方根誤差（圖8b）發(fā)現，SARF 模型的誤差大小及其下降速率始終小于RF 模型，且隨著訓練數據量的增加，該優(yōu)勢愈加明顯；當訓練樣本占比達到60%時，相比于RF 模型，SARF 模型的誤差降低約27%，該結果進一步驗證了將數據源的空間自相關特征納入機器學習框架中，可顯著提升預測結果的準確性。

圖8 訓練樣本量500 時殘差分布直方圖（a）和不同訓練數據占比下的均方根誤差（b）Fig. 8 The histogram of error distribution when the training sample is 500 (a) and root mean square error for different training data shares (b)

本文基于SARF 模型使用350 個訓練數據得到北島海域反演水深，并將其與實測陸上高程點生成反演水深與陸域地形圖（圖9a）。可以看到，珊瑚島周圍海底地形分布表現出復雜的特征。近岸淺水區(qū)域海底主要由沙質沉積物構成，水深變化較為平衡。而在離岸較遠的區(qū)域，受海流的影響，島嶼南北兩側呈現出截然不同的海底地貌細節(jié)特征（圖9b）。顯然，海底地形信息表現出較強的空間相關性，從而驗證了利用空間相關性輔助水深反演的有效性。

圖9 反演的陸域與水下地形圖（a）和地貌細節(jié)圖（b）Fig. 9 Bathymetry retrieval of onshore and inversion bathymetric topographic map (a) and geomorphic details (b)

5 結論與討論

隨著多光譜衛(wèi)星反演水深技術的發(fā)展，如何在小樣本現場實測數據的前提下，有效提高預測精度一直是研究熱點。為此，本研究基于機器學習的框架，提出了結合數據源自相關特征的空間自相關隨機森林水深反演模型，并在北島海域開展多光譜水深反演實驗。通過對比改進方法與對數比值、普通克里金、隨機森林方法的實驗結果，得到以下結論：

（1）Stumpf 模型通過建立波段反射率與水深之間的數學關系進行預測，具有模型構造簡單、對實測數據量要求較低等優(yōu)勢，因此得到了廣泛應用。然而在水深反演過程中，受到來自深度、水質以及底質等因素變化影響，反射率與水深之間的線性關系通常不成立，而Stumpf 這類線性模型受限于模型的表達能力，漸漸無法滿足日益增長的高精度遙感測深需求。對此，以RF 模型為代表的機器學習方法憑借其優(yōu)異的非線性映射能力，可有效提升水深反演精度。

（2）實驗表明，以OK 模型為代表的空間插值方法僅僅利用實測水深的自相關性特征信息，即取得與Stumpf 模型相似的水深預測精度。這一結果客觀上印證了本文通過引入空間自相關特征信息以提高水深反演精度這一思路的可行性。對此，本文在機器學習框架下引入空間自相關特征進行建模，提出了空間自相關隨機森林模型，通過深入挖掘不同數據源之間、相鄰像元內部觀測值的有效信息，將輸入變量的統(tǒng)計互相關性與空間自相關性兩類特征進行融合，得到了更準確的預測模型。該方法比單獨使用兩類特征的反演模型更具優(yōu)勢。SARF 模型在輸入反演因子時引入兩類相關性變量的信息，提高了預測過程中的信噪比，具有比傳統(tǒng)機器學習更出色的抗噪能力。同時，SARF 模型在建模過程中充分利用訓練樣本的觀測值，通過建立訓練樣本與測試樣本之間的聯系，從而提高有效信息的利用率。除此之外，將數據源的空間自相關特征融入機器學習算法中，有效減弱了空間非平穩(wěn)性因素、極淺區(qū)域高輻射亮度值以及深水區(qū)統(tǒng)計值波動的影響，可顯著提升不同實測數據量條件下的多光譜水深反演結果的精度。

（3）鑒于SARF 模型在遙感測深研究中取得了顯著的精度優(yōu)勢，基于相同思路，將自相關性特征引入其他機器學習模型中亦可望取得類似的提升效果。需要注意的是，本研究區(qū)域內，水深主要集中于0～6 m淺水區(qū)，而考慮到不同研究區(qū)數據源的差異，將其應用于更多地區(qū)還需要進一步探討?；诂F有方法，未來的研究方向可考慮挖掘額外的環(huán)境信息，利用機器學習融合多變量的優(yōu)勢，增加更多環(huán)境相關特征（如水體指數、葉綠素濃度）并融合不同時序的多光譜影像數據源進行研究。