ICESat-2水深提取方法及其反演應(yīng)用

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

1 南京大學(xué)地理與海洋科學(xué)學(xué)院, 南京 210023

2 江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室, 南京 210023

3 中國(guó)科學(xué)院地質(zhì)與地球物理研究所, 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室, 北京 100029

4 中國(guó)科學(xué)院大學(xué), 北京 100049

0 引言

淺海水下地形數(shù)據(jù)作為海洋環(huán)境的基礎(chǔ)地理信息數(shù)據(jù),是海洋工程建設(shè)與海洋科學(xué)研究的重要支撐.在海洋工程建設(shè)方面,淺海水下地形廣泛服務(wù)于航行安全(Mavraeidopoulos et al.,2017)、海上救援(Pelot et al.,2015)和資源勘探(Brando et al.,2009)等領(lǐng)域.在海洋科學(xué)研究方面,水下地形數(shù)據(jù)不僅是水動(dòng)力學(xué)和波浪建模(O′Hara Murray and Gallego,2017;Khanarmuei et al.,2020)的重要輸入?yún)⒘亢秃０稁У孛才c沉積物搬運(yùn)研究(Kanari et al.,2020)的關(guān)鍵參考數(shù)據(jù),還可以輔助海底底質(zhì)分類與珊瑚礁棲息地繪圖(Eugenio et al.,2015;Hedley et al.,2018)等.

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

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

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

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)源

1.2.1 ICESat-2 ATL03光子點(diǎn)

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

1.2.2 Sentinel-2多光譜影像

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

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

2 ICESat-2光子點(diǎn)水深提取方法

2.1 ICESat-2 ATL03數(shù)據(jù)預(yù)處理

針對(duì)ICESat-2數(shù)據(jù)量大但有效信息占比并不高的特點(diǎn),對(duì)ICESat-2原始H5格式數(shù)據(jù)做文件重命名、有效信息存儲(chǔ)、空間裁剪、高程截取、無效文件去除等數(shù)據(jù)預(yù)處理.ICESat-2 ATLAS的理論最大測(cè)深深度約為40 m以淺,已有實(shí)例區(qū)域測(cè)深能力為約1個(gè)Secchi深度(Tyler,1968).因此,本文利用Isolation Forest異常點(diǎn)檢測(cè)算法(Liu et al.,2008)初步確定海面點(diǎn)平均高程,截取海面高程下50m至海面高程上5m的光子點(diǎn).預(yù)處理后,兩個(gè)研究區(qū)分別保留了4、18條包含有效海底光子點(diǎn)的波束條帶,如圖1.如表1所示,每條波束條帶軌跡分別提取每個(gè)光子點(diǎn)的時(shí)空定位與信號(hào)水平信息(/gtx/heights： delta_time、lon_ph、lat_ph、h_ph、dist_ph_along、signal_conf_ph)、地球物理改正項(xiàng)信息(/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 基于光子點(diǎn)群高斯分布特征的海面光子點(diǎn)群提取方法

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

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

(1)

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

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

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

2.3 基于深度自適應(yīng)DBSCAN的海底光子點(diǎn)群提取方法

提取海面光子點(diǎn)群后,去除海面光子點(diǎn)候選集,在海面光子點(diǎn)候選集以深的地方設(shè)置0.2m的深度緩沖區(qū),將緩沖區(qū)下的剩余光子點(diǎn)作為海底光子點(diǎn)群提取的輸入數(shù)據(jù).深度自適應(yīng)的DBSCAN海底光子點(diǎn)群提取算法的具體步驟如下：

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

(2)

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

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

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

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

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

(3)

統(tǒng)計(jì)M塊中光子點(diǎn)數(shù)大于ρ的塊數(shù)M1以及M1塊的光子點(diǎn)總數(shù)N1,將這些塊標(biāo)識(shí)為以信號(hào)光子點(diǎn)為主導(dǎo)的塊,則以噪聲光子點(diǎn)為主導(dǎo)的塊數(shù)和光子點(diǎn)總數(shù)為

M2=M-M1,

(4)

N2=N-N1.

(5)

(6)候選min_samples集確定：依次對(duì)epsN中的每個(gè)候選epsk按照下式計(jì)算分段內(nèi)以epsk為半徑的圓內(nèi)的總平均光子點(diǎn)點(diǎn)數(shù)Nall和噪聲為主導(dǎo)的平均光子點(diǎn)數(shù)：

(6)

(7)

其中,h′max,h′min為分段內(nèi)的光子點(diǎn)高程最大、最小值,x′rmax,x′rmin為分段內(nèi)光子點(diǎn)尺度縮放后軌跡距離的最大、最小值.最后,依據(jù)下式即可計(jì)算得到每個(gè)候選epsk對(duì)應(yīng)的候選min_samples：

|1≤k≤N},

(8)

(7)迭代確定最優(yōu)epsopt和min_samplesopt參數(shù)：從k=1開始按照升序依次將對(duì)應(yīng)的epsk、min_samplesk代入DBSCAN算法對(duì)分段內(nèi)光子點(diǎn)進(jìn)行聚類,當(dāng)聚類類別數(shù)連續(xù)三次不變時(shí),最后一次DBSCAN的聚類參數(shù)即為最優(yōu)參數(shù)epsopt和min_samplesopt.

(8)DBSCAN最優(yōu)參數(shù)聚類：按照最優(yōu)參數(shù)閾值epsopt和min_samplesopt對(duì)分段內(nèi)光子點(diǎn)進(jìn)行聚類,得到分段內(nèi)海底光子點(diǎn)群.

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

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

2.4 水深改正方法

2.4.1 地球物理改正

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

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

+ΔHtide_equilibrium,

(9)

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

2.4.2 水體折射改正

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

圖2 折射改正幾何模型

(10)

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

(11)

根據(jù)圖2中的幾何關(guān)系和Snell折射定律,入射角θ1和折射角θ2可由下式計(jì)算：

θ1=Abs(φ-η),

(12)

(13)

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

同樣地,由折射定律可以分別求得入射點(diǎn)到實(shí)際光子點(diǎn)和觀測(cè)光子點(diǎn)間的距離R、S：

φ=θ1-θ2,

(14)

(15)

(16)

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

(17)

(18)

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

(19)

隨即,每個(gè)海底光子點(diǎn)對(duì)應(yīng)的水體折射改正項(xiàng)求解出：

ΔY=Pcosβ,

(20)

ΔZ=Psinβ,

(21)

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

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

(22)

3 主被動(dòng)融合水深反演方法

3.1 Sentinel-2影像預(yù)處理方法

(1)大氣校正

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

(2)耀斑校正與空間域?yàn)V波

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

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為深水區(qū)R(λi)與R(λNIR)線性擬合的斜率,本文選擇兩個(gè)研究區(qū)影像的西北角100×100像元的區(qū)域作為深水區(qū).為抑制遙感影像空間域噪聲,以3×3像元的平滑窗口對(duì)耀斑校正后的遙感影像進(jìn)行均值濾波.

(3)樣本集生成與標(biāo)準(zhǔn)化

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

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

3.2 水深反演模型

(1)單波段模型

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

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

(24)

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

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

(26)

單波段模型由輻射傳輸方程簡(jiǎn)化推導(dǎo)而來,回歸系數(shù)a,b具備鮮明的物理意義.

(2)比值模型

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

(27)

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

(3)多項(xiàng)式模型

將單波段模型和比值模型擴(kuò)展到多個(gè)波段,即得到Lyzenga多項(xiàng)式(Lyzenga Polynomial, LP)模型(Lyzenga,1978,1985)：

(28)

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

(29)

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

(4)支持向量回歸模型

與Lyzenga多項(xiàng)式模型類似,Vojinovic等(2013)對(duì)比值算法作了非線性改進(jìn),利用支持向量機(jī)(Support Vector Machine)構(gòu)建高維特征向量與水深之間的關(guān)系：

(30)

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

(31)

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

(5)多層感知器模型

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

(32)

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

(6)隨機(jī)森林回歸模型

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

3.3 精度評(píng)價(jià)指標(biāo)

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

(33)

(34)

(35)

(36)

(37)

4 結(jié)果分析

4.1 ICESat-2水深提取

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

通過對(duì)兩個(gè)研究區(qū)的所有有效光子點(diǎn)水體折射改正項(xiàng)和改正前水深進(jìn)行線性回歸分析,可以發(fā)現(xiàn)改正項(xiàng)與改正前水深呈強(qiáng)線性關(guān)系,如圖5所示,回歸系數(shù)為-0.2550,這與Parrish等的折射改正模型計(jì)算結(jié)果-0.25416十分接近,二者之間的差別由改進(jìn)模型考慮海面傾斜造成.

4.2 主被動(dòng)融合水深反演

4.2.1 水深反演結(jié)果

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

圖3 Fitzroy_track0268_20191014_gt1l條帶光子點(diǎn)分布和頻數(shù)統(tǒng)計(jì)

圖4 Fitzroy_track0268_20191014_gt1l條帶光子點(diǎn)提取與改正后分布

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

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

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

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

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

4.2.2 模型預(yù)測(cè)表現(xiàn)分析

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

4.2.3 模型定量精度分析

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

表3 John Brewer Reef水深反演整體精度評(píng)價(jià)Table 3 Overall accuracy evaluation of bathymetry inversion in John Brewer Reef

表4 Fitzroy Reef水深反演整體精度評(píng)價(jià)Table 4 Overall accuracy evaluation of bathymetry inversion in Fitzroy Reef

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

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

4.2.4 模型反演深度分析

與John Brewer Reef不同,在Fitzroy Reef,MLP和RF模型反演水深最大值同多項(xiàng)式模型一樣,均達(dá)到了23 m以上,見表5.對(duì)比分析MLP和RF模型在兩個(gè)研究區(qū)的不同表現(xiàn),這是John Brewer Reef研究區(qū)的水深樣本最大深度只有13.907 m導(dǎo)致的,同時(shí)說明MLP和RF模型存在一定程度的過擬合,高度依賴水深樣本分布.雖然3種多項(xiàng)式模型在深水區(qū)域均有不錯(cuò)的反演表現(xiàn),當(dāng)多項(xiàng)式階次過高時(shí),多項(xiàng)式模型的反演結(jié)果嚴(yán)重發(fā)散,極其不穩(wěn)定,如在John Brewer Reef,QP和CP的反演水深最大值分別達(dá)到了27.658 m和45.818 m,與真實(shí)的水下地形嚴(yán)重不符.

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

5 結(jié)論與建議

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

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

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

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