南海珊瑚礁光學淺水海域多光譜遙感水深反演

汪萬智 ，鄒欣慶 , ，李海宇 , *

1. 南京大學 地理與海洋科學學院，南京 210023；2. 南京大學 海岸與海島開發(fā)教育部重點實驗室，南京 210023；3. 中國南海研究協(xié)同創(chuàng)新中心 南京大學，南京 210093

1 引言

南海珊瑚礁海域在全球海洋生物多樣性、生態(tài)系統(tǒng)與海洋經(jīng)濟可持續(xù)發(fā)展方面具有重要價值，是全球最具初級生產(chǎn)力的海洋生態(tài)環(huán)境之一（余克服，2012）。水深數(shù)據(jù)是海洋科學研究的重要基礎(chǔ)資料，對海洋環(huán)境研究、生態(tài)保護、開發(fā)利用和航運安全等具有重要意義。傳統(tǒng)海洋測深使用船載聲納探測方法，具有測量精度高，水深量程大的特點。聲納測深需實地走航實施，對于大范圍海域測量效率較低，測量成本高、周期長、易受施測區(qū)域和天氣海況條件限制。機載激光雷達（LiDAR）測深具有覆蓋面積廣，測量周期短、精度高的優(yōu)點，近年來得到迅速發(fā)展，對于可以低空進入的近岸島礁和淺水海域能實施有效探測。然而南海海域島礁眾多，水下礁盤面積大，遠離大陸，受自然條件和社會經(jīng)濟等多種因素限制，空中及海面實測均存在諸多困難，因此水深數(shù)據(jù)在時間與空間方面非常有限，亟待補充。衛(wèi)星遙感具有覆蓋范圍廣、成本低、周期短、不受區(qū)域限制等明顯優(yōu)勢。在遙感影像中包含微弱的海水深度信息，如何應(yīng)用遙感影像快速獲取可靠的水深數(shù)據(jù)在理論與實踐方面受到各界的廣泛關(guān)注（黃榮永等，2019）。

衛(wèi)星遙感水深反演由于應(yīng)用目的不同有多種方法，包括微波主動遙感和光學被動遙感水深探測等多種途徑（馬毅等，2018）。本項研究采用多光譜光學被動遙感的方法。目前光學水深反演方法主要包括基于物理模型的解析方法和基于統(tǒng)計模型的經(jīng)驗方法兩種主要類型（Kerr et al., 2018），包括查找表法、基于光譜匹配的優(yōu)化算法以及其它多種形式（Mobley et al., 2005; Hu et al., 2014, Hedley et al., 2018;王紀坤等, 2018; Liu et al., 2019; Chu et al., 2019）。物理解析模型基于光線在水中的輻射傳輸模型，考慮水分子、葉綠素、有色溶解有機物（CDOM）、無機懸浮物、以及水流擾動、底質(zhì)等對光的吸收和散射的影響，研究在沒有實測數(shù)據(jù)的情況下反演水深。該模型具有較強的通用性，但所需參數(shù)較多，求解較為困難，一般需要高光譜影像數(shù)據(jù)。統(tǒng)計經(jīng)驗模型主要根據(jù)水深與遙感影像光譜輻射值的統(tǒng)計關(guān)系建立對應(yīng)的統(tǒng)計模型，可以忽略其中的物理機制，相對便捷、容易實施。其中半分析半經(jīng)驗法從輻射傳輸模型出發(fā)，經(jīng)過一定數(shù)學變換減少所需的解析參數(shù)，建立經(jīng)驗模型，求解相對容易，如單波段線性模型（Lyzenga, 1978; Manessa et al., 2018）和對數(shù)比值模型（Stumpf et al., 2003; Ma et al., 2014; Kerr et al., 2018）。在單波段模型的基礎(chǔ)上還發(fā)展出區(qū)分底質(zhì)類型（聶榮娟，2019）、地理加權(quán)回歸（Chybicki et al., 2017）、地理自適應(yīng)（Su et al., 2014; 郭曉雷等，2017）、多波段線性回歸（陳本清等，2017）等多種方法。半分析半經(jīng)驗模型相對于物理解析模型所需參數(shù)大為減少，可以用于多光譜遙感影像，但需要實測數(shù)據(jù)確定模型參數(shù)。由于統(tǒng)計模型以特定區(qū)域的實測數(shù)據(jù)為基礎(chǔ)，不同區(qū)域、不同影像的回歸系數(shù)不同，限制了統(tǒng)計相關(guān)模型的應(yīng)用范圍。

南海珊瑚礁海域遠離大陸，實測數(shù)據(jù)較難取得。針對其島礁淺水海域海水光學特性以葉綠素為主導的特點，本研究以南沙群島道明群礁中的庫歸沙洲為實驗研究區(qū)域（圖1），采用Sentinel-2多光譜遙感影像，從物理解析模型出發(fā)，通過MODIS衛(wèi)星影像計算水體葉綠素濃度，使用半分析模型計算在該葉綠素濃度下不同水深的遙感反射率，采用對數(shù)比值模型進行遙感水深反演。并進一步利用多時相遙感影像進行反演水深融合，去除噪聲，提升反演水深精度。與多波束實測水深的誤差分析顯示，該方法對珊瑚礁光學淺水海域具有較好的水深反演效果。

2 數(shù)據(jù)和方法

2.1 研究區(qū)域

研究區(qū)域位于中國南海庫歸沙洲（10°46′40″N，114°35′40″E），包括兩個不相互聯(lián)系的礁坪，高潮時淹沒，屬于南海南沙群島道明群礁的一部分。該區(qū)域距最近的大陸超過600 km，受大陸和人類活動影響小。海水水質(zhì)清澈，能見度高，屬于典型海洋I類水體，適合遙感水深反演。

2.2 遙感影像數(shù)據(jù)與預處理

本項研究選取成像質(zhì)量良好的三期Sentinel-2A多光譜衛(wèi)星影像和對應(yīng)日期的MODIS影像進行水深提取實驗（表1）。Sentinel-2包括2A和2B兩顆基本相同的衛(wèi)星傳感器，互補重訪周期為五天，搭載的多光譜成像儀（MSI）包括從藍光（442 nm）到短波紅外（2202 nm）13個波段。針對水深反演，我們選取其中第2、3、4和第8波段，分別對應(yīng)藍、綠、紅和近紅外（NIR）波段，空間分辨率為10 m。MODIS影像用于提取海水葉綠素濃度。

表1 Sentinel-2A衛(wèi)星影像信息Table 1 Image Parameters of the satellite Sentinel-2A

我們首先對影像數(shù)據(jù)進行輻射定標、大氣校正、去云層和波浪破碎泡沫處理。其中大氣校正使用的是歐洲航天局官方提供的Sen2Cor，云層和波浪破碎泡沫利用水體在近紅外波段的強吸收特性，使用閾值分割加以區(qū)分。由于本研究區(qū)位于熱帶海域，相對于中高緯度地區(qū)，更容易受耀斑的影響（Kay et al., 2009）。耀斑校正能顯著提高水深反演的精度（Goodman et al., 2008；沈蔚等, 2019)，因此我們進一步利用NIR在水體穿透性極低，且菲涅爾反射對NIR和可見光基本相同的原理，對于研究區(qū)域海洋表面數(shù)據(jù)采用下列公式進行耀斑校正（Hochberg et al., 2003; Hedley et al., 2005）：

其中R和R′分別代表校正前后的遙感反射率，VIS和NIR分別表示可見光和近紅外波段，Rmin代表沒有耀斑時的遙感反射率，ρ是耀斑校正系數(shù)，取值為深水區(qū)VIS和NIR波段線性回歸系數(shù)（Hedley et al.,2005）。

2.3 葉綠素濃度

本項研究使用MODIS影像進行葉綠素濃度提取，研究假設(shè)在研究區(qū)內(nèi)葉綠素濃度空間分布基本一致，并使用反演的葉綠素濃度作為研究區(qū)的背景葉綠素濃度?？紤]到OC2和OC4等用于全球的反演模型（O’ Reilly, 1998）在局部地區(qū)可能誤差較大，本研究采用李新星等（2015）基于MODIS針對南海建立的葉綠素濃度Ca的反演算法，公式為：

2.4 水深測量與校正

本研究采用實測水深數(shù)據(jù)對遙感反演水深結(jié)果進行精度評價分析，以驗證方法的可行性。實測數(shù)據(jù)采用多波束回聲測深儀（Sonic, 2024, R2Sonic Inc., USA）于2017年8月采集自南海庫歸沙洲海域。樣本數(shù)據(jù)經(jīng)重采樣生成與衛(wèi)星影像相同空間分辨率的水深數(shù)據(jù)并校正至理論深度基面。由于實測船只現(xiàn)場數(shù)據(jù)采集條件限制，實測水深最淺處約8 m，絕大多數(shù)數(shù)據(jù)點分布在10～30 m的范圍，故誤差分析針對10～30 m的深度范圍進行。

使用衛(wèi)星影像反演水深得到的是衛(wèi)星成像時刻的水深，為了與多波束水深數(shù)據(jù)進行比較，需要進行潮汐校正，以消除潮汐的影響。本研究使用TOPEX/ Poseidon Global Inverse Solution（TPXO）模型校正潮位，具體方法參考Egbert等（1994,2002）。三個時期研究區(qū)域衛(wèi)星成像時刻的潮位分別為0.21 m、0.01 m和0.19 m（表1）。

2.5 輻射傳輸模型

水深反演半分析模型基于光線在水中輻射傳輸單次散射理論，將衛(wèi)星接收到的水體上行信號近似表達為由海底反射和海水后向散射兩部分組成。在光學淺水區(qū)域，假設(shè)海底反射為漫反射，對于均質(zhì)海水，水面以下的遙感反射率rrs可以表示為海底反射率ρb、總衰減系數(shù)k和水深z的函數(shù)（Lee et al., 1998）：

圖2 遙感水深反演流程圖Fig. 2 The flow chart of water depth inversion model

2.6 對數(shù)比值模型

Stumpf等（2003）提出的對數(shù)比值模型是半理論半經(jīng)驗法常用的模型，表示為：

其中n是常數(shù)，用于保證rrs取對數(shù)后仍為正數(shù)，本研究中n取值為10000，下標1和2表示兩個不同的波段，本研究中分別為藍光和綠光波段。

式（4）中的ln(nrrsl)/ln(nrrs2)通常被稱為水深反演因子，傳統(tǒng)的對數(shù)比值模型使用實測水深控制點與水深反演因子建立回歸模型求解參數(shù)m0、m1，進而對遙感影像進行水深的反演。本研究假設(shè)研究區(qū)葉綠素濃度一致，使用半分析模型分別計算不同波段和深度的rrs，以此建立回歸求解m0、m1，避免了水深控制點的需求。

2.7 遙感反射率計算

考慮光線在水體中向上和向下傳播時，衰減系數(shù)的差異和視角變化引起的路徑延長，式（3）可進一步表示為（Lee et al., 1998）：

式中ρb為海底反射率，θ0和θ1分別為水下太陽天頂角和衛(wèi)星天頂角，分別是來自水體和海底散射的路徑延長因子，k為輻射傳播衰減系數(shù)，包括水體的總吸收系數(shù)at和總后向散射系數(shù)bb，即k=at+bb。這里rrs、Rb、at和bb等變量的取值依賴于波長λ。

對于I類水體，水體的總吸收率和總后向散射率主要受純凈海水和葉綠素的影響。根據(jù)生物光學模型，衰減系數(shù)可以表示為：

其中，kw為光線在純凈海水中傳播的衰減系數(shù)，由海水的吸收系數(shù)aw(λ)和散射系數(shù)bw(λ)構(gòu)成，即kw(λ)=aw(λ)+γbw(λ)，這里γ為后向散射系數(shù)。χ(λ)表征由于葉綠素濃度造成的光線在海水中傳播時衰減系數(shù)的變化（Lee et al., 1998; Morel et al.,1977, 1991; Morel, 1988）。

通過定義u(λ)=bb(λ)/k(λ)，則可以通過水體固有光學特性近似表示為（Gordon et al.,1988, Lee et al., 1998）:

在已知葉綠素濃度和底部反射率的條件下，使用式（5）正向推演得到的rrs的光譜分辨率為1 nm，而多光譜數(shù)據(jù)的光譜分辨率較低，需要借助衛(wèi)星傳感器的光譜響應(yīng)函數(shù)，計算波段平均的采用如下公式：

其中Si為波長為λ時衛(wèi)星光譜響應(yīng)函數(shù)的值。

利用式（10）得到的是水氣界面以下的遙感反射率rrs，遙感影像耀斑校正后得到的是水氣界面以上的遙感反射率Rrs，依據(jù)Fresnel-Snell定理二者在水氣界面之間存在如下關(guān)系：

式中Γ≈1.56，ξ為水氣界面入射光線幾何關(guān)系的函數(shù)（Lee et al., 1998, 1999; Kerr et al., 2018）。

2.8 底部反射率

珊瑚礁生態(tài)系統(tǒng)中鈣質(zhì)沉積的反射率最高，較易區(qū)分，而海草、珊瑚和海藻等生物群落的光譜十分相似（Hochberg et al., 2003，徐京萍等，2019；Petit et al., 2017），水深越深，越難區(qū)分。 Garcia 等（2018）使用高光譜數(shù)據(jù)在10～12 m深度范圍僅能達到53.5%的分類精度，本研究使用多光譜數(shù)據(jù)對此不作細分。采用Lyzenga（1981）的方法，消除輻射傳輸模型的水深參數(shù)，得到表示反射率相對大小的灰度圖像（圖3b），進而使用閾值分割法將底質(zhì)劃分為類別Ⅰ和類別Ⅱ兩類（圖3c），類別Ⅰ的反射率顯著低于類別Ⅱ。同時采用楊君怡（2017）于西沙群島測得的底質(zhì)光譜數(shù)據(jù)集（圖4a），依據(jù)反射率大小將其劃分為“高反射率”和“低反射率”兩類（圖4b），“高反射率”類使用“鈣質(zhì)細沙”的光譜，剩余類別反射率的平均值作為“低反射率”類的光譜，并計算二者的平均光譜。分別使用“低反射率”類和“高反射率”類的光譜作為類別Ⅰ和類別Ⅱ的光譜。

圖3 底質(zhì)分類過程圖。圖中紅線范圍為實測區(qū)域Fig. 3 The classification of bottom types. The red line indicates the area with in situ water depth

圖4 (a)珊瑚礁不同底質(zhì)類型的平均光譜 (b)光譜混合后的平均光譜Fig. 4 (a) Reflectance spectra for different benthic substrates (b) Mixed spectra of classified bottom types

2.9 精度評價

使用均方根誤差（RMSE）和平均相對誤差（MRE）來評價水深反演的總體結(jié)果，計算方法分別為：

式中，Δzi代表數(shù)據(jù)點反演的水深值與實測水深值的差值，n為經(jīng)重采樣實測水深點數(shù)量。相同深度范圍內(nèi)RMSE值越小水深反演精度越高，MRE可用于比較不同深度范圍的相對誤差，值越小精度越高。

3 結(jié)果

3.1 葉綠素濃度與遙感反射率

使用2017年5月19日研究區(qū)Sentinel-2的太陽天頂角、衛(wèi)星天頂角和“高反射率”類的光譜曲線，假設(shè)葉綠素濃度分別為0、0.05、0.1、0.5和1 mg/m3，分別計算藍光波段和綠光波段在不同深度的rrs，結(jié)果見圖 5。

圖5 藍光波段和綠光波段rrs與深度對應(yīng)關(guān)系Fig. 5 The subsurface remote sensing reflectance rrs of blue and green bands depends on water depth for different Chl-a concentrations

對于I類水體，在不同葉綠素濃度下，隨著深度增加，藍光和綠光的遙感反射率都逐漸減小，衰減速率隨深度增加而減小，并最終趨近于0，此時遙感反射率的值即為葉綠素濃度越低，rrs隨水深衰減越慢。相同深度不同葉綠素濃度下，藍光的rrs變化范圍顯著大于綠光，說明藍光的遙感反射率受葉綠素濃度影響更大。

對數(shù)比值模型的水深反演因子與水深z在一定的深度范圍內(nèi)呈線性關(guān)系，葉綠素濃度越低，呈線性關(guān)系的深度范圍越廣，該深度范圍可認為是對數(shù)比值模型的有效深度。圖 6顯示葉綠素濃度為0時有效深度約35 m，葉綠素濃度為0.1 mg/m3時有效深度約為20 m，葉綠素濃度為1 mg/m3時有效深度為15 m左右。

圖6 不同葉綠素濃度下對數(shù)比值模型水深反演因子與深度對應(yīng)關(guān)系Fig. 6 The relationship of depth inversion factor of log-ratio model and water depth for different Chl-a concentrations

3.2 水深反演

利用通過MODIS遙感影像獲取的2017年5月19日研究區(qū)葉綠素濃度（0.06 mg/m3）和半分析模型，分別計算不同深度下藍光波段和綠光波段的rrs。在底質(zhì)分類的基礎(chǔ)上建立對數(shù)比值模型，對經(jīng)過耀斑校正和水氣界面反射校正后的影像進行水深反演。模型反演水深和實測水深的對應(yīng)關(guān)系見圖 7。在不同底質(zhì)和深度范圍的RMSE和MRE見表2。

表2 不同深度范圍水深反演精度評價Table 2 Accuracy evaluation of water depth inversion in different depth ranges

對數(shù)比值模型在不同深度范圍內(nèi)，類別Ⅰ和類別Ⅱ的RMSE與MRE差異不明顯。圖 7中可見對數(shù)據(jù)點均勻分布于y=x兩側(cè)，反演深度和實際深度對應(yīng)較好。雖然“低反射率”類的光譜與類別Ⅰ的光譜不匹配，但從水深反演結(jié)果來看，對數(shù)比值模型不同底質(zhì)的數(shù)據(jù)點分布并無顯著差別，說明對數(shù)比值模型對不同底質(zhì)類型的適應(yīng)性較好。在不使用實測水深數(shù)據(jù)調(diào)整模型參數(shù)的情況下，對于單幅影像在10～30 m深度范圍內(nèi)總體水深反演精度RMSE為2.88 m，MRE為10.7%。

圖7 反演水深和實測水深結(jié)果對比散點圖Fig. 7 Scattering plot of retrieved depth vs. in situ depth

4 討論

4.1 底質(zhì)類型影響

海底反射率在不同底質(zhì)類型下有顯著差異，是影響水深反演結(jié)果的重要因素（Manessa et al.,2018）。對數(shù)比值法反演水深能夠在一定程度上消除底質(zhì)的影響。本項研究將庫歸沙洲海域在影像上劃分為高反射率與低反射率兩種不同底質(zhì)類型。結(jié)合單波段水深反演模型，分析兩種底質(zhì)類型的遙感反射率rrs，顯示在不同的底部反射率光譜下，各模型的水深反演因子都與水深高度相關(guān)（R2≥0.98）。其中對數(shù)比值模型在約30 m水深范圍內(nèi)呈現(xiàn)良好的線性關(guān)系，而且對于高低兩種底質(zhì)類的線性擬合關(guān)系高度一致（圖 8）。需要指出的是不同底質(zhì)類型反映了不同的底部生物群落環(huán)境，受生物作用影響，在潟湖、珊瑚礁邊緣部分水深小于10 m范圍內(nèi)海水IOP受葉綠素與CDOM共同作用明顯，從而影響水深反演結(jié)果（Russell et al., 2019）。

圖8 不同底質(zhì)光譜下單波段模型和對數(shù)比值模型水深反演因子與深度對應(yīng)關(guān)系Fig. 8 The relationship between water depth inversion factors and water depth in the single-band model and ratio algorithm model under different bottom types

4.2 多時相反演水深平均

單幅影像水深反演存在一定的隨機誤差，同時在數(shù)據(jù)預處理過程中，建立掩膜去除影像中的云層和礁平臺附近淺水區(qū)的波浪破碎等噪聲信息，會導致水深反演結(jié)果存在少量像素缺失數(shù)據(jù)。因此，在假設(shè)海底地形不變的情況下，通過多時相影像序列反演水深，可以填補缺失，獲取完整的研究區(qū)水深反演數(shù)據(jù) (Chu et al., 2019; Barnes et al., 2018), 而且可以減小隨機誤差，在一定程度上提高反演精度。本項研究使用多時相反演水深算術(shù)平均的方法，用于填補數(shù)據(jù)缺失和減少隨機誤差，進而提高水深反演精度和效果。

通過對研究區(qū)域三期Sentinel-2遙感影像水深反演比較（表3），2018年3月15日和2018年4月4日影像水深反演結(jié)果在不同深度范圍內(nèi)誤差都較大，2017年5月19日影像水深反演誤差較小。但通過三期影像的平均融合，不同深度范圍內(nèi)的RMSE和MRE都有所減小，整體RMSE從2.88 m降至2.68 m，MRE從10.7%降至9.99%。從平均后的水深圖上（圖 9）可以看到云及陰影等多種噪音信號得到抑制和消除，顯著提高了水深反演效果。對于水深較淺的區(qū)域，反演水深與實際深度總體對應(yīng)較好，而深度較大的區(qū)域，反演水深總體偏低，與前文分析相一致，即對數(shù)比值模型在水深反演因子脫離線性增長階段（約水深25 m）, 隨著深度增加，反演深度會低于實際深度。

圖9 多時相水深反演結(jié)果Fig. 9 The retrieved water depth using multi-temporal images

表3 三期影像反演水深平均融合前后精度變化Table 3 The accuracy improvement after averaging water depth retrieved from three images

5 結(jié)論

本研究使用Sentinel-2多光譜遙感衛(wèi)星影像數(shù)據(jù)，結(jié)合MODIS衛(wèi)星數(shù)據(jù)，以中國南海南沙群島庫歸沙洲為例，通過半分析模型和對數(shù)比值水深反演方法，嘗試在不依靠實測水深資料的情況下，進行該地區(qū)光學淺水海域水深測量，并通過多波速實測數(shù)據(jù)進行驗證分析。該方法通過MODIS數(shù)據(jù)反演研究區(qū)域內(nèi)的背景葉綠素濃度，使用半分析模型計算遙感反射率，據(jù)此建立水深反演模型。研究分析了不同光譜波段和底質(zhì)分類條件下對反演水深的影響，表明該方法在本地區(qū)對不同底質(zhì)有較好的適應(yīng)性。使用三期遙感影像進行反演水深并進行均值處理消除噪音信號，在不同深度范圍精度都得到提高。與實測數(shù)據(jù)相較，整體誤差RMSE為2.68 m，MRE為9.99%。本方法可以快速獲取初步海洋水深數(shù)據(jù)，并應(yīng)用于南海其它相似島礁海域，為相關(guān)海洋領(lǐng)域的研究和應(yīng)用提供數(shù)據(jù)支持。

致謝：本項研究得到中國科學院學部咨詢評議項目“海洋絲路地質(zhì)與海洋環(huán)境戰(zhàn)略分析”支持。編輯和匿名審稿專家對論文提出了十分寶貴和建設(shè)性意見，謹致謝忱。