習(xí)曉環(huán)，王子家，3，王 成，3

（1. 中國科學(xué)院空天信息創(chuàng)新研究院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室，北京 100094；2. 可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心，北京 100094；3. 中國科學(xué)院大學(xué)資源與環(huán)境學(xué)院，北京 100049）

近海岸水深測量是海洋測繪的基礎(chǔ)工作，傳統(tǒng)的借助于船載平臺的單波束或多波束測深儀操作簡單、直接，但人力物力消耗大，而且在淺水區(qū)域容易擱淺，測量效率低，難以滿足長時(shí)間序列的大面積動(dòng)態(tài)變化監(jiān)測需要。光學(xué)遙感數(shù)據(jù)源豐富、覆蓋范圍廣，在水深反演中發(fā)揮了一定的作用，但受成像技術(shù)等多種綜合因素影響，測量精度較低［1-2］。機(jī)載測深LiDAR（Light Detection And Ranging）可以兼顧精度與效率，實(shí)現(xiàn)對島礁周邊的水深探測，但數(shù)據(jù)獲取成本高，且常受空域管控等限制。星載LiDAR覆蓋范圍廣、重復(fù)頻率高、受外界因素影響小，在湖泊、冰川和海洋遙感研究中得到廣泛應(yīng)用［3-4］。

ICESat-2 （Ice， Cloud and land Elevation Satellite-2）是美國NASA 于2018 年9月發(fā)射的新一代星載LiDAR 衛(wèi)星，搭載的ATLAS（Advanced Topographic Laser Altimeter System）傳感器采用了多波束、微脈沖、光子計(jì)數(shù)技術(shù)［5］，具有更高的脈沖重復(fù)頻率，可以獲取小光斑、高密度的光子點(diǎn)云數(shù)據(jù)（光斑直徑約17m、同軌光斑間隔僅0.7m）。同時(shí)，該系統(tǒng)還采用了532nm 的單脈沖激光器，水體在該波段散射弱、衰減系數(shù)小，激光可穿透1倍的塞克板深度（即透明度），其最大測深能力可達(dá)38m［6］。由于ATLAS 系統(tǒng)發(fā)射脈沖是弱信號，受太陽背景、系統(tǒng)自身、大氣散射、地表覆蓋等因素的影響，接收的回波中包含大量噪聲信號，特別是白天獲取的數(shù)據(jù)更為嚴(yán)重。同時(shí)受水中懸浮物散射與吸收、風(fēng)浪等影響，水體垂直方向不同區(qū)域的噪聲分布也存在一定差異。Ma 等［7］利用機(jī)載模擬光子數(shù)據(jù)（Multiple Altimeter Beam Experimental Lidar，MABEL）信號光子和噪聲光子的分布特征，基于改進(jìn)的DBSCAN（Density-Based Spatial Clustering of Applications with Noise）方法提取地面和水下的表面輪廓；樓樂樂［8］提出一種基于光子點(diǎn)云密度統(tǒng)計(jì)并結(jié)合差分回歸高斯自適應(yīng)的去噪算法，在淺水區(qū)域取得了較高的精度，但對復(fù)雜水底地形的適用性有待進(jìn)一步探索；Chen 等［9］針對光子在水面和水柱的分布提出了一種自適應(yīng)橢圓濾波測深方法（Adaptive Variable Ellipse Filtering Bathymetric Method，AVEBM），需要依據(jù)一定的規(guī)則確定濾波閾值，且深水處易出現(xiàn)提取錯(cuò)誤，濾波結(jié)果仍受橢圓濾波器大小和分塊大小的參數(shù)影響。

由于水面光子點(diǎn)分布密集而水下光子噪聲多，且隨水深變化，信號光子密度在沿軌方向和高程方向也發(fā)生變化，如何有效剔除背景噪聲是ICESat-2/ATLAS 數(shù)據(jù)獲取近岸水深數(shù)據(jù)的關(guān)鍵。但目前光子點(diǎn)云去噪的研究多針對陸地類型，水域方面的研究相對較少。本研究利用南海島礁淺海區(qū)域的ATLAS光子數(shù)據(jù)開展近海岸水深提取方法研究，重點(diǎn)解決復(fù)雜水底地形及水體較深位置處水深信息的提取難題，并利用機(jī)載LiDAR測深數(shù)據(jù)驗(yàn)證本文方法的有效性和精度。

1 研究方法

ICESat-2/ATLAS 光子數(shù)據(jù)存在大量噪聲，在水深信息提取應(yīng)用時(shí)首先要進(jìn)行去噪處理，然后將去噪后的光子數(shù)據(jù)分為水面光子和水底光子；基于兩者的密度分布差異，分別采用區(qū)間估計(jì)和改進(jìn)OPTICS 算法進(jìn)行水體信號光子提取，并結(jié)合RANSAC 算法獲取水面高程信息。由于光子在水體傳輸中會(huì)產(chǎn)生折射以及海面受潮力產(chǎn)生周期性波動(dòng)，為了提高水深反演精度，需進(jìn)行折射校正和潮汐校正。最后，對實(shí)驗(yàn)結(jié)果進(jìn)行精度評價(jià)，包括利用人工標(biāo)注水體信號光子，對本文精去噪結(jié)果的定性和定量評價(jià)，并與ATL03 內(nèi)置去噪算法、AVEBM 算法結(jié)果進(jìn)行對比，以及利用機(jī)載測深LiDAR數(shù)據(jù)對提取的水深信息進(jìn)行的精度評定?？傮w技術(shù)流程如圖1所示。

圖1 基于ICESat-2/ATLAS光子數(shù)據(jù)的水深提取流程Fig.1 Bathymetry workflow based on ICESat-2/ATLAS data

1.1 光子粗去噪

通過兩步完成光子去噪，即粗去噪和精去噪。粗去噪即去除原始光子數(shù)據(jù)中比較明顯的噪聲，利用光子點(diǎn)云文件中的置信度參數(shù)“signal_conf_ph”去除參數(shù)為0 的光子，但去噪后仍然存在大量噪聲光子。

1.2 光子精去噪

考慮到水面和水底數(shù)據(jù)的分布密度差異較大，大部分信號光子集中在水面，因此通過構(gòu)建高程分布直方圖，利用混合高斯分布模型對其進(jìn)行兩次高斯函數(shù)擬合，以兩個(gè)高斯曲線的交點(diǎn)來確定水面光子的范圍，實(shí)現(xiàn)水面和水底光子的分離，其中二者精去噪的方法有所不同。

1.2.1 水面光子精去噪

由于水面光子分布較為集中，可利用區(qū)間估計(jì)法去除水底光子和其他噪點(diǎn)。統(tǒng)計(jì)水面光子的高程分位數(shù)，將［0，0.02］、［0.98，1］對應(yīng)的光子點(diǎn)作為噪點(diǎn)去除，從而得到水面光子數(shù)據(jù)，然后利用RANSAC（random sample consensus）直線擬合確定水面高程，計(jì)算水面光子的最低高程點(diǎn)并將其設(shè)為閾值，低于該高程點(diǎn)的光子認(rèn)為是水底光子。

1.2.2 水底光子精去噪

由于水底光子分布密度不均，隨水深增加，信號光子密度減小。通過改變?yōu)V波參數(shù)進(jìn)行兩次精去噪，最大程度提取水底信息。具體步驟如下：

（1）基于改進(jìn)密度聚類算法（OPTICS）對水底光子精去噪［10］，采用試錯(cuò)法，最終選定橢圓搜索半徑a和b分別為11和1，利用式（1）計(jì)算輸入點(diǎn)的閾值參數(shù)MinPts［8］：

式中：S1是橢圓搜索區(qū)域內(nèi)噪聲光子和信號光子的總數(shù)；N1是給定樣本的總光子數(shù)；h和l分別是樣本光子點(diǎn)云中垂直距離和沿軌距離的最大差值；S2是橢圓搜索區(qū)域內(nèi)包括的噪聲光子數(shù)。假設(shè)最低水深處向上5m的范圍為噪聲光子，N2為該范圍內(nèi)的光子數(shù)，則h2為5m。一般將沿軌方向10 000光子點(diǎn)定為樣本點(diǎn)，實(shí)際應(yīng)用中水下光子總數(shù)往往小于10 000，故通常是將所有光子數(shù)據(jù)參與計(jì)算。

（2）進(jìn)一步區(qū)分信號和噪聲光子。OPTICS 算法引入了兩個(gè)特征，即核心距離（Core Distance，CD）和可達(dá)距離（Reachability Distance，RD）。Cd表示核心距離，是對于輸入的MinPts 使該點(diǎn)p成為核心點(diǎn)的最小距離，而當(dāng)不滿足輸入?yún)?shù)時(shí)Cd(p)無意義。Rd表示可達(dá)距離，針對核心對象，若p為核心對象，Rd（p，q）表示樣本點(diǎn)q與樣本點(diǎn)p之間核心距離，即dist(p，q)的最大值；當(dāng)p不為核心對象時(shí)，Rd(p)無意義。

式中：ε是鄰域半徑；Nε(p)是以p為圓心；ε為半徑的鄰域內(nèi)的點(diǎn)數(shù)。

Rd值越大表明該光子離其他光子越遠(yuǎn)，即為噪聲光子，反之為信號光子。本文采用最大類間方差法（OTSU）確定信號和噪聲光子，通過對光子數(shù)據(jù)進(jìn)行分塊并統(tǒng)計(jì)每個(gè)分塊中光子數(shù)N。假設(shè)前t個(gè)光子為信號光子，剩下N-t個(gè)光子為噪聲；當(dāng)類間方差g最大時(shí)，將此時(shí)的Rd設(shè)為閾值Rdth；當(dāng)光子Rd＞Rdth時(shí)，此光子記為噪聲光子，反之記為信號光子。

（3）利用改進(jìn)的OPTICS算法進(jìn)行第一次濾波，然后對濾波結(jié)果進(jìn)行縱向分塊，構(gòu)建高程統(tǒng)計(jì)直方圖。假設(shè)最低水深處向上5m的范圍為噪聲光子，并以此為最低高度確定二次濾波范圍。以直方圖峰值為中心，向兩邊搜索高程點(diǎn)數(shù)小于平均高程點(diǎn)數(shù)的分塊，以一次濾波確定的最小鄰域光子數(shù)為基礎(chǔ)，擴(kuò)大橢圓搜索半徑進(jìn)行二次濾波。

1.3 水深信息提取

ICESat-2/ATLAS 數(shù)據(jù)產(chǎn)品（ATL03）僅考慮了激光在空氣介質(zhì)中的傳播。而實(shí)際應(yīng)用中，激光脈沖從空氣進(jìn)入水體會(huì)發(fā)生折射，導(dǎo)致光子位置偏移，可根據(jù)斯涅爾定律對水底光子點(diǎn)進(jìn)行折射校正，獲取光子點(diǎn)在沿軌方向和垂直方向的偏移量，然后基于水底信號光子計(jì)算各光子點(diǎn)位置的水深數(shù)據(jù)。折射校正時(shí)，需要根據(jù)激光單位指向矢量方位角將光子坐標(biāo)的經(jīng)緯度轉(zhuǎn)換到其對應(yīng)的位置。此外，由于海水受到潮力影響會(huì)產(chǎn)生周期性波動(dòng)，且光子數(shù)據(jù)提取的水深值是獲取時(shí)刻的瞬時(shí)水深，為了得到基于深度基準(zhǔn)面的水深值，本文利用T_TIDE 模型［11］選取8 分潮（M2、S2、K2、N2、K1、O1、P1、Q1）的潮汐調(diào)和常數(shù)對水深數(shù)據(jù)進(jìn)行潮汐校正。

1.4 精度評價(jià)

采用定性分析和定量指標(biāo)對精去噪結(jié)果進(jìn)行評價(jià)。定性評價(jià)主要基于目視識別，即通過噪聲和信號光子點(diǎn)的分布情況，人為判斷是否存在錯(cuò)分現(xiàn)象；定量指標(biāo)是利用人工標(biāo)注的水體信號光子，用召回率（Recall，R）、準(zhǔn)確率（Precision，P）和F值來評價(jià)。其中召回率表示被正確提取的有效信號光子個(gè)數(shù)占原有信號光子總數(shù)的比例，準(zhǔn)確率是指被正確提取的有效信號光子個(gè)數(shù)與提取到的有效信號光子總數(shù)的比值，F(xiàn)值是召回率和正確率的調(diào)和平均值。

最后，為了驗(yàn)證光子數(shù)據(jù)提取的水深精度，將機(jī)載測深LiDAR獲取的水深數(shù)據(jù)作為真值進(jìn)行比較。實(shí)際中由于光子位置與機(jī)載點(diǎn)云位置無法完全重合，本文將距離光子位置最近且距離值小于2m的機(jī)載測深LiDAR點(diǎn)云水深值作為真值，通過繪制散點(diǎn)圖直觀地顯示驗(yàn)證結(jié)果。

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)區(qū)與數(shù)據(jù)

選取2019—2021 年南海永樂環(huán)礁的ICESat-2/ATLAS數(shù)據(jù)（共19個(gè)條帶）進(jìn)行實(shí)驗(yàn)，包括羚羊礁、甘泉礁、晉卿島、甘泉島、珊瑚島等島礁。該區(qū)域水質(zhì)清澈含沙量低，退潮時(shí)會(huì)露出大面積的礁石，臺風(fēng)大潮時(shí)會(huì)被海水淹沒。研究區(qū)的ATL03 光子數(shù)據(jù)分布如圖2所示，每一條ATL03數(shù)據(jù)包括6個(gè)條帶，分 別 為GT1L、GT1R、GT2L、GT2R、GT3L 和GT3R，并記錄了每個(gè)光子的經(jīng)緯度、高程、光子往返時(shí)間、激光器位置和姿態(tài)角等信息，并對大多數(shù)誤差進(jìn)行了糾正（如固體潮汐、大氣延遲等），可獲取光子高度的最佳估計(jì)［12］。

圖2 實(shí)驗(yàn)區(qū)ICESat-2/ATLAS光子數(shù)據(jù)分布情況Fig.2 ICESat-2/ATLAS data distribution in the study areas

機(jī)載LiDAR 測深數(shù)據(jù)由中國科學(xué)院上海光學(xué)精密機(jī)械研究所于2018年8月28日獲取，覆蓋范圍為甘泉島和羚羊礁區(qū)域，所用設(shè)備為機(jī)載雙頻雷達(dá)激光測深系統(tǒng)Mapper5 000，最大和最小測深分別為51m 和0.25m，測深精度0.23m［13］。由于機(jī)載數(shù)據(jù)與ATLAS 數(shù)據(jù)的坐標(biāo)系不同，在精度評價(jià)前對其進(jìn)行坐標(biāo)系統(tǒng)轉(zhuǎn)換。

2.2 三種去噪方法結(jié)果分析

ICESat-2/ATLAS光子點(diǎn)云噪聲點(diǎn)分布范圍較廣，通過粗去噪可以去除明顯的噪聲點(diǎn)（圖3）。其中20190524GT3L數(shù)據(jù)水底地勢平坦，20200419GT3R數(shù)據(jù)水底地形復(fù)雜。圖4分別顯示了人工標(biāo)注水體信號光子、ATL03 內(nèi)置精去噪（置信度參數(shù)“signal_conf_ph=4”光子數(shù)據(jù)）、AVEBM 和本文方法精去噪結(jié)果。其中，圖4 c～f 中的圓圈表示ATL03 去噪和AVEBM 方法對光子錯(cuò)提和漏提位置?？梢钥闯?，針對平緩的水底地形，兩種方法均可有效去除噪聲點(diǎn)，但ATL03內(nèi)置去噪算法難以提取較深區(qū)域的信號光子，且容易在水底地形起伏較大的區(qū)域出現(xiàn)信號光子缺失；AVEBM 對較深處的噪聲光子易產(chǎn)生錯(cuò)提取。

圖3 光子粗去噪結(jié)果Fig.3 Coarse noise photons removal results

圖4 光子精去噪結(jié)果比較Fig.4 Comparison of fine noise photons removal results

進(jìn)一步進(jìn)行定量分析（表1），本文方法提取的水體信號光子可獲得更高的F值，平均為97.98%，而ATL03 內(nèi)置精去噪算法為92.55%，AVEBM 為94.78%，分別提高了約5.87%和3.38%。

表1 三種方法精去噪統(tǒng)計(jì)指標(biāo)結(jié)果Tab.1 Statistical indicators of fine noise photons removal results based on the three methods

2.3 水深信息精度評價(jià)

利用機(jī)載LiDAR 測深數(shù)據(jù)對本文光子數(shù)據(jù)提取的水深信息進(jìn)行精度評價(jià)。通過獲取的水面高程和折射校正后水底信號光子高程值，計(jì)算水底信號光子對應(yīng)位置的水深值。以20190524GT3L數(shù)據(jù)為例，對獲取的水面信號光子和經(jīng)過折射校正前后的水底信號光子采用B 樣條曲線對其進(jìn)行擬合，得到連續(xù)的水面曲線和水底地形（圖5）。

圖5 20190524GT3L數(shù)據(jù)折射校正前后結(jié)果圖Fig.5 Signal photons before and after refraction correction of 20190524GT3L data

由于淺海區(qū)域的細(xì)沙、碎沙較多，且機(jī)載測深受底質(zhì)物以及透明度的影響，分析中對較淺區(qū)域（水深值≤0.5m）的數(shù)據(jù)進(jìn)行了剔除。本文以羚羊礁區(qū)域20190524GT3L 和20200820GT3L 數(shù) 據(jù) 提 取 的 水 深結(jié)果為例進(jìn)行精度評價(jià)（圖6）。結(jié)果顯示，最大探測深度7.83m，最大誤差約1.69m，最小約0.001m，差值較大，可能因機(jī)載數(shù)據(jù)和ICESat-2/ATLAS 數(shù)據(jù)位置無法完全重合所導(dǎo)致；決定系數(shù)R2為0.91，均方根誤差RMSE為0.53m。

圖6 ATLAS測深結(jié)果與機(jī)載LiDAR測量驗(yàn)證Fig.6 Accuracy assessment of the bathymetric results obtained by ATLAS and airborne LiDAR

3 結(jié)論

ICESat-2/ATLAS光子計(jì)數(shù)數(shù)據(jù)應(yīng)用的關(guān)鍵是去除原始數(shù)據(jù)中大量的噪聲光子，本文以水體清澈、受雜質(zhì)影響較小的南海為實(shí)驗(yàn)區(qū)，將改進(jìn)OPTICS算法用于水底光子信息精去噪，降低了濾波參數(shù)的敏感性。通過對一次水底濾波結(jié)果分塊、擴(kuò)大搜索半徑等，快速、高精度提取水體深處更密集的信號光子，提高了復(fù)雜水底地形區(qū)域的水體光子提取精度。由于光子受水體及其他懸浮物散射的影響，獲取的水面數(shù)據(jù)存在一定的高程差異，本文在水面精去噪后重新計(jì)算了水底數(shù)據(jù)的范圍，雖一定程度上降低對淺水區(qū)域水底信息錯(cuò)分，但受高斯擬合曲線的影響，在水面波動(dòng)較大的區(qū)域仍存在部分水面光子錯(cuò)分為水底光子的現(xiàn)象；同時(shí)二次精去噪的濾波參數(shù)不能實(shí)現(xiàn)自動(dòng)化設(shè)置，還需進(jìn)一步提高方法的自動(dòng)化程度。此外，考慮到大部分近海區(qū)域水質(zhì)并非如南海島礁水質(zhì)好、透明度高，本文方法的適應(yīng)性還有待進(jìn)一步研究，這也是目前激光雷達(dá)測深應(yīng)用的難點(diǎn)。

作者貢獻(xiàn)聲明：

習(xí)曉環(huán)：學(xué)術(shù)指導(dǎo)，提出選題，設(shè)計(jì)論文框架，論文審閱與修改。

王子家：整理文獻(xiàn)，方法實(shí)現(xiàn)，數(shù)據(jù)整理，數(shù)據(jù)處理及精度驗(yàn)證，論文撰寫與修改。

王成：學(xué)術(shù)指導(dǎo)，論文審閱與修改。