閔鈺魁，柯櫻海，韓月，尹小嵐，周德民

1.城市環(huán)境過(guò)程與數(shù)字模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地,北京 100048;

2.水資源安全北京實(shí)驗(yàn)室,北京 100048;

3.首都師范大學(xué) 資源環(huán)境與旅游學(xué)院,北京 100048

1 引言

濱海濕地位于海洋和陸地的交錯(cuò)帶，具有極高的生態(tài)服務(wù)價(jià)值，在維持生物多樣性和調(diào)節(jié)氣候方面發(fā)揮著不可替代的作用（雷昆和張明祥，2005）。互花米草原產(chǎn)自北美西海岸，是生長(zhǎng)在潮間帶的一種鹽沼植物，具有耐鹽、耐淹的特性，可促進(jìn)泥沙沉降和淤積（陳一寧 等，2005）。1979年出于海岸保護(hù)的目的引入中國(guó)（Chung，2006；于冬雪 等，2022）。由于具有強(qiáng)大的適應(yīng)能力和繁殖能力，互花米草在中國(guó)濱海濕地迅速擴(kuò)張（Zuo等，2012），造成土著物種退化，底棲動(dòng)物減少，生物多樣性受到嚴(yán)重威脅。為了遏制互花米草入侵，近年來(lái)，全國(guó)范圍如上海崇明東灘（湯臣棟，2016）、山東黃河口（山東省自然資源廳，2020）、福建閩江河口（周冬，2020）、江蘇鹽城和廣西北海（李飛飛，2021）等濱海濕地先后開(kāi)展了互花米草的清除和試點(diǎn)工程。2020 年6 月國(guó)家發(fā)改委和自然資源部聯(lián)合發(fā)布的《全國(guó)重要生態(tài)系統(tǒng)保護(hù)和修復(fù)重大工程總體規(guī)劃（2021—2035）》將加強(qiáng)互花米草等外來(lái)入侵物種災(zāi)害防治納入海岸帶生態(tài)保護(hù)和修復(fù)重大工程（中華人民共和國(guó)國(guó)家發(fā)展和改革委員會(huì)和中華人民共和國(guó)自然資源部，2020）。

黃河三角洲濕地是中國(guó)三大河口濕地之一和重要的候鳥(niǎo)棲息地（陳建 等，2011）。1989年互花米草引入黃河口，2008 年之后在潮間帶呈現(xiàn)爆發(fā)式擴(kuò)張（Ren 等，2019；Wang 等，2021），擠占本土物種（堿蓬、蘆葦、海草）生態(tài)位，嚴(yán)重威脅濱海濕地生態(tài)系統(tǒng)健康（任武陽(yáng) 等，2019；李昱蓉 等，2021）。2020 年山東省自然資源廳等五部門(mén)聯(lián)合印發(fā)《山東省互花米草防治實(shí)施方案》（山東省自然資源廳，2020），要求在3 年內(nèi)清除省內(nèi)互花米草，對(duì)黃河口等海岸區(qū)域開(kāi)展大規(guī)模的互花米草治理工程，通過(guò)刈割移除和二次翻耕，實(shí)現(xiàn)互花米草清除。及時(shí)準(zhǔn)確地了解互花米草清除動(dòng)態(tài)和清除面積，對(duì)于互花米草治理工程實(shí)施效果評(píng)價(jià)以及濱海濕地管理十分重要。

遙感技術(shù)已成為濱海濕地監(jiān)測(cè)的重要工具，國(guó)內(nèi)外學(xué)者在濱海濕地長(zhǎng)時(shí)序監(jiān)測(cè)方面開(kāi)展了一系列研究（溫慶可 等，2011；張興贏 等，2015；Xu 等，2018；Fernández-Guisuraga 等，2020）。如Lymburner 等（2020）應(yīng)用Landsat 數(shù)據(jù)，監(jiān)測(cè)了1984 年—2016 年澳大利亞濱海濕地紅樹(shù)林動(dòng)態(tài)變化；Wang等（2021）利用Landsat影像，獲取并分析了1984 年—2018 年中國(guó)濱海濕地的年際變化。然而，對(duì)于濱海濕地年內(nèi)動(dòng)態(tài)監(jiān)測(cè)，尤其是入侵植物清除動(dòng)態(tài)監(jiān)測(cè)的研究較為缺乏。目前，植物清除的遙感監(jiān)測(cè)研究主要集中于作物收獲（張峰等，2004；Gao 等，2017；陸俊 等，2018；Shang等，2020）、牧草收割（Schwieder等，2022），森林?jǐn)_動(dòng)（López-Amoedo等，2021；Han等，2021；Zhao等，2022），通過(guò)植物清除后的地表響應(yīng)特征檢測(cè)清除事件，提取清除時(shí)間和空間范圍。如，陸俊等（2018）融合MODIS 和HJ-1B 數(shù)據(jù)，構(gòu)建冬小麥?zhǔn)斋@前后NDVI 時(shí)序，通過(guò)NDVI 曲率閾值法檢測(cè)收割日期；Schwieder 等（2022）基于Sentinel-2和Landsat 影像，通過(guò)牧草時(shí)序NDVI 曲線(xiàn)特征檢測(cè)牧草清除事件，對(duì)德國(guó)2017 年—2020 年間的牧草收割范圍進(jìn)行制圖。與內(nèi)陸地表環(huán)境不同，濱海濕地具有復(fù)雜的水文條件，下墊面長(zhǎng)期處于動(dòng)態(tài)變化的狀態(tài)。尤其在潮間帶，頻繁的潮汐淹沒(méi)會(huì)對(duì)刈割清除事件的特征識(shí)別造成較大干擾，給時(shí)序動(dòng)態(tài)監(jiān)測(cè)帶來(lái)挑戰(zhàn)。同時(shí)，濱海濕地多云多雨，單一平臺(tái)遙感影像無(wú)法構(gòu)建密集的時(shí)序數(shù)據(jù)，難以及時(shí)捕捉清除事件導(dǎo)致的地表快速變化。因此，亟需融合多平臺(tái)衛(wèi)星影像，構(gòu)建融合的密集時(shí)序數(shù)據(jù)集，實(shí)現(xiàn)清除事件的準(zhǔn)確識(shí)別。

針對(duì)以上問(wèn)題，本文旨在提出一種濱海濕地互花米草刈割清除動(dòng)態(tài)的遙感監(jiān)測(cè)方法。以黃河口濕地為研究區(qū)，基于Sentinel-2 和國(guó)產(chǎn)GF-1 影像，融合構(gòu)建密集時(shí)序數(shù)據(jù)，充分考慮植被光譜指數(shù)的時(shí)序變化和潮間帶潮汐淹沒(méi)動(dòng)態(tài)，實(shí)現(xiàn)互花米草清除范圍的快速提取和清除時(shí)間的精準(zhǔn)識(shí)別。

2 研究區(qū)和數(shù)據(jù)

2.1 研究區(qū)概況

黃河口濕地位于山東省東營(yíng)市黃河入?？冢?7°38'45″N—37°55'15″N，119°01'30″E—119°19'45″E（圖1）。屬暖溫帶地區(qū)，以半濕潤(rùn)大陸性季風(fēng)氣候?yàn)橹?。互花米草主要分布于黃河口濕地潮間帶，12月—4月為非生長(zhǎng)季，5月進(jìn)入拔節(jié)期，6月中旬開(kāi)始開(kāi)花，延續(xù)到9 月；10 月—11 月種子成熟，11 月進(jìn)入生長(zhǎng)后期，11 月底—12 月初逐漸枯死（石東里 等，2009；Ren等，2019）。自2008年，互花米草面積從39.91 ha 擴(kuò)張到2019 年的4672.38 ha（Wang 等，2021）。為解決互花米草生態(tài)入侵問(wèn)題，東營(yíng)市政府于2021 年9—12 月在黃河口開(kāi)展了大規(guī)模的互花米草治理工程，利用大型機(jī)械對(duì)互花米草草甸進(jìn)行刈割、翻耕等處理（圖1（b））。本文選取2021 年黃河口濕地互花米草定植區(qū)作為研究區(qū)，監(jiān)測(cè)其清除動(dòng)態(tài)。2021 年互花米草分布范圍邊界基于Wang 等（2021）互花米草邊界的基礎(chǔ)上結(jié)合遙感影像目視解譯獲得（圖1（a））。

圖1 研究區(qū)位置和實(shí)地勘探圖Fig.1 Location of study area and photos of field investigations

2.2 數(shù)據(jù)與數(shù)據(jù)預(yù)處理

2.2.1 Sentinel-2和GF-1數(shù)據(jù)

Sentinel-2 是歐洲航天局研發(fā)的寬幅、高分辨率、多光譜成像的對(duì)地觀測(cè)系統(tǒng)，該系統(tǒng)由Sentinel-2A 和Sentinel-2B 兩顆衛(wèi)星組成，分別于2015 年6 月23 日和2017 年3 月7 日發(fā)射，二者結(jié)合可以實(shí)現(xiàn)5天一次的對(duì)地觀測(cè)。兩顆衛(wèi)星都搭載了多光譜成像儀（MSI），包括13 個(gè)波段，涵蓋可見(jiàn)光、紅邊、近紅外和短波紅外光譜，其中藍(lán)、綠、紅和近紅外4 個(gè)波段為10 m 分辨率，其他波段分別為20 m和60 m分辨率。

高分一號(hào)（GF-1）是中國(guó)高分辨率對(duì)地觀測(cè)系統(tǒng)的第一顆衛(wèi)星，于2013 年4 月26 日發(fā)射。該衛(wèi)星搭載2 個(gè)全色和多光譜傳感器（PMS）和4 個(gè)寬視場(chǎng)傳感器（WFV），其中PMS相機(jī)可以獲取分辨率為2 m 的全色波段和8 m 的多光譜影像，WFV相機(jī)可以獲取由藍(lán)、綠、紅和近紅波段組成的16 m分辨率的多光譜影像。

互花米草治理工程作業(yè)時(shí)間在2021 年9 月初至12 月中旬。根據(jù)衛(wèi)星過(guò)境時(shí)間，本文選取了2020 年和2021 年9 月至12 月之間云量小于30%的19幅Sentinel-2A/B影像，其中2020年8幅，2021年11幅。表1所示，2021年9月13日至10月8日期間沒(méi)有滿(mǎn)足條件的Sentinel-2影像，10月13日至11月12 日期間僅有1 景Sentinel-2 影像。為提高觀測(cè)頻率，選取2021年9—10月5幅云量小于30%的GF-1影像（3 幅GF-1 PMS 和2 幅GF-1 WFV 多光譜影像）作為補(bǔ)充。將9 月8 日作為起始日期，其他日期以此為基準(zhǔn)用天數(shù)來(lái)表示。GF-1 Level-1A 和Sentinel-2 Level-C 數(shù)據(jù)分別來(lái)自中國(guó)資源衛(wèi)星應(yīng)用中心和歐洲航天局網(wǎng)站。

表1 哨兵2號(hào)和高分一號(hào)影像參數(shù)Table1 Sentinel-2 and GF-1 image parameters

2.2.2 數(shù)據(jù)預(yù)處理

Sentinel-2 數(shù)據(jù)和GF-1 數(shù)據(jù)的預(yù)處理包括：輻射校正、云掩膜、重采樣、鑲嵌裁剪和幾何校正。利用Sen2Cor_2.90 軟件對(duì)Sentinel-2A/B Level-C 數(shù)據(jù)進(jìn)行大氣校正，同時(shí)利用QA60 波段對(duì)云像元進(jìn)行掩膜。從中國(guó)資源衛(wèi)星應(yīng)用中心獲取GF-1 PMS/WFV 絕對(duì)輻射定標(biāo)系數(shù)，對(duì)GF-1 數(shù)據(jù)進(jìn)行輻射定標(biāo)，并利用FLAASH 模型進(jìn)行大氣校正；通過(guò)設(shè)置閾值（紅波段地表反射率大于0.2）并結(jié)合目視解譯去除GF-1 影像上的云像元。為與Sentinel-2 數(shù)據(jù)空間分辨率保持一致，將PMS 和WFV 數(shù)據(jù)重采樣到10 m，選取地面控制點(diǎn)并通過(guò)多項(xiàng)式校正模型和三次卷積法內(nèi)插法將GF-1 影像數(shù)據(jù)配準(zhǔn)到Sentinel-2 影像，最后對(duì)每一景影像計(jì)算歸一化植被指數(shù)NDVI（Normalized Difference Vegetation Index）。

2.2.3 實(shí)地調(diào)查數(shù)據(jù)

研究人員于2021 年10 月初和2021 年12 月初開(kāi)展了兩次野外踏勘，利用高精度手持RTK 進(jìn)入研究區(qū)并對(duì)互花米草清除狀態(tài)進(jìn)行實(shí)地調(diào)查，記錄互花米草清除狀況（圖1（b）上面4幅圖片拍攝于2021 年12 月07 日）。實(shí)地勘測(cè)點(diǎn)分布見(jiàn)圖1（a）。對(duì)于難以進(jìn)入的區(qū)域，使用大疆精靈4多光譜無(wú)人機(jī)P4M（Phantom4-M）進(jìn)行勘探（圖1（b）下面2幅圖片拍攝于2021年10月11日），飛行航高30 m，航向重疊90%，旁向重疊80%，采集的圖像空間分辨率為2—4 cm。

3 方法和實(shí)現(xiàn)

本文方法主要包括3個(gè)步驟：（1）基于Sentinel-2 和GF-1 的NDVI 校正和時(shí)序NDVI 構(gòu)建；（2）基于2020 年、2021 年時(shí)序NDVI 差異特征的潛在清除時(shí)段提??；（3）潛在清除時(shí)段中潮汐淹沒(méi)時(shí)刻識(shí)別和清除日期提取。

3.1 時(shí)序NDVI構(gòu)建

由于GF-1 PMS、GF-1 WFV相機(jī)與Sentinel-2A/B MSI傳感器的波譜響應(yīng)和成像條件不同，因此即使同一時(shí)間過(guò)境，GF-1 和Sentinel-2 計(jì)算的NDVI 也會(huì)存在差異。本研究盡量選取同一天（或日期相近）獲取的GF-1 WFV 和Sentinel-2 MSI 無(wú)云影像對(duì)，以及GF-1 PMS 和Sentinel-2 MSI 無(wú)云影像對(duì)，采用線(xiàn)性回歸方法對(duì)GF-1 PMS 和WFV 獲得的NDVI進(jìn)行校正，使其與Sentinel-2的NDVI盡可能一致。經(jīng)過(guò)篩選，選取2020 年9 月18 日獲取的GF-1 WFV 和Sentinel-2 MSI影像；選取2019年8月30日獲取的GF-1 PMS 和同年9 月3 日獲取的Sentinel-2 MSI影像構(gòu)成影像對(duì)。

在研究區(qū)內(nèi)均勻抽取4000個(gè)像元，分別以像元GF-1 PMS、WFV獲取的NDVI為自變量，Sentinel-2 MSI NDVI 為因變量，進(jìn)行最小二乘線(xiàn)性回歸擬合（式（2）），最后利用解算出的回歸系數(shù)對(duì)表1 中GF-1 PMS、WFV影像得到的NDVI圖層進(jìn)行校正。

式中，NDVIS2為Sentinel-2 MSI NDVI，a、b為回歸系數(shù)，NDVIGF-1為GF-1（PMS 或WFV）的NDVI。校正后的GF-1 PMS和WFV NDVI與Sentinel-2 MSI獲得的NDVI 聯(lián)合使用，構(gòu)建密集的NDVI 影像數(shù)據(jù)集合。由于表1中部分影像存在云覆蓋，會(huì)造成部分像元的時(shí)序NDVI 出現(xiàn)缺失。為彌補(bǔ)缺失值，同時(shí)保留NDVI 時(shí)序變化特征，逐像元地利用缺失值前后的無(wú)云觀測(cè)值進(jìn)行線(xiàn)性插值。然后，利用Savitzky-Golay（S-G）濾波逐像元對(duì)NDVI 時(shí)序去噪平滑，最后通過(guò)三次樣條插值方法分別生成各像元2021年9月8日至12月17日逐日的NDVI時(shí)序曲線(xiàn)（NDVI2021）。

3.2 潛在清除時(shí)段提取

3.2.1 突變時(shí)段提取

互花米草刈割移除發(fā)生在9—12月份，此時(shí)互花米草處于成熟季到枯萎季的變化階段，NDVI 呈現(xiàn)逐漸下降的特征。因此，2021 年刈割時(shí)，互花米草NDVI 時(shí)序變化特征同時(shí)受到自然枯落和刈割移除的雙重影響。為了排除前者的影響，引入2020 年的時(shí)序NDVI 作為標(biāo)準(zhǔn)NDVI 時(shí)序曲線(xiàn)。具體而言，首先根據(jù)2020 年各影像計(jì)算互花米草定植區(qū)NDVI 的均值，各時(shí)期NDVI 均值進(jìn)行S-G 濾波和三次樣條插值，生成標(biāo)準(zhǔn)NDVI 時(shí)序曲線(xiàn)（NDVI2020）。NDVI2020代表自然衰落的時(shí)序變化特征（圖2）。對(duì)NDVI2020和NDVI2021做差值運(yùn)算，生成NDVIDIF曲線(xiàn)，通過(guò)該曲線(xiàn)中相鄰兩個(gè)時(shí)間的值做差，生成變化率曲線(xiàn)F，公式如下：

圖2 潛在清除時(shí)段提取示意圖Fig.2 Illustration of potential cutting period extraction

式中，t代表時(shí)間，NDVI2021(t)代表像元2021 年t時(shí)刻的NDVI，NDVI2020(t) 代 表2020 年t 時(shí) 刻N(yùn)DVI，NDVIDIF(t)上升表明該像元在2021 年t時(shí)刻相較于自然衰落具有更強(qiáng)烈的NDVI 下降趨勢(shì)，F(xiàn)是NDVIDIF的變化率曲線(xiàn)。F(t)反映t時(shí)刻N(yùn)DVIDIF的變化趨勢(shì)，F(xiàn)(t)大于0，表明NDVIDIF(t)在（t，t+1）時(shí)間段具有上升的趨勢(shì)，F(xiàn)(t)值越大，上升趨勢(shì)越顯著。F(t)中連續(xù)大于0 的時(shí)段代表NDVIDIF(t)具有持續(xù)的上升趨勢(shì)，將這些時(shí)段作為單獨(dú)的突變時(shí)段提取出來(lái)，表示為（Si，Ti），其中，Si表示第i個(gè)突變時(shí)段的起始日期，Ti表示終止日期（圖2）。通過(guò)識(shí)別突變時(shí)段，去除NDVIDIF時(shí)序變化平穩(wěn)和具有下降趨勢(shì)的時(shí)段，縮小提取范圍。

3.2.2 潛在清除時(shí)段提取

通過(guò)突變時(shí)段提取，獲取了NDVIDIF時(shí)序中所有具有上升趨勢(shì)的時(shí)段。為了表征各突變時(shí)段上升趨勢(shì)的強(qiáng)度，對(duì)各突變時(shí)段F(t)進(jìn)行積分，積分結(jié)果用momentum表示（式（4））?；セ撞葚赘钋宄录l(fā)生會(huì)導(dǎo)致NDVI2021(t)大幅下降，對(duì)各突變時(shí)段起止日期的NDVI2021(t)做差值，計(jì)算下降幅度magnitude（式（5））。

式中，i表示突變時(shí)段的序號(hào)，i=1，2，…，n，n為突變時(shí)段數(shù)量，Si和Ti分別代表第i個(gè)突變時(shí)段的起止日期，momentum(i)反映第i個(gè)突變時(shí)段NDVIDIF(t)（t∈（Si，Ti））上升趨勢(shì)的強(qiáng)度，magnitude(i)表示第i個(gè)各突變時(shí)段內(nèi)的NDVI2021(t)（t∈（Si，Ti））的下降幅度。將滿(mǎn)足momentum(i)＞0.1和magnitude(i)＞0.2的突變時(shí)段認(rèn)定為潛在清除時(shí)段，表明該時(shí)段內(nèi)發(fā)生了引起NDVI 時(shí)序大幅下降的極端事件。如果沒(méi)有滿(mǎn)足該條件的突變時(shí)段，表明該像元NDVI 時(shí)序變化較為平穩(wěn)，未出現(xiàn)劇烈下降趨勢(shì)，我們將該像元認(rèn)定為未發(fā)生刈割清除事件。

3.3 潮汐淹沒(méi)監(jiān)測(cè)與清除日期提取

3.3.1 潮汐淹沒(méi)監(jiān)測(cè)

互花米草生長(zhǎng)在潮間帶，頻繁的潮汐淹沒(méi)會(huì)造成時(shí)序NDVI 產(chǎn)生較大波動(dòng)。黃河口屬于弱潮多沙型河口，互花米草刈割前，由于互花米草植株較高且密集，漲潮水面多被互花米草植株遮蓋，NDVI 會(huì)出現(xiàn)小幅下降?；セ撞萸宄?，漲潮時(shí)潮水會(huì)淹沒(méi)潮灘，此時(shí)在遙感影像上反映的是水體反射率，NDVI 劇烈下降至負(fù)值，這與刈割清除事件反映的時(shí)序特征是相似的，因此會(huì)對(duì)清除日期的識(shí)別產(chǎn)生干擾。此外，互花米草清除后，可能發(fā)生多次潮汐淹沒(méi)。因此，有必要對(duì)潮汐淹沒(méi)進(jìn)行監(jiān)測(cè)識(shí)別。

潮汐淹沒(méi)監(jiān)測(cè)是基于2021 年NDVI 真實(shí)觀測(cè)值，用NDVIObs(tj)表示，j代表觀測(cè)值的序號(hào)（1-16），tj表示第j個(gè)觀測(cè)值對(duì)應(yīng)的天數(shù)，tj=1，6，9，15，…，101（表1）。針對(duì)每個(gè)像元，遍歷從3.2小節(jié)識(shí)別的所有潛在清除時(shí)段中的NDVIObs(tj)，若NDVIObs(tj)＜0，表明像元在tj時(shí)刻發(fā)生了潮汐淹沒(méi)，將該像元?jiǎng)澐譃檠蜎](méi)像元，并將該時(shí)刻命名為淹沒(méi)時(shí)刻；若存在多個(gè)淹沒(méi)時(shí)刻，取日期最早的潮汐淹沒(méi)時(shí)刻，如果不存在淹沒(méi)時(shí)刻，表明該像元在研究時(shí)段內(nèi)未發(fā)生潮汐淹沒(méi)事件，將該像元?jiǎng)澐譃榉茄蜎](méi)像元。

3.3.2 清除日期提取

對(duì)于非淹沒(méi)像元，遍歷潛在清除時(shí)段上所有NDVIObs(tj)，判斷相鄰NDVI 觀測(cè)值的最大差值（NDVIObs(tj)-NDVIObs(tj+1)）是否 大于清除閾值0.15，若滿(mǎn)足則選取差值最大的相鄰觀測(cè)值時(shí)段作為清除時(shí)段（t1，t2），否則判斷為未清除。

對(duì)于淹沒(méi)像元，記錄潮汐淹沒(méi)的次數(shù)和時(shí)刻。遍歷所有潛在清除時(shí)段，在其中尋找相鄰觀測(cè)時(shí)刻N(yùn)DVI 相差大于0.15 的時(shí)段（如圖3），即NDVIObs(tj)-NDVIObs(tj+1)＞0.15，這些時(shí)段可能存在3 種情況（圖3 綠色陰影表示）：（1）時(shí)段內(nèi)未發(fā)生潮汐淹沒(méi)（如圖3（a）中的①時(shí)段）；（2）第一次潮汐淹沒(méi)發(fā)生在該時(shí)段（如圖3（a）中的②時(shí)段）；（3）發(fā)生了潮汐淹沒(méi)，但不是第一次（如圖3（a）中的③時(shí)段）。如果第一次潮汐淹沒(méi)之前已檢測(cè)到符合第1種情況的時(shí)段，則在其中尋找差值最大的時(shí)段，作為清除時(shí)段（如圖3（a）、（b）中紅色標(biāo)注）。否則，說(shuō)明清除和潮汐淹沒(méi)在同一時(shí)段被監(jiān)測(cè)到，符合這種情況的時(shí)段則被判斷為清除時(shí)段（如圖3（c）中紅色標(biāo)注）。由于潮汐淹沒(méi)導(dǎo)致的NDVI＜0 只能出現(xiàn)在互花米草清除后，因此不考慮第3種情況。為了最大程度接近真實(shí)清除日期，取清除時(shí)段中點(diǎn)作為最終的清除日期。

圖3 清除日期提取示意圖Fig.3 Schematic diagram of clearing date identification

4 結(jié)果與分析

4.1 GF-1 NDVI校正系數(shù)

從圖4 可以看出，Sentinel-2 MSI與GF1 WFV，GF-1 PMS 的NDVI 都具有較強(qiáng)的線(xiàn)性相關(guān)性，線(xiàn)性回歸的決定系數(shù)R2分別為0.82 和0.87，通過(guò)線(xiàn)性回歸擬合公式對(duì)GF-1 WFV、PMS 的NDVI 進(jìn)行校正后，可以與Sentinel-2 MSI 的NDVI 具有較好的一致性。

圖4 Sentinel-2 NDVI和GF-1 WFV/PMS NDVI散點(diǎn)圖Fig.4 Scatterplots of Sentinel-2 NDVI and GF-1 WFV/PMS NDVI

經(jīng)過(guò)最小二乘線(xiàn)性回歸擬合，得到GF-1 WFV、PMS的NDVI校正公式分別為

4.2 清除日期精度驗(yàn)證

本文通過(guò)在研究區(qū)內(nèi)建立300 m×300 m的漁網(wǎng)，選取漁網(wǎng)中點(diǎn)作為驗(yàn)證樣點(diǎn)（圖1（a）），根據(jù)實(shí)地觀測(cè)建立目視解譯標(biāo)志，并逐幅對(duì)比影像的光譜、紋理特征，人工解譯記錄樣本點(diǎn)上互花米草的清除時(shí)間。與互花米草清除日期提取結(jié)果進(jìn)行對(duì)比，得到混淆矩陣。如表2所示，互花米草清除日期識(shí)別的總體精度為88.24%，Kappa系數(shù)為0.87；平均制圖精度和用戶(hù)精度分別為86.73%和87.29%。2021 年10 月4 日用戶(hù)精度較低，僅為64.18%，這是由于2021 年9 月30 日和2021 年10 月8 日的影像中，云量集中分布在北岸的互花米草定植區(qū)域，清除日期提取受到較大干擾。2021年12月1 日清除面積小，樣本點(diǎn)少，因而制圖精度較低。其余時(shí)段清除日期監(jiān)測(cè)的制圖精度和用戶(hù)精度均高于73%。

表2 清除日期識(shí)別精度Table 2 Accuracy of the clearance date detection

4.3 互花米草清除過(guò)程

2021年，互花米草定植區(qū)總面積為5189.26 ha，互花米草清除工程從2021年9月8日至2021年12月17 日，總清除面積為4816.35 ha，未清除面積為372.91 ha。未清除的區(qū)域主要分布在黃河口北岸中部近海區(qū)域（圖5（a）中區(qū)域Ⅱ）和南岸沿河道區(qū)域，這些地區(qū)潮溝交錯(cuò)、水文條件復(fù)雜，泥沙較為松軟，不利于大型工程器械進(jìn)入作業(yè)，推測(cè)可能為該區(qū)域未及時(shí)清除的原因。

圖5 互花米草清除日期的空間分布以及相應(yīng)的時(shí)序影像Fig.5 Spatial distribution of Spartina Alterniflora clearance dates and the corresponding time-series imagery

從互花米草清除時(shí)間的空間分布來(lái)看（圖5（a）），清除工程采用“多點(diǎn)分布式”的方式來(lái)進(jìn)行，在清除早期，清除區(qū)域均勻分散在黃河口區(qū)域，再向四周擴(kuò)散清除。清除方向主要是從內(nèi)陸向海岸區(qū)域。同時(shí)，各個(gè)時(shí)段清除區(qū)域呈現(xiàn)出以潮溝為界線(xiàn)分布，表明潮溝對(duì)于互花米草清除工程的空間連通性帶來(lái)較大影響，潮溝的分布應(yīng)該作為清除工程的重要考量因素。

從各個(gè)時(shí)段的互花米草清除面積分布來(lái)看，清除工程主要分兩個(gè)時(shí)段進(jìn)行清除（圖6）。9月8日至10月4日進(jìn)行了第一階段清除，共2333.46 ha。受降水影響，2021 年10 月5 日黃河口突發(fā)特大洪水（中國(guó)水利部，2021），打斷了互花米草清除進(jìn)程。10 月11 日至12 月17 日進(jìn)行第二階段清除，共2482.89 ha，其中11月15日清除面積最大，清除了788.50 ha，12月前已完成的絕大部分區(qū)域的清除。

圖6 各時(shí)期互花米草清除面積Fig.6 Area of the cleared Spartina Alterniflora on each date

4.4 潮汐淹沒(méi)影響分析

本文通過(guò)判斷在NDVI 時(shí)序上是否存在負(fù)值來(lái)判斷像元是否發(fā)生潮汐淹沒(méi)。圖7展示了淹沒(méi)區(qū)的空間分布和淹沒(méi)次數(shù)，淹沒(méi)區(qū)主要分布在靠近海岸和潮溝的區(qū)域，顏色越深表明發(fā)生淹沒(méi)越頻繁，淹沒(méi)區(qū)的總面積為1848.98 ha，占互花米草定植區(qū)面積的35.63%，表明潮汐淹沒(méi)范圍較大。

圖7 淹沒(méi)區(qū)分布Fig.7 Spatial distribution of inundated area

為探究潮汐淹沒(méi)事件對(duì)于清除監(jiān)測(cè)的影響，論證潮汐淹沒(méi)監(jiān)測(cè)的必要性，本文對(duì)于淹沒(méi)區(qū)進(jìn)行對(duì)比試驗(yàn)，分別運(yùn)用潮汐監(jiān)測(cè)和未運(yùn)用潮汐監(jiān)測(cè)算法（即按照非淹沒(méi)區(qū)處理）提取淹沒(méi)區(qū)的清除日期。如圖8所示，若未進(jìn)行潮汐淹沒(méi)監(jiān)測(cè)，識(shí)別的清除日期普遍推后。選取分布淹沒(méi)區(qū)內(nèi)的191個(gè)樣本點(diǎn)分別對(duì)二者提取結(jié)果進(jìn)行精度驗(yàn)證，結(jié)果表明，經(jīng)過(guò)潮汐監(jiān)測(cè)的總體精度為87.95%，未進(jìn)行潮汐監(jiān)測(cè)的結(jié)果總體精度為27.75%，表明潮汐淹沒(méi)會(huì)嚴(yán)重干擾清除日期的提取。通過(guò)對(duì)樣本點(diǎn)上監(jiān)測(cè)清除日期與清除日期參考值做差，將樣本點(diǎn)分為準(zhǔn)確點(diǎn)、延遲點(diǎn)（監(jiān)測(cè)清除日期比參考日期靠后）和提前點(diǎn)。從圖9（a）驗(yàn)證點(diǎn)分布可以看出，準(zhǔn)確提取的樣本點(diǎn)分布在內(nèi)陸區(qū)域，錯(cuò)誤提取的樣本點(diǎn)主要以延遲點(diǎn)為主，多分布在近海岸區(qū)域，表明潮汐淹沒(méi)對(duì)于近海岸影響更為顯著。這主要是因?yàn)楫?dāng)發(fā)生大規(guī)模漲潮，潮汐淹沒(méi)的NDVI 時(shí)序變化特征比清除事件更強(qiáng)烈，從而導(dǎo)致識(shí)別的清除日期延后。

圖8 進(jìn)行潮汐淹沒(méi)監(jiān)測(cè)和未進(jìn)行潮汐淹沒(méi)監(jiān)測(cè)的清除日期提取結(jié)果對(duì)比Fig.8 Comparison of clearance dates with or without consideration of tidal inundation monitoring

圖9 淹沒(méi)區(qū)樣本點(diǎn)分布圖和典型樣本點(diǎn)時(shí)序NDVI示意圖Fig.9 Distribution of sample points in the tidal inundated area and illustration of NDVI time series at a sample

4.5 考慮NDVI2020對(duì)互花米草清除監(jiān)測(cè)的作用

本文以2020 年時(shí)序NDVI（NDVI2020）為基準(zhǔn)，代表互花米草未清除時(shí)NDVI 變化特征，目的是為了排除互花米草自然枯萎對(duì)清除識(shí)別的影響。為探討使用2020 年基準(zhǔn)NDVI 對(duì)于清除監(jiān)測(cè)的重要性，本研究設(shè)置一個(gè)實(shí)驗(yàn)，單獨(dú)使用2021 年的時(shí)序NDVI 提取互花米草清除日期，即用NDVI2021代替式（2）中的NDVIDIF，實(shí)驗(yàn)結(jié)果如圖10所示。

圖10 考慮NDVI2020和未考慮NDVI2020進(jìn)行清除日期識(shí)別的結(jié)果對(duì)比Fig.10 Comparison of clearance dates with or without consideration of NDVI2020

圖10（a）為研究區(qū)未清除互花米草的主要分布區(qū)域，Sentinel-2 影像為2021 年12 月17 日遙感假彩色合成影像。結(jié)果表明，若不考慮NDVI2020，未被清除的互花米草容易被誤判為已清除，此時(shí)研究區(qū)未清除的面積僅為151.59 ha。使用NDVI2020作為輔助數(shù)據(jù)的主要作用是抵消互花米草自然枯萎對(duì)NDVI 時(shí)序的影響，凸顯清除事件的特征。僅利用NDVI2021，會(huì)導(dǎo)致提取的潛在清除時(shí)段范圍過(guò)大。同時(shí)由于觀測(cè)點(diǎn)時(shí)間分布不均勻，部分時(shí)間間隔較大的NDVI 差值可能會(huì)大于清除閾值，從而造成誤判。對(duì)于未被清除的互花米草區(qū)域，NDVI2021與NDVI2020十分相近，因此NDVIDIF曲線(xiàn)變化平穩(wěn)，誤判為清除的概率大大降低。

4.6 引入GF-1數(shù)據(jù)的必要性分析

本文融合Sentinel-2 和GF-1 數(shù)據(jù)構(gòu)建密集時(shí)序，對(duì)互花米草清除治理進(jìn)行動(dòng)態(tài)監(jiān)測(cè)。為了探討GF-1 數(shù)據(jù)在本方法必要性，本文假定單獨(dú)使用Sentinel-2 時(shí)序數(shù)據(jù)，應(yīng)用該算法提取互花米草清除日期。實(shí)驗(yàn)結(jié)果如圖11 所示。為了便于二者對(duì)比，本實(shí)驗(yàn)取3.3.2 中清除時(shí)段的t2時(shí)刻，即該時(shí)段第二景影像獲取日期作為清除日期。

圖11 單獨(dú)使用Sentinel-2數(shù)據(jù)和使用Sentinel-2、GF-1融合數(shù)據(jù)的清除日期提取結(jié)果Fig.11 Comparison of data extraction results based on Sentinel-2 and fusion data

圖11（Ⅰ）、（Ⅲ）為融合數(shù)據(jù)的結(jié)果，（Ⅱ）、（Ⅳ）為僅使用Sentinel-2 數(shù)據(jù)的結(jié)果，虛線(xiàn)黑框內(nèi)代表GF-1 影像，其余影像為Sentinel-2 影像。結(jié)果表明，若僅使用Sentinel-2 數(shù)據(jù)，有部分區(qū)域識(shí)別日期推后，如圖11黃色橢圓邊框內(nèi)，在GF-1影像顯示已經(jīng)清除。Sentinel-2 和GF-1 結(jié)合使用較好地規(guī)避了由于單一數(shù)據(jù)源時(shí)間分布不均的問(wèn)題，通過(guò)在Sentinel-2 時(shí)間間隔較大的時(shí)段填補(bǔ)GF-1 數(shù)據(jù)，使時(shí)序信息更加完整，可以更加及時(shí)地監(jiān)測(cè)清除過(guò)程，提高監(jiān)測(cè)精度。

5 結(jié)論

本文以黃河口濕地為研究區(qū)，基于Sentinel-2和GF-1 時(shí)序影像，提出了一種濱海濕地互花米草清除動(dòng)態(tài)監(jiān)測(cè)方法，能夠克服潮汐淹沒(méi)的影響，快速準(zhǔn)確地提取互花米草刈割清除的時(shí)間。主要結(jié)論如下：（1）線(xiàn)性擬合表明，GF-1 WFV/PMS NDVI 與Sentinel-2 MSI NDVI 具有較好的線(xiàn)性關(guān)系（R2分別為0.82，0.87）；通過(guò)回歸系數(shù)，可以使GF-1 與Sentinel-2 NDVI 保持較好的一致性，有利于時(shí)序一致性NDVI 曲線(xiàn)的構(gòu)建。（2）本方法提取的清除日期總體精度達(dá)到88.24%，Kappa 系數(shù)為0.87。研究區(qū)屬于潮間帶，頻繁的潮汐淹沒(méi)導(dǎo)致NDVI 劇烈波動(dòng)，對(duì)互花米草清除監(jiān)測(cè)造成較大影響。本文提出的潮汐淹沒(méi)識(shí)別方法可以有效避免該問(wèn)題，顯著提高互花米草清除監(jiān)測(cè)精度。（3）融合Sentinel-2 和GF-1 數(shù)據(jù)相較僅使用Sentinel-2 數(shù)據(jù)源，可以構(gòu)建密集、均勻的時(shí)序，更準(zhǔn)確地監(jiān)測(cè)互花米草清除動(dòng)態(tài)變化。（4）2021 年黃河口互花米草清除面積為4816.35 ha，占互花米草分布總面積的92.81%。未清除區(qū)域主要分布于水文復(fù)雜，潮溝交錯(cuò)的北岸靠海區(qū)域。（5）黃河口互花米草清除工程分兩階段完成，第一段在9 月初到10 月初，第二段在10月中到12月中，10月上旬黃河下游洪峰打斷了施工進(jìn)度。本文提出的方法對(duì)互花米草清除監(jiān)測(cè)可以達(dá)到較高精度，但仍依賴(lài)于本地運(yùn)算資源，在大尺度的互花米草治理監(jiān)測(cè)的應(yīng)用中受到限制。隨著中國(guó)各沿海省份陸續(xù)開(kāi)展互花米草治理，下一步的研究將致力于改進(jìn)算法并應(yīng)用于更大范圍治理工程的動(dòng)態(tài)監(jiān)測(cè)和治理效果的綜合評(píng)估中。