孟 黎，孟 靜

（1.山東城市建設(shè)職業(yè)學(xué)院，山東 濟南 250103；2.山東省國土測繪院，山東 濟南 250102）

1 引言

在內(nèi)陸水域中，湖泊富營養(yǎng)化是影響其自身及周邊生態(tài)環(huán)境最嚴(yán)重的問題之一，葉綠素a（Chla）作為浮游植物攜帶的重要色素[1]，其濃度含量會影響浮游植物生物量，也可能改變水庫、湖泊、河流等水體的初級生產(chǎn)力與富營養(yǎng)化程度，因此定量評價內(nèi)陸水體中Chla 濃度含量，描述水體光學(xué)性質(zhì)與水質(zhì)參數(shù)之間的關(guān)系具有重要意義[2]。相較于傳統(tǒng)手段，遙感技術(shù)可以在宏觀尺度上實時對目標(biāo)進行監(jiān)測，節(jié)省了大量人力物力財力[3]。對于Chla 濃度反演，國內(nèi)外學(xué)者目前已作出重要突破，如朱廣偉等[4]基于站點數(shù)據(jù)反演了太湖長時序Chla 濃度變化，并對其驅(qū)動因素作出分析。馬榮華等[5]所提出的OCx算法，目前已經(jīng)成為NASA 水色估算中的默認(rèn)算法。徐京萍等[6]使用Modis 數(shù)據(jù)，對太湖藍藻水華情況進行監(jiān)測，為相關(guān)部門提供決策支持。此類研究中，數(shù)據(jù)源往往集中于Modis、Landsat、HJ1 等中低分辨率衛(wèi)星影像，這些數(shù)據(jù)在Chla 反演過程中時空分辨率過低，導(dǎo)致在中小型湖泊上不適用。而于2015 年發(fā)射的哨兵二號衛(wèi)星具有高時空分辨率[7]，在中小型湖泊Chla 反演中的應(yīng)用前景探討相對較少[8]。隨著計算機技術(shù)的發(fā)展，機器學(xué)習(xí)通過算法本身改善自變量與因變量的關(guān)系從而可以解決非線性問題，為水質(zhì)參數(shù)反演提供了新思路。對于一些常見的機器學(xué)習(xí)算法，如決策樹、支持向量機、多層感知機等，國內(nèi)外眾多學(xué)者已作出相應(yīng)研究[9]，而鮮有研究考慮集成學(xué)習(xí)算法在內(nèi)陸水體Chla 反演中的適用性。鑒于此，本文使用哨兵二號多光譜數(shù)據(jù)聯(lián)合2 種集成學(xué)習(xí)算法，開展內(nèi)陸水體Chla反演算法研究。

2 研究區(qū)概況與數(shù)據(jù)來源

2.1 研究區(qū)概況

南四湖是淮河流域第二大淡水湖，位于中國山東省南部微山縣，全湖面積1 266 km2，由微山湖、昭陽湖、獨山湖及南陽湖串聯(lián)而成，是山東省第一大湖，也是中國南水北調(diào)東線工程中的重要調(diào)水區(qū)。南四湖水環(huán)境健康運營對其周邊濕地系統(tǒng)生態(tài)穩(wěn)定至關(guān)重要，因此本文以對南四湖水體進行精準(zhǔn)大范圍遙感水質(zhì)監(jiān)測為手段，開展南四湖水體Chla 濃度反演算法研究。研究區(qū)概況，如圖1所示。

圖1 南四湖研究區(qū)概況

2.2 哨兵二號數(shù)據(jù)

哨兵二號A 星于2015 年發(fā)射，B 星于2017 年發(fā)射，實現(xiàn)衛(wèi)星雙星組網(wǎng)，重訪周期為5 d，搭載MSI傳感器共有13個波段，空間分辨率分10、20、60 m。本文采用哨兵二號A 星L2A 級產(chǎn)品，該產(chǎn)品于歐空局哥白尼數(shù)據(jù)中心下載（https：//scihub.copernicus.eu），已經(jīng)過了嚴(yán)格的幾何校正與Sen2Res 大氣校正，可直接用于一些遙感定量反演研究。

2.3 光譜及水質(zhì)濃度數(shù)據(jù)

2021 年9 月10—12 日 和9 月21 日在南 四 湖布設(shè)126 個采樣點，其中9 月10—12 日僅測定Chla 濃度，9 月21 日測量Chla 濃度與水面光譜信息。采樣期間，天氣狀況良好，天空無云、無風(fēng)，野外光譜測定條件達標(biāo)。實測點位，如圖1 所示。水面光譜信息通過美國ASD 公司生產(chǎn)的FieldSpec 4 Hi-Res 便攜式地物光譜儀測定，每個采樣點分別測量5 次光譜曲線求平均值以減小誤差。將實測數(shù)據(jù)帶入以下公式計算，獲取湖面遙感反射率光譜曲線[10]。

式中：Rrs為湖面遙感反射率；Lw為純水遙感反射率；r為水與天空光的反射比常數(shù)，一般取0.025；Lsky為天空光遙感反射率，取0.99；π取3.14；ES為灰板輻射信號；ρp為灰板反射率。

Chla 濃度采用美國安諾ChloroTech 121A 型手持式葉綠素測定儀測定。126 個采樣點Chla 濃度分布，如圖2所示。微山湖采樣點較多，濃度與獨山湖相比較低，總體數(shù)據(jù)均值為32.86 μg/L，標(biāo)準(zhǔn)差為17.65 μg/L，無異常值出現(xiàn)。

圖2 采樣點Chla濃度分布

2.4 衛(wèi)星數(shù)據(jù)處理

下載2021 年9 月11—21 日 哨 兵二號L2A 級產(chǎn)品，影像無云，質(zhì)量良好。首先，在SNAP 軟件中對所有波段進行重采樣處理，將其空間分辨率采樣至10 m；然后，輸出為ENVI 格式，在ENVI 軟件中對所有波段進行合成，使用中國山東省微山縣矢量數(shù)據(jù)對影像進行裁剪，并作反射率歸一化處理；最后，進行影像拼接，以NDWI 法提取水體邊界[11]，并基于實測點位提取影像光譜信息。由于9月10—12日實測點位較多，因此使用9月11日影像反演，9月21日光譜信息僅用于增加樣本點數(shù)量。

3 研究方法和模型因素選擇

3.1 遙感反射率校正

由于哨兵二號L2A級反射率產(chǎn)品本質(zhì)上屬于地表反射率產(chǎn)品，對于水色遙感而言，嚴(yán)格意義上應(yīng)當(dāng)使用遙感反射率進行計算。因此，參考劉瑤等[12]的方法進行遙感反射率校正，由于認(rèn)為內(nèi)陸水體在短波紅外（SWIR）的信號很小，所以從可見光和近紅外波段減去短波紅外波段的最小值，然后除以π，實現(xiàn)地表反射率到遙感反射率的轉(zhuǎn)換。

式中：Rrs為影像地表反射率；Rsr為轉(zhuǎn)換后的遙感反射率；Rswir為所有短波紅外波段；π取3.14。

3.2 GBDT算法

GBDT（Gradient Boosting Decision Tree）是Boosting中的代表性算法，它既是當(dāng)代強力算法XGBoost、LGBM 等算法的基石[13]，也是實際應(yīng)用場景中最穩(wěn)定的算法之一。GBDT 中上一個弱評估器的輸出結(jié)果會影響下一個弱評估器的計算過程，其基本核心思想為：依據(jù)上一個弱評估器的結(jié)果，計算損失函數(shù)，并使用損失函數(shù)自適應(yīng)地影響下一個弱評估器的構(gòu)建。集成模型輸出的結(jié)果，受到整體所有弱評估器的影響。

3.3 XGBoost算法

XGBoost（EXtreme Gradient Boosting）是2014 年由中國學(xué)者陳天奇[14]提出的，是基于GBDT 升級的新一代算法。XGBoost 使用估計貪婪算法、平行學(xué)習(xí)算法、分位數(shù)草圖算法創(chuàng)造了全新的建樹流程；使用感知緩存訪問技術(shù)與核外計算技術(shù)提升算法在硬件上的計算性能；引入Dropout 技術(shù)，為整體建樹流程增加隨機性，其基層樹模型可以很好地擬合非線性數(shù)據(jù)。

3.4 模型輸入特征選擇

選取校正后的哨兵二號L2A數(shù)據(jù)的可見光及近紅外共9個波段作為集成學(xué)習(xí)算法的輸入變量。此外，加入4 種波段反射率比值，分別為藍/綠、紅/綠、近紅/綠、近紅/紅，其中藍、綠、紅、近紅分別對應(yīng)哨兵二號的第2、3、4、8 波段。共計13 個輸入變量，輸出為Chla濃度。

3.5 評估指標(biāo)選取

選取決定系數(shù)（R2）與均方根誤差（RMSE）評估所有模型在全部波段選擇策略上的泛化能力。2 個評價指標(biāo)計算公式如下：

式中：n為樣本數(shù)量；yi為實測數(shù)據(jù);yj為模型預(yù)測值為實測數(shù)據(jù)平均值為模型預(yù)測值平均值。R2越大，RMSE越小，模型精度越高。

4 結(jié)果與分析

4.1 遙感反射率校正結(jié)果

采用式（1）計算各實測點位遙感反射率，取400～900 nm 形成光譜曲線。按上文提到的方法，對哨兵二號實測點位提取光譜進行遙感反射率校正。由于水質(zhì)反演是使用可見光及近紅外波段，提取400～900 nm 范圍內(nèi)的實測光譜、哨兵二號原始光譜及校正后的前9個波段對比，如圖3所示。

圖3 實測光譜與遙感反射率校正光譜對比

從圖3 可以看出，在藍光波段及670 nm 處均有吸收峰，570 nm 附近的反射峰是由于葉綠素和胡蘿卜素的弱吸收以及細(xì)胞散射形成的，該反射峰值與色素組成有關(guān)，可以作為葉綠素定量的標(biāo)志。685～715 nm 處反射峰的出現(xiàn)是含藻類水體最明顯的光譜特征，該反射峰的位置和數(shù)值是Chla 濃度的指示，其出現(xiàn)原因是由于水體和Chla 在此處的吸收系數(shù)達到最小。所以，經(jīng)校正后的反射率不僅保留了原始的Chla 濃度反射特征，而且更加貼合于實測光譜。因此，可以認(rèn)為所選擇的遙感反射率校正方法是有效的。

4.2 模型反演結(jié)果

本文模型的構(gòu)建、訓(xùn)練、調(diào)參、測試均在Python與Anaconda的集成開發(fā)環(huán)境中完成，GBDT模型已在Scikit-learn庫中提供方法，而XGBoost模型使用其原生代碼所提供的Scikit-learn API接口實現(xiàn)，主要調(diào)試參數(shù)包括n_estimators、max_depth、learning_rate、subsample等。使用KFold五折交叉驗證的平均得分，評估模型理論泛化能力。2種模型反演結(jié)果，如圖4所示。

圖4 基于實測點位的模型反演結(jié)果

從圖4 可以看出，當(dāng)13 個特征輸入2 種模型時，2 種模型均具有較強的魯棒性。五折交叉驗證的決定系數(shù)R2在XGBoost 模型達到最高（0.723 5）；均方根誤差出現(xiàn)類似情況，在XGBoost 模型達到最低（9.168 1 μg/L）。經(jīng)觀察發(fā)現(xiàn)，Chla 濃度值為20～40 μg/L 時，2 種模型均產(chǎn)生了高估，結(jié)合圖2，認(rèn)為這是由于處于這個濃度的訓(xùn)練數(shù)據(jù)較少，模型學(xué)習(xí)不充分，從而產(chǎn)生了Chla 濃度值的高估?？傮w而言，XGBoost 模型的精度在基于實測數(shù)據(jù)建模中達到最高，后文將把2種模型應(yīng)用于遙感影像，進一步探討二者在遙感影像上的泛化能力。

4.3 哨兵影像反演結(jié)果

遙感影像上的反演結(jié)果，如圖5 所示。Chla 濃度高低分布狀況大體一致，獨山湖明顯高于微山湖。通過實地考察得知，在實地測量前，微山湖經(jīng)過了大量放水，流向為自南向北，所以獨山湖高于微山湖。我們的結(jié)果與Zhang 等[15]的研究結(jié)果一致，因此可以認(rèn)為本文結(jié)果是準(zhǔn)確的。在GBDT 模型反演結(jié)果中，獨山湖Chla 濃度反演值幾乎全處于40 μg/L 以上，這不符合實測數(shù)據(jù)的情況，所以GBDT 模型存在明顯的Chla 濃度高估情況。此外，結(jié)合圖2 與圖1，就下級湖微山湖而言，GBDT 模型也高估了Chla 濃度，因此可以認(rèn)為XGBoost 模型在哨兵二號數(shù)據(jù)上更具魯棒性，反演結(jié)果更加符合實際情況。

圖5 哨兵二號遙感影像反演結(jié)果

5 結(jié)論

本文以山東省微山縣南四湖為研究背景，使用歐空局提供的Sentinel-2A 影像數(shù)據(jù)及實測數(shù)據(jù)，選取影像前9 個波段及4 種波段比值構(gòu)建了Chla 濃度反演的13 個特征波段，在此基礎(chǔ)上使用NDWI 法進行水體提取、光譜提取及反演模型構(gòu)建，得到以下結(jié)論：經(jīng)過遙感反射率校正的哨兵二號影像與實測光譜更具一致性，更適合用于水質(zhì)參數(shù)反演。XGBoost、GBDT 模型可以用于南四湖水質(zhì)參數(shù)反演，XGBoost 模型在實測數(shù)據(jù)及影像反演上均具有較強的魯棒性，反演結(jié)果與實際情況更加一致。后期研究將會嘗試將該模型應(yīng)用于長時序水質(zhì)參數(shù)反演。