基于無(wú)人機(jī)與衛(wèi)星圖像的大型綠藻孔石莼(Ulva pertusa)遙感估算研究

孟苗苗, 鄭向陽(yáng), 邢前國(guó), 劉海龍

1. 中國(guó)科學(xué)院海岸帶環(huán)境過(guò)程與生態(tài)修復(fù)重點(diǎn)實(shí)驗(yàn)室(煙臺(tái)海岸帶研究所), 山東 煙臺(tái) 264003;

2. 中國(guó)科學(xué)院海洋大科學(xué)研究中心, 山東 青島 266071;

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

近二十年來(lái), 由大型藻類暴發(fā)而導(dǎo)致的海洋災(zāi)害在世界范圍內(nèi)呈增長(zhǎng)趨勢(shì)(Morand et al, 2005; Ye et al, 2011; Smetacek et al, 2013), 不僅危害海洋生態(tài)環(huán)境(Valiela et al, 1997), 還影響海上交通、沿海地區(qū)的觀光旅游等活動(dòng)(邢前國(guó) 等, 2011)。在中國(guó),大型藻類災(zāi)害受到重視始自2008 年夏季, 時(shí)值奧運(yùn)會(huì)帆船比賽前, 青島沿岸暴發(fā)了以漂浮大型綠藻滸苔(U. prolifera)為主的綠潮災(zāi)害。國(guó)外, 類似的海洋大型藻類災(zāi)害也在增加, 如法國(guó) Brittany 的綠藻(Schreyers et al, 2021)、墨西哥灣與大西洋的漂浮馬尾藻帶(Wang et al, 2019; Gower et al, 2019)。

在我國(guó)近海與海岸帶地區(qū), 黃海夏季由紫菜養(yǎng)殖導(dǎo)致的大規(guī)模漂浮滸苔(Liu et al, 2009; Xing et al,2019)以及冬春季黃海、東海的漂浮馬尾藻銅藻(Xing et al, 2017; Qi et al, 2017)因規(guī)模較大受到了較多關(guān)注; 以底棲為主的大型藻也常在本地局部地區(qū)形成災(zāi)害或環(huán)境污染(Xing et al, 2016a; Han et al,2019), 由于在沿海分布廣泛, 其生態(tài)環(huán)境效應(yīng)與資源管理問(wèn)題不容忽視。

孔石莼(U. pertusa)是一種潮間帶與淺海常見(jiàn)的綠色大型藻, 其出現(xiàn)與水體的高度富營(yíng)養(yǎng)化密切相關(guān)?？资辉跍\海以底棲為主, 能在非附著的條件下生長(zhǎng), 附著氣泡較多時(shí)亦可在水面呈漂浮狀?？资簧L(zhǎng)迅速, 生物量大, 常因暴發(fā)、過(guò)度聚集而形成生態(tài)災(zāi)害或污染?？资灰灿兄娴纳鷳B(tài)功能與資源價(jià)值, 在生長(zhǎng)過(guò)程中可吸收大量營(yíng)養(yǎng)鹽凈化水體、吸收二氧化碳, 也可作為食材、動(dòng)物飼料, 或者深加工作為藥物、保健品等原材料(劉巖 等,2010)。

使用衛(wèi)星影像能夠大范圍、同步監(jiān)測(cè)綠藻在海面上的分布狀態(tài), 為綠潮的防控提供先決條件。微波影像可以根據(jù)粗糙度信息識(shí)別出漂浮在海面的綠藻(蔣興偉 等, 2009; Shen et al, 2014), 而熱紅外影像可以根據(jù)不同的熱發(fā)射率將綠藻從海水中識(shí)別出來(lái)(Xing et al, 2015)。監(jiān)測(cè)綠藻最常用的是光學(xué)影像, 因?yàn)榫G藻與陸上植被具有相似的光譜特性, 可見(jiàn)光波段反射率低, 紅光波段和藍(lán)光波段出現(xiàn)吸收谷, 近紅外波段出現(xiàn)反射峰(樊彥國(guó) 等, 2015), 具有典型的“紅邊效應(yīng)”(安德玉 等, 2018)。早期識(shí)別綠藻的遙感算法有差值植被指數(shù)(difference vegetation index, DVI)、 比值植被指數(shù)(ratio vegetation index, RVI) 、 歸一化植被指數(shù)(normalized difference vegetation index,NDVI)(Cui et al, 2012; Xing et al, 2018)， 其中NDVI 在光學(xué)影像綠藻提取中應(yīng)用最為廣泛(Li et al, 2018b)。由于植被指數(shù)具有一定局限性, 又有學(xué)者提出適用于綠藻提取的指數(shù), 例如漂浮藻類指數(shù)(floating algae index, FAI)(Hu, 2009; Wang et al,2016)和虛擬基線漂浮藻類指數(shù)(virtual-baseline floating macroalgae height,VB-FAH)(Xing et al,2016a)， 這些指數(shù)對(duì)太陽(yáng)光、氣溶膠變化等的影響敏感性較低, 更適用于海面綠藻的提取。

由于海面綠藻斑塊的大小可以從幾厘米到幾公里不等(Li et al, 2018a), 傳統(tǒng)的閾值提取方法屬于“硬分類”, 由于混合像元的影響, 在估算綠藻面積時(shí)存在一定誤差。針對(duì)混合像元分解的綠潮面積精細(xì)化提取, 有學(xué)者基于同步的高分辨率和低空間分辨率衛(wèi)星影像間分析了混合像元分解(丁一 等,2015)， 還有學(xué)者通過(guò)實(shí)測(cè)綠藻光譜, 分析綠藻提取指數(shù)對(duì)綠藻表面覆蓋度的響應(yīng), 并建立混合像元綠藻覆蓋度估算模型(Qi et al, 2016; Li et al, 2018b)。

本文針對(duì)孔石莼展開(kāi)研究, 利用實(shí)測(cè)同步厘米級(jí)分辨率無(wú)人機(jī)數(shù)據(jù)提取的綠藻面積, 建立基于Landsat-8 數(shù)據(jù)計(jì)算的植被指數(shù)(NGRDI、NDVI和 VB-FAH)以及單波段反射率對(duì)綠藻亞像素覆蓋度的多個(gè)反演模型, 并對(duì)模型進(jìn)行了驗(yàn)證與評(píng)價(jià)。

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

1.1 數(shù)據(jù)及數(shù)據(jù)預(yù)處理

1.1.1 衛(wèi)星與無(wú)人機(jī)數(shù)據(jù)采集與處理

本文研究區(qū)域位于渤海灣沿岸, 渤海為陸表海,平均水深18m (Xing et al, 2016b), 研究區(qū)域現(xiàn)場(chǎng)照片如圖1a, 整個(gè)研究區(qū)海面分布有大量孔石莼, 孔石莼在海里淺水區(qū)呈現(xiàn)沉水非定生狀態(tài), 退潮時(shí)呈干出狀態(tài)(圖1b、c)。本文采用準(zhǔn)同步的無(wú)人機(jī)影像和Landsat TM 衛(wèi)星影像開(kāi)展研究, 無(wú)人機(jī)影像采集時(shí)間為2020 年10 月23 日11 時(shí)8 分, Landsat 衛(wèi)星過(guò)境時(shí)間為10 時(shí)49 分。

無(wú)人機(jī)正射圖像采集使用大疆御Mavic 2 為搭載平臺(tái), 其搭載的傳感器為1/2.3 英寸CMOS, 有R、G、B 三個(gè)波段。使用Pix4Dmapper 軟件對(duì)自帶GPS定位信息的無(wú)人機(jī)影像進(jìn)行快速校正、拼接處理,對(duì)于無(wú)法拼接的無(wú)特征點(diǎn)圖像進(jìn)行手動(dòng)拼接處理,輸出滿足遙感解譯需求的影像產(chǎn)品, 拼接后的無(wú)人機(jī)影像如圖3c。對(duì)Landsat-8 OLI 數(shù)據(jù)首先進(jìn)行輻射定標(biāo)處理, 再使用FLAASH 大氣校正模型以獲取反射率數(shù)據(jù), 最后進(jìn)行裁剪處理。為保證地理配準(zhǔn)的精度, 將處理后的無(wú)人機(jī)影像與Landsat-8 影像進(jìn)行比較, 選取 15 個(gè)控制點(diǎn), 發(fā)現(xiàn)配準(zhǔn)結(jié)果較好,RMSE 為0.51。

1.1.2 高光譜數(shù)據(jù)采集與處理

為了解孔石莼的光譜特征, 幫助選擇參與建模的指數(shù), 于2020 年11 月11 日測(cè)量了研究區(qū)域不同狀態(tài)下孔石莼的反射率。采用的光譜儀為 ASD HandHeld2, 波長(zhǎng)范圍為325～1075nm, 光譜分辨率為1nm; 先測(cè)參考板的輻射亮度值, 再測(cè)海水與不同狀態(tài)下孔石莼的輻射亮度值, 每組重復(fù)測(cè)量3 次取平均值, 反射率計(jì)算公式為

式中R( )λ代表被測(cè)物體的反射率,L( )λ代表地物輻射亮度值,LP( )λ代表參考板輻射亮度值,ρ( )λ是參考板的反射率。

所測(cè)得的海水與孔石莼的反射光譜曲線如圖2所示。綠藻的光譜特征為在可見(jiàn)光波段反射率低,在近紅外波段反射率高, 與海水光譜的差異明顯。將綠藻光譜按照Landsat-8 OLI 傳感器第二到四波段的光譜響應(yīng)函數(shù)進(jìn)行重采樣處理, 發(fā)現(xiàn)其在藍(lán)、綠、紅、近紅外波段的反射率低于海水, 在近紅外波段反射率高于海水。

1.2 方法

1.2.1 無(wú)人機(jī)影像綠藻提取

在無(wú)人機(jī)影像上選取一塊大小為417×435 像素的樣本區(qū), 使用動(dòng)態(tài)閾值法對(duì)樣本區(qū)進(jìn)行綠藻提取,以人工目視解譯結(jié)果為基準(zhǔn), 對(duì)比提取結(jié)果, 進(jìn)行精度驗(yàn)證。

由于無(wú)人機(jī)影像只有紅、綠、藍(lán)三個(gè)波段, 因此選取4 種典型植被指數(shù): 歸一化綠紅差異指數(shù)(normalized green-red difference index, NGRDI)、歸一化綠藍(lán)差異指數(shù)(normalized green-blue difference index,NGBDI)、紅綠比值指數(shù)(red green ratio index,RGRI)、過(guò)綠指數(shù)(excess green, ExG)做綠藻提取研究, 其計(jì)算公式如下:

式中R、G、B分別表示紅波段、綠波段和藍(lán)波段像素值。

1.2.2 衛(wèi)星影像綠藻覆蓋面積與植被指數(shù)關(guān)系建模

(1) 植被指數(shù)構(gòu)建

將無(wú)人機(jī)影像上綠藻的提取結(jié)果以30m×30m的矩形為單位建立網(wǎng)格, 使之對(duì)應(yīng) Landsat-8 衛(wèi)星影像的每個(gè)像素。統(tǒng)計(jì)每個(gè)像素內(nèi)綠藻覆蓋面積除以該像素的總面積(900m2), 作為每個(gè)混合像元的大型綠藻覆蓋度(portion of macroalgae, POM)。根據(jù)研究區(qū)特點(diǎn)選取POM接近1 的像素作為典型綠藻像素,POM 接近0 的像素作為典型海水像素, 端元所在位置如圖3a 彩色標(biāo)注的像素(綠色表示典型的綠藻像素, 紅色表示典型的海水像素)。分析其光譜曲線(圖4)發(fā)現(xiàn), 在可見(jiàn)光波段, 海水的反射率值高于綠藻, 在近紅外波段, 海水反射率值降低, 而綠藻反射率值增高。有明顯區(qū)分特征的波段是藍(lán)、綠、紅、近紅外波段。

建立Landsat 衛(wèi)星影像混合像元POM 與光譜值之間的線性關(guān)系:y=ax+b,y為POM,x選取7 種指數(shù)參與建模分析, 分別為藍(lán)波段、綠波段、紅波段、近紅外波段的反射率與NGRDI、NDVI、VB-FAH，其計(jì)算公式如下:

式中R為反射率, 下標(biāo)NIR、RED、GREEN 分別表示近紅外、紅、綠波段。

(2) 模型建立與驗(yàn)證

在衛(wèi)星圖像上選取兩塊同樣大小的15×12 個(gè)像素組成的區(qū)域作為樣本區(qū)與驗(yàn)證區(qū), 如圖3a, 為避免配準(zhǔn)時(shí)的精度偏差, 以3×3 個(gè)像素為窗口, 計(jì)算每個(gè)窗口內(nèi)的平均POM 與每個(gè)窗口內(nèi)包含的9 個(gè)像素的平均光譜反射率,在樣本區(qū)域建立 POM 與NGRDI、NDVI、VB-FAH、藍(lán)波段反射率、綠波段反射率、紅波段反射率、近紅波段反射率的回歸模型, 選取有較好線性關(guān)系的模型公式在驗(yàn)證區(qū)域進(jìn)行驗(yàn)證, 對(duì)反演模型擬合的POM 結(jié)果與無(wú)人機(jī)影像提取結(jié)果擬合的POM 進(jìn)行對(duì)比分析, 選取最好的反演模型應(yīng)用到整片無(wú)人機(jī)覆蓋區(qū)域。

2 結(jié)果與討論

2.1 基于無(wú)人機(jī)影像的綠藻提取結(jié)果

在無(wú)人機(jī)影像上均勻選取一定數(shù)量的綠藻與海水樣本點(diǎn), 其地表像元亮度值(digital number, DN)均值曲線如圖5a 所示。發(fā)現(xiàn)綠藻在紅波段的DN 值明顯低于海水, 因此使用動(dòng)態(tài)閾值法提取綠藻時(shí)選用紅波段DN 值與4 種植被指數(shù)。為了確定基于此研究區(qū)域的無(wú)人機(jī)影像提取綠藻的最佳方法, 選取一塊樣本區(qū), 以目視解譯的結(jié)果作為綠藻提取的真實(shí)結(jié)果, 對(duì)比分析5 種指數(shù)的動(dòng)態(tài)閾值法提取的結(jié)果。

不同提取方法的提取結(jié)果如圖5b 所示, 從圖中可以看出, 使用ExG 指數(shù)提取的結(jié)果噪點(diǎn)最多。從表1 可以看出, 監(jiān)測(cè)精度最高的是紅波段DN 值, 提取一致性達(dá)到94%， Kappa 為0.88, F1 score 為0.92。且紅波段灰度圖(圖3d)可以反映出無(wú)人機(jī)影像單個(gè)像素的POM, 為便于之后進(jìn)行更精細(xì)的提取, 使用紅波段DN 值對(duì)無(wú)人機(jī)影像進(jìn)行綠藻提取, 最終提取結(jié)果如圖3e 所示。

表1 不同綠藻提取方法精度對(duì)比Tab. 1 Comparation on accuracy of green macroalgae extraction with different methods

2.2 衛(wèi)星影像建模結(jié)果討論

樣本區(qū)域的建模結(jié)果表明, POM 與可見(jiàn)光三個(gè)波段的反射率線性關(guān)系最好, 其R2都大于0.9, 與近紅外波段的線性關(guān)系最差。3 個(gè)指數(shù)中, POM 與NDVI、VB-FAH 的線性關(guān)系較好, 與NGRDI 的關(guān)系較差(圖6)。

選用可見(jiàn)光三個(gè)波段所建立的模型在驗(yàn)證區(qū)域進(jìn)行驗(yàn)證。以模型的反演結(jié)果作為估測(cè)的POM, 以無(wú)人機(jī)影像提取的結(jié)果作為真實(shí)的POM, 分析擬合結(jié)果(圖7)。從表2 可知, 使用綠波段進(jìn)行反演,R2最大, 為0.92; 均方根誤差最小, 平均相對(duì)誤差也最小, 表明使用綠波段所建立的模型效果最好。

表2 藍(lán)、綠、紅波段反演模型結(jié)果比較Tab. 2 Comparison of retrieval results for blue, green and red bands’ retrieval models

2.3 模型適用性分析

對(duì)Landsat 衛(wèi)星影像的陸地部分進(jìn)行掩膜處理, 將綠波段反演公式應(yīng)用于研究區(qū)域, 計(jì)算每個(gè)像素的POM, 衛(wèi)星影像反演的綠藻覆蓋度如圖8a 所示, 無(wú)人機(jī)影像提取結(jié)果反演的綠藻覆蓋度如圖8b 所示。

計(jì)算得到研究區(qū)域的衛(wèi)星反演綠藻覆蓋面積為0.4883km2, 無(wú)人機(jī)影像提取的真實(shí)綠藻覆蓋面積為0.5729km2, 反演面積與真實(shí)面積的比值為0.85。反演結(jié)果雖然對(duì)部分像素綠藻覆蓋度有所低估, 但對(duì)于整體綠藻的分布有較好的反映, 可以應(yīng)用于此區(qū)域綠藻覆蓋面積的估測(cè)。

3 結(jié)論

結(jié)合無(wú)人機(jī)低空航拍影像和Landsat-8 衛(wèi)星影像, 提出了基于NGRDI、NDVI、VB-FAH 和藍(lán)、綠、紅、近紅外波段反射率的綠藻亞像素覆蓋度估算線性模型。

在參與新模型建模的眾多參數(shù)中, 藍(lán)、綠、紅波段反射率與綠藻亞像素覆蓋度呈現(xiàn)出更好的線性關(guān)系, 其R2分別為0.95、0.96、0.93。在驗(yàn)證區(qū)域?qū)诖巳N模型擬合的POM 結(jié)果進(jìn)行驗(yàn)證,其中綠波段反射率所建立的模型擬合結(jié)果最好,其R2為0.92, 均方根誤差為0.07, 平均相對(duì)誤差為10.85%。使用綠波段反射率所建立的模型應(yīng)用于整個(gè)研究區(qū)域所反演的綠藻覆蓋面積與真實(shí)綠藻覆蓋面積(無(wú)人機(jī)圖像所提取的綠藻)比值為0.85。此模型可以較好地反映此研究區(qū)域的綠藻分布情況, 而對(duì)于更廣泛的海水區(qū)域, 則需要進(jìn)一步研究。