ICESat-2水深提取方法及其反演應用

胡琪鑫, 程亮, 楚森森, 程儉, 徐亞

1 南京大學地理與海洋科學學院, 南京 210023

2 江蘇省地理信息技術重點實驗室, 南京 210023

3 中國科學院地質與地球物理研究所, 中國科學院油氣資源研究重點實驗室, 北京 100029

4 中國科學院大學, 北京 100049

0 引言

淺海水下地形數據作為海洋環境的基礎地理信息數據,是海洋工程建設與海洋科學研究的重要支撐.在海洋工程建設方面,淺海水下地形廣泛服務于航行安全(Mavraeidopoulos et al.,2017)、海上救援(Pelot et al.,2015)和資源勘探(Brando et al.,2009)等領域.在海洋科學研究方面,水下地形數據不僅是水動力學和波浪建模(O′Hara Murray and Gallego,2017;Khanarmuei et al.,2020)的重要輸入參量和海岸帶地貌與沉積物搬運研究(Kanari et al.,2020)的關鍵參考數據,還可以輔助海底底質分類與珊瑚礁棲息地繪圖(Eugenio et al.,2015;Hedley et al.,2018)等.

當前常用的水深測量方法主要包括船載聲吶(單波束、多波束、側掃聲吶等)測深,機載LiDAR(Light detection and ranging)測深和基于衛星影像的水深反演.其中,船載聲吶測深精度高,測深能力強但費用昂貴,且在復雜地形或極淺水域難以開展工作;機載LiDAR測深精度和效率高,適用于極淺水域,但依舊存在費用昂貴、不適合大范圍測深的缺點(Kotilainen and Kaskela,2017);基于衛星影像的水深反演研究興起于1972年美國發射的第一顆陸地資源衛星Landsat-1,其優點是不受地域限制,覆蓋面積大,費用低且獲取方便;缺點是被動光學信號微弱,易受環境影響,精度較低,且大多數經驗統計反演模型仍然需要其他手段的測深數據作為訓練樣本.近年來,隨著搭載主動激光測高系統的衛星陸續升空以及對其水下測深能力的驗證(Neumann et al.,2019;Parrish et al.,2019),結合主動激光測深與被動光學影像的主被動融合水深反演成為大范圍、低成本獲取淺海水下地形的主流研究方向.

ICESat-2光子點水深提取主要包括海面、海底光子點群提取和光子點水深值改正三部分內容.目前,ICESat-2海面、海底光子點群提取方法大多借鑒于傳統的機載光子云數據處理算法,主要包括基于柵格圖像處理的提取方法,基于空間濾波的提取方法,基于分布統計特征的提取方法和基于密度空間聚類的提取方法.(1)基于柵格圖像處理的提取方法將光子點云轉換為柵格圖像然后對其作去噪或濾波處理,原理簡單,但在轉換過程中容易損失部分精度(Chen and Pang,2015);(2)基于空間濾波的提取方法如基于體素的空間濾波(Tang et al.,2016),主要適用高信噪比的光子點云數據,難以適用于高背景噪聲的ICESat-2數據;(3)基于分布統計特征的提取方法對光子點云的空間分布進行局部或全局的特征統計,利用統計特征差異區分信號點和光子點,如局部距離統計(夏少波等,2014)、窗格統計(Popescu et al.,2018)、直方圖統計(Hsu et al.,2021;Ranndal et al.,2021)等,基于分布統計特征的提取方法精度高,效果明顯,但仍然不足以應對海底光子點群信噪比低、分布場景多樣化的特點;(4)基于密度空間聚類的提取方法主要有貝葉斯、OPTICS(Ordering Points to Identify the Clustering Structure)、DBSCAN(Density-Based Spatial Clustering of Applications with Noise)三種,目前應用最為廣泛的是DBSCAN及其改進算法,DBSCAN算法對掃描半徑和最小包含點數兩個參數極其敏感,基于DBSCAN的改進算法主要解決參數尋優的問題,如基于距離矩陣的候選集尋優(Ma et al.,2020;Xie et al.,2021)、基于點云自身分布特性的K-平均最近鄰尋優(孟文君等,2021)、基于聚類分析的自適應尋優等(秦磊等,2020).目前,基于DBSCAN及其改進的提取方法在海底光子點群提取中具有較好的提取結果,但均未考慮光子點密度隨深度的變化特點,且在不同軌跡長度、海面起伏度的海面光子點提取場景中適用性較差.

提取海面和海底光子點后,需要對光子點水深值進行地球物理改正和水體折射改正.目前,應用最廣泛的水體折射改正模型是Parrish等(2019)基于Snell折射定律提出的水體折射改正幾何模型,其他模型均是基于Parrish模型的改進或簡化模型,如Ma 等(2020)在Parrish模型的基礎上考慮海面波浪造成的高程差異、Hsu等(2021)基于Parrish模型統計分析結果得到簡化經驗模型等.水體折射改正是光子點水深改正項中數量級最大的改正項,而現有的改正模型均未考慮海面傾斜造成的入射角偏移.

1 研究區與數據源

1.1 研究區概況

如圖1所示,本研究選取兩個典型的珊瑚礁作為研究區,分別為位于澳大利亞GBR(Great Barrier Reef)中部的John Brewer Reef 和南部的Fitzroy Reef.研究區水質清澈,葉綠素和懸浮物濃度低,含有豐富的珊瑚、沙石、藻類等底質.John Brewer Reef實驗區位于澳大利亞昆士蘭州Townsville市東北約38海里處,地理范圍為18.66°S—18.61°S、147.02°E—147.08°E,總面積約38.97 km2.John Brewer Reef是一個約6 km寬的巨大裸露中陸架礁石,由西南向東北方向延伸,水深呈現從東南到西北逐漸變淺的趨勢,最深處達≥18 m.受氣候變化和人類活動的影響,John Brewer Reef珊瑚白化現象嚴重,現已處于地質年齡的衰老階段.為此,英國雕塑家和海洋保護主義者Jason DeCaires Taylor在其東北角修筑了深達18 m的人工構筑物群——Museum of Underwater Art(MOUA),用以警示人類保護海洋環境.Fitzroy Reef實驗區位于澳大利亞昆士蘭州Seventeen Seventy鎮東北約32海里處,地理范圍為23.64°S—23.60°S、152.12°E—152.19°E,總面積約38.04 km2.Fitzroy Reef是Bunker Group中最大的干燥封閉環礁,擁有極其豐富的珊瑚、魚和海龜種類.Fitzroy Reef礁盤內存在大量水深低于6 m的區域,周圍礁脊至深海區域水深變化劇烈,中部為GBR南段唯一自然形成的全潮汐入口瀉湖,最深處達≥10 m.

1.2 數據源

1.2.1 ICESat-2 ATL03光子點

美國國家航空航天局(National Aeronautics and Space Administration, NASA)于2018年9月發射成功的冰、云與陸地高程衛星二號(Ice, Cloud, Land Elevation Satellite-2, ICESat-2)搭載了先進地形激光測高系統(Advanced Topographic Laser Altimeter System, ATLAS),ATLAS的主要構成是一臺激光波長為藍綠波段(532 nm)的光子計數激光雷達,可實現全球海拔高度監測.ATL03數據集,即全球定位光子點(Global Geolocated Photons),是ICESat-2數據集的一個子集,數據級別為Level 2.研究選取ATL03 2021年11月發布的最新Version 5產品,ATL03 V5更新了與光子點往返時間相關的測距誤差改正.本研究從EARTHDATA(https：∥search.earthdata.nasa.gov/search)獲取了兩個研究區自2018年10月13日至2022年6月30日的所有數據,為獲取經過研究區的所有軌跡條帶,設置范圍經緯度方向分別外擴0.2°,分別獲得96、109個項目文件,數據量分別為5.06、5.73 GB.

1.2.2 Sentinel-2多光譜影像

Sentinel系列衛星是歐盟委員會和歐洲航天局(European Space Agency, ESA)提出的全球環境與安全監測系統項目(“哥白尼計劃”)的專用衛星系列.其中,Sentinel-2用于陸地監測,可提供植被、土壤和水覆蓋、內陸水路及海岸區域等圖像,還可用于緊急救援服務.Sentinel-2現階段在軌運行衛星包括Sentinel-2A和Sentinel-2B兩顆同高度軌道、相位相差180°的衛星,分別于2015年6月23日和2017年3月07日發射成功,Sentinel-2C計劃于2024年發射.Sentinel-2的主要有效載荷為一枚多光譜成像儀(Multispectral Imaging ,MSI),在可見光/近紅外(VNIR)和短波紅外光譜范圍(SWIR)中具有13個光譜通道.Sentinel-2衛星提供的數據產品包括Level-1B、Level-1C、Level-2A,其中,Level-1B為原始數據,需要專業的校正技術;Level-1C是經過輻射校正和幾何校正的大氣上層表觀反射率產品.本研究基于GEE(Google Earth Engine：https：∥earthengine.google.com)獲取了兩個研究區自2018年1月1日至2022年6月30日所有云量像素百分比低于10%的Sentinel-2 L1C影像并進行空間裁剪和波段選擇,最終通過重新計算更為精確的云量分數并排序,兩個研究區最佳影像分別選定為John Brewer Reef： 20190822時相和Fitzroy Reef： 20180902時相,影像大小分別為2.70 MB、2.54 MB.

圖1 研究區與數據源(a)John Brewer Reef 和Fitzroy Reef 的地理位置;(b)、(c)John Brewer Reef、Fitzroy Reef的Sentinel-2影像(彩色背景)與ICESat-2軌跡(綠色條帶)((b)中黃色星號代表MOUA)

2 ICESat-2光子點水深提取方法

2.1 ICESat-2 ATL03數據預處理

針對ICESat-2數據量大但有效信息占比并不高的特點,對ICESat-2原始H5格式數據做文件重命名、有效信息存儲、空間裁剪、高程截取、無效文件去除等數據預處理.ICESat-2 ATLAS的理論最大測深深度約為40 m以淺,已有實例區域測深能力為約1個Secchi深度(Tyler,1968).因此,本文利用Isolation Forest異常點檢測算法(Liu et al.,2008)初步確定海面點平均高程,截取海面高程下50m至海面高程上5m的光子點.預處理后,兩個研究區分別保留了4、18條包含有效海底光子點的波束條帶,如圖1.如表1所示,每條波束條帶軌跡分別提取每個光子點的時空定位與信號水平信息(/gtx/heights： delta_time、lon_ph、lat_ph、h_ph、dist_ph_along、signal_conf_ph)、地球物理改正項信息(/gtx/geophys_corr： dac、tide_earth、tide_load、tide_ocean、tide_equilibrium、tide_oc_pole、tide_pole)和分段地理信息(/gtx/geolocation： segment_id、ref_azimuth、ref_elev)并以波束條帶為基本單元保存為CSV文件.

表1 ATL03 V5 字段含義Table 1 Field meanings of ATL03 V5 dataset

2.2 基于光子點群高斯分布特征的海面光子點群提取方法

(1)滑動窗口設置：ATL03 V5 波束條帶軌跡長度從幾千米到幾萬米不等,當處理窗口大小過小時,光子點數量較少可能造成高斯分布不明顯;當處理窗口大小過大時,海面高程可能出現較大起伏造成雙峰擬合偏離;經試驗,設置滑動窗口大小和步長均為10000 m時,可以較好地提取出代表海面光子點的高斯主峰.

(2)分段雙重高斯擬合：統計計算滑動窗口內的所有光子點高程分布直方圖,直方圖中高程區間寬度采用Freedman-Diaconis準則(Freedman and Diaconis,1981).Freedman-Diaconis準則基于最小化直方圖和理論概率分布密度之間平方差的積分和,可以根據數據分布確定合適的區間寬度,具備良好的抗噪性.對光子點高程分布直方圖進行雙重高斯擬合,雙重高斯模型概率密度函數如下：

(1)

其中,α1,α2為振幅,μ1,μ2為期望,σ1,σ2為標準差.設定當雙重高斯擬合標準差最大值≥0.95或期望值之差<0.3時,則判定雙重高斯分布不明顯,用單高斯模型代替進行重新擬合.高斯擬合后的主峰期望即為平均海面高程,選取高程范圍[μ-3σ,μ+3σ]的光子點作為該分段的海面候選光子點.滑動循環重復此步驟,即可得到波束條帶軌跡內的所有海面候選光子點,得到海面候選光子點集.

(3)尺度縮放：海面候選光子點軌跡方向點間距和高程方向點間距存在明顯不平衡的尺度差異,為利于下一步DBSCAN聚類,將海面候選光子點軌跡方向距離作1/3縮放處理.

(4)DBSCAN聚類去噪：與基于距離的聚類算法不同,DBSCAN作為基于密度空間聚類算法中的經典算法,具備區分出非球狀復雜形狀簇并找出不屬于任何簇的噪聲點的能力和優勢(Ester et al.,1996;Khan et al.,2014).本文采取滑動窗口分段處理的策略對海面候選光子點作DBSCAN去噪,設置滑動窗口大小和步長分別為100 m和50 m.DBSCAN聚類涉及鄰域的距離閾值eps和樣本點數量閾值min_samples兩個參數,本文利用sklearn.optimize模塊實現該算法,設置eps=1,min_samples=4.

2.3 基于深度自適應DBSCAN的海底光子點群提取方法

提取海面光子點群后,去除海面光子點候選集,在海面光子點候選集以深的地方設置0.2m的深度緩沖區,將緩沖區下的剩余光子點作為海底光子點群提取的輸入數據.深度自適應的DBSCAN海底光子點群提取算法的具體步驟如下：

(1)尺度縮放：同海面光子點分布類似,海底光子點空間分布軌跡方向和高程方向存在尺度差異,為與DBSCAN算法的聚類原理一致,將軌跡方向距離范圍縮放至與高程方向范圍一致,縮放比例通過如下公式計算：

(2)

其中,xmax,xmin分別為沿軌跡方向距離最大、最小值,hmax,hmin分別為高程最大、最小值,r為軌跡方向距離縮放因子,xr為縮放后的軌跡方向距離.

(2)深度自適應分段：引入hwin、hstep兩個參數分別控制深度分段窗口大小和步長,對每個深度段光子點的軌跡距離分布做Kolmogorov-Smirnov(KS test)(Massey,1951)均勻分布檢驗,KS test假設檢驗的置信度水平閾值設為：返回值statistic<0.1且pvalue>0.05時,則認定該深度段光子點軌跡距離分布符合均勻分布,該段深度段無海底光子點,對其不做處理.

(3)軌跡方向分段：對不符合均勻分布假設深度段的光子點做軌跡距離方向分段處理,引入xwin、xstep兩個參數分別控制軌跡方向距離分段窗口大小和步長.同(2)對每個分段做KS test 非參數檢驗,置信度水平閾值不變,對符合假設的分段不做任何處理.

(4)距離矩陣和DBSCAN候選eps集確定：對不符合均勻分布的分段,計算分段內光子點間軌跡方向的距離,形成一個N×N的距離矩陣：DN×N={dist(i,j)|1≤i≤N,1≤j≤N},其中N為分段內的光子點總數.對距離矩陣的每行分別升序排列,則每行的第1列元素均為0,第k列元素代表每個點的k近鄰距離Dk,所有行Dk的均值即組成候選的eps集epsN={epsk|1≤k≤N}.

(5)信號/噪聲分塊：將步驟(3)后的分段進一步在軌跡方向上分成M塊,則分段光子點的塊密度為

(3)

統計M塊中光子點數大于ρ的塊數M1以及M1塊的光子點總數N1,將這些塊標識為以信號光子點為主導的塊,則以噪聲光子點為主導的塊數和光子點總數為

M2=M-M1,

(4)

N2=N-N1.

(5)

(6)候選min_samples集確定：依次對epsN中的每個候選epsk按照下式計算分段內以epsk為半徑的圓內的總平均光子點點數Nall和噪聲為主導的平均光子點數：

(6)

(7)

其中,h′max,h′min為分段內的光子點高程最大、最小值,x′rmax,x′rmin為分段內光子點尺度縮放后軌跡距離的最大、最小值.最后,依據下式即可計算得到每個候選epsk對應的候選min_samples：

|1≤k≤N},

(8)

(7)迭代確定最優epsopt和min_samplesopt參數：從k=1開始按照升序依次將對應的epsk、min_samplesk代入DBSCAN算法對分段內光子點進行聚類,當聚類類別數連續三次不變時,最后一次DBSCAN的聚類參數即為最優參數epsopt和min_samplesopt.

(8)DBSCAN最優參數聚類：按照最優參數閾值epsopt和min_samplesopt對分段內光子點進行聚類,得到分段內海底光子點群.

(9)重復(2—8)步驟,即完成一次深度自適應DBSCAN算法聚類去噪.

(10)對每條波束軌跡條帶做2次大尺度深度自適應DBSCAN聚類和1次小尺度深度自適應DBSCAN聚類(次數依據實際光子點分布狀況調節),大尺度深度自適應DBSCAN聚類旨在去除遠離有效海底光子點的大量噪聲點,參數設置為：hwin=5,hstep=2.5,xwin=1,xstep=100000,M=100,而小尺度深度自適應DBSCAN聚類旨在去除有效海底光子點鄰近的少量噪聲點,提取更可靠的有效海底光子點,參數設置為hwin=2,hstep=1,xwin=100,xstep=50,M=50.

2.4 水深改正方法

2.4.1 地球物理改正

不同地球物理改正項作用的光子點群有所不同,且大部分改正已在ATL03中實現.因此,本文對光子點的地球物理改正實際包括(1)對海表光子點作動態大氣改正(Dynamic Atmospheric Correction, DAC)與逆氣壓(Inverted Barometer, IB)效應改正;(2)對折射改正后的初步水深值作海洋潮汐改正.改正公式如下：

Dgc=R(Hbot-(Hsur-ΔHDAC))+Δd+ΔHtide_ocean

+ΔHtide_equilibrium,

(9)

其中,Dgc為水體折射改正和地球物理改正后的水深值,R(·)為折射改正函數,Hbot為提取的海底光子點高程,Hsur為提取的海面光子點高程,Δd為區域平均海平面與瞬時海面光子點的高程差.ΔHDAC,ΔHtide_ocean,ΔHtide_equilibrium分別為動態大氣改正與逆氣壓效應改正項和海洋潮汐、均衡潮改正項.

2.4.2 水體折射改正

如圖2所示,當波束以入射高度角φ入射時,由于海面波浪等影響,折射界面偏離水平,導致入射角θ1≠φ,且當瞬時海面與軌跡方向之間的傾角符號不同時,實際光子點與觀測光子點之間的相對位置發生變化,最終導致折射改正計算公式不一.入射高度角φ可以從ATL03 V5數據集中的“ref_elev”字段計算得到,每個海底光子點對應海面傾角η利用軌跡方向鄰近的5個點擬合直線求取.公式如下：

圖2 折射改正幾何模型

(10)

k=polyfit(neighbor(5)),η=arctan(k).

(11)

根據圖2中的幾何關系和Snell折射定律,入射角θ1和折射角θ2可由下式計算：

θ1=Abs(φ-η),

(12)

(13)

其中,na,nw分別為大氣和海水折射率,na=1.00029,nw=1.34116.

同樣地,由折射定律可以分別求得入射點到實際光子點和觀測光子點間的距離R、S：

φ=θ1-θ2,

(14)

(15)

(16)

其中,φ為ΔRSP中邊R和S之間的夾角,D=Hbot-(Hsur-ΔHDAC)為根據海面與海底光子點高程得到的初步水深值.求出φ,R,S后,改正距離P和R對應的角α可由余弦定理求出：

(17)

(18)

針對不同的η值,改正距離P的傾角β不同：

(19)

隨即,每個海底光子點對應的水體折射改正項求解出：

ΔY=Pcosβ,

(20)

ΔZ=Psinβ,

(21)

其中,由于入射高度角非常小,在水深30 m時,ΔY≈9 cm,這相對于17 m波束足跡而言,可以忽略不計(Parrish et al.,2019),即只需加上水深方向的折射改正分量,得到折射改正后的光子點水深值Drc：

Drc=R(D)=D+ΔZ.

(22)

3 主被動融合水深反演方法

3.1 Sentinel-2影像預處理方法

(1)大氣校正

為消除或減少大氣分子、氣溶膠的散射和臭氧、水汽等氣體吸收對地物反射率的影響(Vermote et al.,2002),需要對下載的原始L1C級產品進行大氣校正,將大氣上層表觀反射率轉換為大氣下層表觀反射率.本文選用歐洲航天局官方最新發布的Sen2Cor 2.9插件對兩個研究區的Sentinel-2 L1C 影像作批量化大氣校正,Sen2Cor核心算法基于大氣輻射傳輸模型libRadtran,相比于大氣校正簡化模型和6S(Second Simulation of a Satellite Signal in the Solar Spectrum)模型,擁有更高的精度和可靠性(蘇偉等,2018).

(2)耀斑校正與空間域濾波

為削弱太陽耀斑噪聲,本文采取Hedley等(2005)提出的耀斑校正方法對耀斑現象嚴重的影像進行校正,公式如下：

R(λi)sg=R(λi)-ki[R(λNIR)-min(R(λNIR))],(23)

其中,R(λi),R(λi)sg分別為耀斑校正前后可見光波段λi(i=blue,green,red)的遙感反射率值,R(λNIR)為近紅外波段的遙感反射率值,ki為深水區R(λi)與R(λNIR)線性擬合的斜率,本文選擇兩個研究區影像的西北角100×100像元的區域作為深水區.為抑制遙感影像空間域噪聲,以3×3像元的平滑窗口對耀斑校正后的遙感影像進行均值濾波.

(3)樣本集生成與標準化

將提取的ICESat-2的水深點與處理后的遙感影像疊加,逐個遍歷遙感影像10 m空間分辨率的所有像元,取單個像元內所有水深中值及其藍、綠、紅波段的反射率值作為水深反演實驗的樣本集,得到兩個研究區水深樣本總數分別為835、4701個,并以3∶7的比例在水深范圍內隨機均勻地抽取訓練集和驗證集,最終分別得到訓練點、驗證點250、585個和1410,3291個,如表2所示.為使水深反演不受樣本集不同輸入變量數量級不一致的影響,實驗利用sklearn的StandardScalers對反射率值和水深值均進行標準差和取均值標準化處理.

表2 水深反演樣本集Table 2 Sample set of bathymetry inversion

3.2 水深反演模型

(1)單波段模型

單波段(Single Band, SB)模型以輻射傳輸方程為基礎,利用光輻射通量沿水深的變化呈指數衰減的規律,建立影像可見光波段與水深的相對關系(Polcyn and Sattinger,1969),即：

R(λi)=R∞(λi)+airbie-kifz,

(24)

其中,R(λi)為可見光波段λi(i=blue,green,red)的反射率值,R∞(λi)為λi的深水區反射率值,ai,rbi,ki分別為反映波段λi太陽輻射、水體底部反射率和水體衰減對反射率影響的系數,f,z分別為水體路徑長度和水深.由上式,可以得到水深關于反射率的關系：

z=a(ln[R(λi)-R∞(λi)])+b,

(26)

單波段模型由輻射傳輸方程簡化推導而來,回歸系數a,b具備鮮明的物理意義.

(2)比值模型

Stumpf等(2003)在單波段模型的基礎上,證明了利用不同波段之間的比值可以有效削弱水體和底質對水深反演的影響,并推導出水深反演的比值(Band Ratio, BR)表達式：

(27)

其中,R(λi),R(λj)為不同波段λi,λj的反射率值,一般取藍、綠波段;m1,m0為模型參數,n為保證對數運算而設定的固定常數,一般取1000或10000(本文實驗取1000).

(3)多項式模型

將單波段模型和比值模型擴展到多個波段,即得到Lyzenga多項式(Lyzenga Polynomial, LP)模型(Lyzenga,1978,1985)：

(28)

其中,a0,ai(i=1,2,…,N)為擬合系數.基于Lyzenga一次多項式推廣到更高階的多項式即得到K階多項式模型(滕惠忠等,2009;Paredes and Spero,1983)：

(29)

ak1,k2,…,kN(i=1,2,…,N)為擬合系數,Xi=ln[R(λi)-R∞(λi)](i=1,2,…,N),通常K取2或3.

(4)支持向量回歸模型

與Lyzenga多項式模型類似,Vojinovic等(2013)對比值算法作了非線性改進,利用支持向量機(Support Vector Machine)構建高維特征向量與水深之間的關系：

(30)

其中,f為非線性函數.對于非線性問題,支持向量回歸(Support Vector Regression, SVR)模型使用核函數將輸入變量隱式地轉換到線性空間,核函數包括線性核函數、多項式核函數、sigmoid核函數、徑向基核函數(Radial Basis Function, RBF)等,常用的徑向基核函數如下式：

(31)

其中,xi,xj為特征向量.

(5)多層感知器模型

ANN根據神經元是否具有記憶功能分為前饋神經網絡和后饋神經網絡,其中,前饋神經網絡也叫多層感知器(Multi-Layer Perceptron, MLP),它由輸入層、若干隱藏層和輸出層依次級聯組成(Ceyhun and Yal?n, 2010),可由下式表達：

(32)

其中,oi,j+1為第j+1層的第i個神經元節點,om,j為第j層的第m個神經元節點,m=1,2,…,n,n為第j層的神經元節點數,wm,i為oi,j+1與上一次每個神經元節點om,j的連接權重,bi為偏移量,f為每一層的激活函數.在水深反演中,MLP的輸入不需要對數運算處理,為可見光波段的反射率值.

(6)隨機森林回歸模型

隨機森林(Random Forest, RF)是通過集成學習的Bagging思想將多棵決策樹集成的一種算法,常用于多特征變量監督分類,同時還可以用來解決回歸問題(Breiman,2001).同多層感知器模型一樣,將隨機森林回歸模型運用至水深反演時,輸入變量無需任何處理.隨機森林回歸模型參數設置為：決策樹數量n_estimators=[50,100,150,200],分割準則選擇均方誤差和平均絕對誤差criterion=[‘mse’,‘mae’],最大深度max_depth=[None,3,5,7,9,11],最小分割樣本數min_samples_split=[2,4,6,8,10],最小葉子節點樣本數min_samples_leaf=[1,2,3,4,5],同時采用bootstrap采樣和袋外驗證.

3.3 精度評價指標

為定量評價水深反演精度,本文采用均方根誤差(Root Mean Squared Error, RMSE)、平均絕對誤差(Mean Absolute Error, MAE)、平均絕對百分比誤差(Mean Absolute Percentage Error, MAPE)、決定系數(coefficient of determination, r2_score)和Person相關系數(Pearson correlation coefficient)的平方R2共5個精度評價指標.各評價指標計算公式如下：

(33)

(34)

(35)

(36)

(37)

4 結果分析

4.1 ICESat-2水深提取

圖3為Fitzroy研究區某軌跡條帶上根據ATL03數據集提供的光子點置信度水平所作的空間分布和頻數統計圖,從中可以看到ICESat-2光子點分布具有如下幾個特征：(1)高程方向明顯的雙高斯分布和有效光子點密度隨深度減小特征;(2)海面光子點密度聚集的同時存在大量噪聲點,且海面局部傾斜;(3)光子點置信度水平能夠部分區分有效光子點與噪點,但依舊無法成為完整精確提取有效光子點群的有力憑據.而基于光子點高斯分布特征的海面光子點提取和基于深度自適應DBSCAN的海底光子點提取方法,能夠依次提取出有效的海面和海底光子點群,去除大量噪聲點.經過地球物理改正和折射改正后,光子點水深整體變淺,且越深的光子點水深值改正越大,如圖4所示.

通過對兩個研究區的所有有效光子點水體折射改正項和改正前水深進行線性回歸分析,可以發現改正項與改正前水深呈強線性關系,如圖5所示,回歸系數為-0.2550,這與Parrish等的折射改正模型計算結果-0.25416十分接近,二者之間的差別由改進模型考慮海面傾斜造成.

4.2 主被動融合水深反演

4.2.1 水深反演結果

基于提取并生成的ICESat-2光子點水深與Sentinel-2影像反射率樣本集,利用單波段(SB)、比值(BR)、Lyzenga多項式(LP)、二次多項式(Quadratic Polynomial, QP)、三次多項式(Cubic Polynomial, CP), 支持向量回歸(SVR), 多層感知器(MLP)和隨機森林回歸 (RF)等8種水深反演模型在兩個研究區分別開展主被動融合水深反演實驗.基于以上8種不同模型,制作水下地形反演DEM(Digital Elevation Model)如圖6、圖7所示,除了SB模型和BR模型,其他6種模型均能反演出較為精細的水下地形,而SB模型和BR模型由于模型缺陷,水深反演結果整體偏淺,難以恢復深水區域的水下地形,如John Brewer Reef的東北角區域和Fitzroy Reef的西南角深水區域.

圖3 Fitzroy_track0268_20191014_gt1l條帶光子點分布和頻數統計

圖4 Fitzroy_track0268_20191014_gt1l條帶光子點提取與改正后分布

圖5 折射改正-水深線性回歸圖

圖6 基于8種反演模型的John Brewer Reef水下地形圖(a) 單波段模型; (b) 比值模型; (c) Lyzenga多項式模型; (d) 二次多項式模型; (e) 三次多項式模型; (f) 支持向量回歸模型; (g) 多層感知器模型; (h) 隨機森林回歸模型.

圖7 基于8種反演模型的Fitzroy Reef水下地形圖(a) 單波段模型; (b) 比值模型; (c) Lyzenga多項式模型; (d) 二次多項式模型; (e) 三次多項式模型; (f) 支持向量回歸模型; (g) 多層感知器模型; (h) 隨機森林回歸模型.

圖8 John Brewer Reef實測水深與反演水深散點圖(a) 單波段模型; (b) 比值模型; (c) Lyzenga多項式模型; (d) 二次多項式模型; (e) 三次多項式模型; (f) 支持向量回歸模型; (g) 多層感知器模型; (h) 隨機森林回歸模型.

圖9 Fitzroy Reef實測水深與反演水深散點圖(a) 單波段模型; (b) 比值模型; (c) Lyzenga多項式模型; (d) 二次多項式模型; (e) 三次多項式模型; (f) 支持向量回歸模型; (g) 多層感知器模型; (h) 隨機森林回歸模型.

4.2.2 模型預測表現分析

對驗證集水深值和模型預測的水深值分布進行線性回歸并繪制散點圖,如圖8和圖9所示.綜合兩個研究區的水深反演結果發現：MLP、高階多項式、RF模型的水深反演效果最好,在淺水和深水水域都擁有很好的反演表現,反演水深可達近30 m,其中MLP的表現最穩定;LP和SVR模型能較好地反演出0～30 m的水深,但在0～5 m的極淺水域和20～30 m的更深水域反演效果欠佳,這主要是受到模型復雜度不夠或基于比值的基本原理所限制,復雜的非線性項可以更好地替代比值處理對底質影響的削弱作用;BR模型反演效果不穩定,在水深分布范圍較窄的淺水區域反演結果尚可,但不具備深水反演能力,在樣本點數更多水深范圍更廣的區域如Fitzroy Reef出現明顯的深水區域預測更淺現象;SB模型反演表現最差,對水深與反射率關系的擬合能力最差,在實際的水深反演和水下地形產品生產中可用性較弱.

4.2.3 模型定量精度分析

進一步對兩個研究區的訓練集和驗證集分別計算RMSE、MAE、MAPE、r2_score,定量評價不同反演模型的整體精度.如表3和表4所示,在John Brewer Reef研究區,RF模型訓練集的擬合效果最好,訓練集的RMSE=0.34 m,MAE=0.21 m,MAPE=3.19%,r2_score=0.99,但測試集的預測精度并非最高,這說明RF模型存在一定程度的過擬合;MLP、RF、多項式模型的整體反演精度最高,訓練集和測試集精度均達到了RMSE<0.7 m,MAPE<0.5 m,MAPE<7%,r2_score≥0.97,其中高階多項式模型比LP一次多項式模型精度更高;SVR、BR模型的反演精度次之,訓練集和測試集精度均達到了RMSE<0.9 m,MAPE<0.65 m,MAPE<13.5%,r2_score≥0.95,且SVR模型精度比BR模型更高;SB模型的精度最差,訓練集和驗證集的RMSE>2 m,MAPE>1.5 m,MAPE>38%,r2_score<0.7.

表3 John Brewer Reef水深反演整體精度評價Table 3 Overall accuracy evaluation of bathymetry inversion in John Brewer Reef

表4 Fitzroy Reef水深反演整體精度評價Table 4 Overall accuracy evaluation of bathymetry inversion in Fitzroy Reef

由于水深范圍和數據量不同,Fitzroy Reef研究區的整體反演精度比John Brewer Reef研究區略差.在Fitzroy Reef研究區,與John Brewer Reef研究區一致,RF模型擬合效果最好,訓練集RMSE=0.55 m,MAE=0.26 m,MAPE=7.53%,r2_score=0.99,但存在一定程度的過擬合;MLP、高階多項式、RF的測試集精度最高,其中MLP和CP模型的測試集精度達到了RMSE=0.79 m,MAE=0.39 m,r2_score=0.98;SVR、LP模型精度次之,訓練集和測試集精度RMSE<1.2 m,MAE 0.75 m,MAPE<31%,r2_score≥0.96,與John Brewer Reef研究區不同的是,SVR模型的精度高于LP模型,原因是當數據量增多,水深范圍變大時,雖然基于BR模型,非線性的SVR模型相比于線性的低階多項式模型具有更好的反演能力.

綜合兩個研究區的水深反演精度：MLP、高階多項式、RF的水深反演精度最高,其中,RF模型擬合效果最好,但容易出現一定程度的過擬合現象;SVR模型、LP模型精度次之,且當數據量增多,水深范圍變大時,非線性的SVR模型相比于線性的低階多項式模型具有更好的反演能力;BR模型精度較低,SB模型的精度最差,均不適合用于實際水深反演.

4.2.4 模型反演深度分析

與John Brewer Reef不同,在Fitzroy Reef,MLP和RF模型反演水深最大值同多項式模型一樣,均達到了23 m以上,見表5.對比分析MLP和RF模型在兩個研究區的不同表現,這是John Brewer Reef研究區的水深樣本最大深度只有13.907 m導致的,同時說明MLP和RF模型存在一定程度的過擬合,高度依賴水深樣本分布.雖然3種多項式模型在深水區域均有不錯的反演表現,當多項式階次過高時,多項式模型的反演結果嚴重發散,極其不穩定,如在John Brewer Reef,QP和CP的反演水深最大值分別達到了27.658 m和45.818 m,與真實的水下地形嚴重不符.

表5 8種模型反演的水深范圍Table 5 Water depth range of bathymetric inversion based on eight models

5 結論與建議

(1)基于光子點群的高斯分布特征進行海面光子點提取,可以有效處理具有多種多樣光子點分布的波束軌跡條帶,對不同軌跡長度、不同海面起伏度的情況均有較好的適應性.海底光子點群密度隨深度增加而減少的先驗信息有助于海底光子點群提取,在DBSCAN聚類算法基礎上引入深度自適應和滑動窗口可以有效提高海底光子點群的提取效率和準確率.

(2)以Parrish折射改正幾何模型為基礎,考慮海面傾斜造成的入射角偏移,可以進一步改正水體折射效應對水深測量造成的影響.折射改正后光子點水深呈整體變淺,且水深越深,變淺越多的線性規律,改正項約為改正前水深的0.2550倍.

(3)結合ICESat-2數據和Sentinel-2影像在水質清澈的淺海域進行主被動融合水深反演效果較好.相比其他水深實測數據,ICESat-2測深數據具有高度一致性,這直接導致利用ICESat-2水深值進行水深反演的精度更高,可靠性更好.在眾多水深反演模型中,SB模型和BR模型反演精度較差,多項式模型和MLP、RF等機器學習模型反演效果相比于BR模型和SVR模型,反演效果更好,但MLP和RF模型的反演效果對水深樣本分布具有高度依賴性,且RF模型往往出現過擬合.當區域ICESat-2水深樣本豐富、深度范圍寬廣時,MLP模型是較為理想的可靠反演模型,而當區域ICESat-2水深深水樣本不足或先驗信息匱乏時,LP模型是相對穩妥的選擇.

致謝感謝審稿專家提出的修改意見和編輯部的大力支持!