林騰 陳哲鋒 楊義煒

（福建省地質(zhì)測繪院 福建福州 350011）

0 引言

葉綠素a 作為水體中藻類植物和浮游植物最豐富的色素之一，其濃度大小可以在一定程度上反映水體的富營養(yǎng)化程度，是進行水體營養(yǎng)狀態(tài)評價的重要參數(shù)之一［1-2］。對葉綠素a濃度進行監(jiān)測是保護水生態(tài)系統(tǒng)的有效手段之一，是人類社會與自然相互協(xié)調(diào)發(fā)展的必然趨勢。傳統(tǒng)葉綠素a 濃度監(jiān)測主要以地面水站監(jiān)測或?qū)嵉厮畼硬杉O(jiān)測為主，該方法可以獲得監(jiān)測站點或水樣采集點的較高精度的水質(zhì)情況，但受限于監(jiān)測范圍，無法對整個水面進行監(jiān)測，且該方法以人工為主，工作量大，時效性不高，難以滿足實時、動態(tài)和廣泛的需求［3］。遙感技術(shù)的實時性、連續(xù)性、大范圍、低成本等優(yōu)勢，正好彌補了傳統(tǒng)監(jiān)測手段的缺點，成為目前葉綠素a 濃度監(jiān)測最有效的手段之一［4］。

隨著遙感傳感器技術(shù)的進步和水體光譜特征研究的不斷深入，葉綠素a 濃度遙感監(jiān)測從定性發(fā)展到定量，模型算法不斷改進，監(jiān)測精度也不斷提高［5］。國內(nèi)外學(xué)者將遙感技術(shù)應(yīng)用于葉綠素a 濃度監(jiān)測的研究開始于20 世紀70 年代，GORDON等［6］利用蒙特卡羅模擬構(gòu)建固有光學(xué)特性與表觀光學(xué)特性之間的一般關(guān)系，為水質(zhì)遙感的發(fā)展奠定基礎(chǔ)；KEINER 等［7］基于Landsat TM 數(shù)據(jù)，利用神經(jīng)網(wǎng)絡(luò)對葉綠素和懸浮物濃度進行建模反演，取得不錯的效果；郭宇龍等［8］以標準三波段算法為基礎(chǔ)，結(jié)合實測葉綠素a 濃度和光譜數(shù)據(jù)，構(gòu)建了三峽水庫、巢湖、洞庭湖和太湖的葉綠素a 模型，結(jié)果表明COCI 三波段算法具有較強的應(yīng)用潛力；朱利等［9］利用GF-1 WFV 和HJ-1A CCD 數(shù)據(jù)，對太湖進行葉綠素a、懸浮物、透明度監(jiān)測和和富營養(yǎng)化評價，研究表明了2 種數(shù)據(jù)的水質(zhì)反演結(jié)果的一致性，驗證了GF-1 WFV 數(shù)據(jù)的水質(zhì)監(jiān)測潛力。

古田水庫是古田縣人工湖，又名“翠屏湖”，是福建省級風景名勝之一，以發(fā)電為主，兼顧旅游業(yè)，與古田縣經(jīng)濟發(fā)展和生態(tài)環(huán)境息息相關(guān)［10］。近年來，水庫周邊城鎮(zhèn)化進程日益加快，工農(nóng)漁業(yè)生產(chǎn)迅速發(fā)展，生活污水、工業(yè)廢水、農(nóng)業(yè)污染的排放嚴重影響了水體水質(zhì)，導(dǎo)致水庫水質(zhì)每況愈下。目前，針對古田水庫水環(huán)境遙感監(jiān)測的研究還相對較少，本文基于Sentinel-2 衛(wèi)星遙感影像數(shù)據(jù)，結(jié)合同步實測葉綠素a 濃度數(shù)據(jù)對2021 年古田水庫每季度的葉綠素a 參數(shù)開展遙感反演建模，并對結(jié)果進行時空特征分析，為古田水庫水質(zhì)監(jiān)測和治理提供參考依據(jù)。

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

1.1 研究區(qū)概況

古田水庫位于福建省古田縣中部閩江支流古田溪上，水庫控制流域面積1 325 km2，總庫容6.47 億m3，地理位置：118°46′E～118°52′E，26°34′N～26°43′N。研究區(qū)屬于亞熱帶海洋性季風氣候，氣候溫暖濕潤，年平均氣溫在14～20 ℃之間，雨量充沛，3～6 月為梅雨季節(jié)，7～9 月多發(fā)臺風。庫區(qū)水資源豐富，水質(zhì)保護目標為Ⅲ類，主要有玉源溪、柏源溪、芹溪、達才溪、橫洋溪等支流流入水庫。水庫上游分布有平湖、新舫、橋洋、大橋、吉巷等多個鄉(xiāng)鎮(zhèn)，既有以食用菌生產(chǎn)加工為主的工業(yè)生產(chǎn)，又有農(nóng)業(yè)耕種和漁業(yè)養(yǎng)殖，多方面因素對庫內(nèi)水體水質(zhì)產(chǎn)生較大影響，水環(huán)境保護形勢嚴峻。

1.2 數(shù)據(jù)來源與預(yù)處理

本文采用的數(shù)據(jù)為歐空局的5 景L1C 級哨兵2 號（Sentinel-2）多光譜衛(wèi)星遙感影像，覆蓋光譜波段13 個，空間分辨率分別為10、20、60 m，獲取時間分別為2021 年3、6、9、12 月和2022 年1 月，幅寬290 km，完整覆蓋古田水庫。

由于L1C 級數(shù)據(jù)是經(jīng)過幾何精校正的正射影像，并沒有進行輻射定標和大氣校正，在開展古田水庫葉綠素a 濃度反演之前，要先對Sentinel-2 數(shù)據(jù)進行輻射定標、大氣校正和影像裁切等預(yù)處理。輻射校正的精度直接影響影像的質(zhì)量，本文借助SNAP 平臺和歐空局（ESA）發(fā)布的Sen2Cor 插件對原始數(shù)據(jù)進行輻射定標和大氣校正處理，并將空間分辨率統(tǒng)一重采樣為10 m，再對處理后的數(shù)據(jù)進行裁切，最終獲得研究區(qū)5 個時相的Sentinel-2 多光譜影像數(shù)據(jù)。

1.3 實測葉綠素a 濃度數(shù)據(jù)及處理

本文通過實地水質(zhì)采樣來測得古田水庫葉綠素a 濃度數(shù)據(jù)，采樣時間為2022 年1 月13 日—14 日的上午10 時至下午13 時，與Sentinel-2 影像（時相2022 年1 月13 日）實現(xiàn)準同步。研究區(qū)共采集樣本20 組，剔除3 組異常數(shù)據(jù)，在剩余數(shù)據(jù)中隨機選取13 組數(shù)據(jù)用于建立模型，其他4 組數(shù)據(jù)用作模型驗證（圖1）。

圖1 古田水庫采樣點分布圖

2 研究方法

目前，遙感技術(shù)手段在葉綠素a 濃度監(jiān)測領(lǐng)域的應(yīng)用已經(jīng)相對成熟，各種類型的反演算法相繼被提出，常用的可歸為3類：經(jīng)驗?zāi)Ｐ退惴ā虢?jīng)驗/半分析模型算法和分析模型算法。經(jīng)驗?zāi)Ｐ退惴ㄊ窃诖罅康臉颖净A(chǔ)上，基于多光譜遙感數(shù)據(jù)的波段反射率與實測水質(zhì)參數(shù)的相關(guān)性統(tǒng)計分析，篩選最優(yōu)的波段或波段組合來構(gòu)建水質(zhì)參數(shù)的回歸反演模型；半經(jīng)驗?zāi)Ｐ退惴ㄊ且劳袑崪y水面光譜和水質(zhì)參數(shù)的基礎(chǔ)上，利用數(shù)理統(tǒng)計等手段對兩者進行反演建模；分析模型算法是在輻射傳輸理論的基礎(chǔ)上，計算水質(zhì)參數(shù)對光譜特定波段的吸收系數(shù)和散射系數(shù)，進而模擬水質(zhì)參數(shù)據(jù)濃度［11］?？紤]到研究區(qū)僅有少量實測葉綠素a 參數(shù)數(shù)據(jù)，缺少同步實測光譜數(shù)據(jù)和吸收、散射等光學(xué)數(shù)據(jù)，本文采用支持向量回歸機（Support Vactor Regression，SVR）模型進行葉綠素a 參數(shù)反演。SVR 是支持向量機在函數(shù)回歸中的應(yīng)用，它是1 個非線性擬合函數(shù)的過程，有效解決了樣本數(shù)量不足的問題，一定程度上也可以提高葉綠素a 濃度反演精度［12］。SVR 葉綠素a 定量反演是以葉綠素a 參數(shù)敏感波段為輸入樣本，利用交叉驗證法尋找最佳模型參數(shù)來構(gòu)建模型，并利用驗證點對模型進行驗證。

2.1 水體信息提取

水質(zhì)遙感反演的對象主要為水體水域。為了提高反演效率，避免植被、建筑、陸地等非水體數(shù)據(jù)參與反演計算，需要對古田水庫進行水域信息提取。目前，常用的遙感水體信息提取方法較多，包括單波段法、多波段法、植被指數(shù)法和水體指數(shù)法等。本文采用歸一化水體指數(shù)方法來提取水域邊界信息。歸一化水體指數(shù)（NDWI）［13］主要依據(jù)的原理是可見光至近紅外波段內(nèi)水體與非水體地物間存在的反射光譜差異，利用水體反射最強的近紅外波段與水體反射最弱的綠波段的比值運算，進一步加大兩者間的差異，最大程度地削減非水體信息而增強水體信息。NDWI 的計算如式（1）。

式中：Green 為綠波段反射率，對應(yīng)Sentinel-2 影像第3 波段；NIR 為近紅外波段反射率，對應(yīng)Sentinel-2 影像第8 波段。

經(jīng)上述運算處理后，可根據(jù)NDWI 圖像直方圖設(shè)置合適的閾值進行水體與非水體的劃分，以提取水體信息。

2.2 相關(guān)性分析

在水質(zhì)遙感監(jiān)測中，一般通過分析單波段或波段組合與水質(zhì)參數(shù)濃度的相關(guān)性，尋找各水質(zhì)參數(shù)的敏感性波段。相關(guān)性分析是探究2 個變量之間的相互關(guān)系而建立的分析方法，分析結(jié)果可以反映2 個變量之間的密切程度，相關(guān)性系數(shù)的數(shù)值范圍是從-1 到1 之間，絕對值越接近1，則兩者的相關(guān)度越大［14］。本文采用Pearson 相關(guān)系數(shù)法進行相關(guān)性分析，其計算如式（2）。

式中：xi和分別為葉綠素a 的濃度及其算術(shù)平均值，yi和y分別為波長λ 的光譜反射率及其算術(shù)平均值；R（λ）為X 和Y之間在波長λ 處的相關(guān)系數(shù)。

2.3 模型構(gòu)建

SVR 模型的構(gòu)建是在Matlab 軟件平臺中利用林智仁教授開發(fā)的libsvm 程序庫，以訓(xùn)練樣本的水質(zhì)敏感波段和水質(zhì)參數(shù)濃度為輸入變量，選擇徑向基核函數(shù)為核函數(shù)，通過調(diào)節(jié)懲罰系數(shù)c、核函數(shù)系數(shù)g 和不敏感系數(shù)p 等模型參數(shù)來尋找效果最優(yōu)的模型［15］。擬合度和均方根誤差是判斷模型是否最優(yōu)的主要依據(jù)。

3 反演模型的建立

3.1 波段敏感性分析

本文采用Pearson 相關(guān)分析法，對實測葉綠素a 濃度與2022 年1 月13 日獲取的Sentinel-2 中400～900 nm 波段范圍內(nèi)各波段（對應(yīng)Sentinel-2 的B1～B8 波段）的光譜反射率進行相關(guān)性分析，分析結(jié)果見表1。從表1 可以看出，葉綠素a 濃度與B3、B5-B8 的相關(guān)系數(shù)在0.6 以上，與B5 的相關(guān)系數(shù)最大，為0.720 2。

表1 Sentinel-2 單波段與葉綠素a 參數(shù)的相關(guān)性

根據(jù)黃靈光等［16］研究發(fā)現(xiàn)，遙感影像各波段之間存在信息冗余的情況，不能很好地反映其與水質(zhì)參數(shù)的關(guān)系。采用波段組合可以在一定程度上消除各波段間的干擾，有效減少水中其他雜質(zhì)的影響，突出水質(zhì)參數(shù)與水體反射率的相關(guān)性。結(jié)合前人經(jīng)驗總結(jié)和波段相關(guān)性分析結(jié)果，選擇Sentinel-2 中與葉綠素a 相關(guān)的B3、B4 和B5 波段進行比值、差值等組合運算，再通過波段組合與葉綠素a 濃度進行相關(guān)性分析，分析結(jié)果見表2。從表2 可以看出，葉綠素a 濃度與波段組合間的相關(guān)性明顯高于單波段，葉綠素a 與（B5-B4）/B3 的相關(guān)性最高，相關(guān)系數(shù)達0.876 2。最終本文選擇（B5-B4）/B3 波段組合為葉綠素a 參數(shù)的敏感波段。

表2 Sentinel-2 波段組合與葉綠素a 參數(shù)的相關(guān)性

3.2 SVR 反演模型建立

本文根據(jù)葉綠素a 參數(shù)敏感性分析結(jié)果，選擇（B5-B4）/B3作為葉綠素a 參數(shù)反演模型構(gòu)建的輸入變量，以徑向基核函數(shù)為核函數(shù)，并采用交叉驗證法來尋找擬合效果最佳的模型參數(shù)，最終確定葉綠素a 參數(shù)的最優(yōu)SVR 反演模型。葉綠素a 參數(shù)的SVR 模型擬合效果如圖2 所示。從圖2 中可以看出，葉綠素a 反演模型的擬合效果較好，擬合度為0.874，均方根誤差為2.691。

圖2 葉綠素a 參數(shù)的SVR 模型擬合圖

3.3 模型精度檢驗

為了驗證模型的準確性，本文利用SVR 反演模型對測試樣本進行葉綠素a 濃度估算，采用相對誤差指標對得到的估算值與實測值進行對比分析，從而對模型的精度進行檢驗（表3）。從表3 可以看出，模型整體精度較好，4 個測試點的相對誤差均控制在30%以內(nèi)，平均相對誤差為13.32%，可認為滿足精度要求。

表3 模型精度檢驗結(jié)果

4 結(jié)果與分析

為研究古田水庫葉綠素a 濃度的時空變化，本文分別選取2021 年4 個不同季節(jié)的月份的Sentinel-影像數(shù)據(jù)，利用已建立的葉綠素a 參數(shù)的SVR 反演模型對古田水庫進行遙感反演，并在ArcGIS 軟件平臺下進行制圖輸出，得到古田水庫葉綠素a 濃度的時空變化分布圖（圖3）。

圖3 古田水庫葉綠素a 濃度2021 年不同月份時空變化分布圖

從空間變化來看，古田水庫葉綠素a 濃度在不同季節(jié)的4 個月份存在著顯著的差異。3 月和6 月的葉綠素a 濃度整體偏高，庫區(qū)內(nèi)部葉綠素a 濃度略高于上游，呈現(xiàn)出南高北低的趨勢；9 月和12 月的葉綠素a 濃度整體偏低，較3 月和6 月呈現(xiàn)的是下降的趨勢，其中庫區(qū)內(nèi)葉綠素a 濃度下降明顯，使得空間分布上葉綠素a 濃度從南向北遞增，呈現(xiàn)上游明顯高于下游的特征。

從季節(jié)性變化來看，古田水庫葉綠素a 濃度春、夏2 季明顯高于秋、冬2 季，且夏季整體達到最高，主要是由于溫度對葉綠素a 的調(diào)節(jié)作用。古田水庫底部的氮磷營養(yǎng)物質(zhì)在上游城鎮(zhèn)生活污水、畜禽農(nóng)業(yè)污染和工業(yè)廢水的長期影響下不斷堆積，一旦天氣變暖，水溫升高，庫內(nèi)底泥中的氮、磷會被釋放出來，水庫中藻類大量繁殖，將導(dǎo)致庫內(nèi)葉綠素a 濃度上升。

5 結(jié)論

本文選用大寬幅的歐空局Sentinel-2 多光譜衛(wèi)星遙感數(shù)據(jù)，結(jié)合實測葉綠素a 參數(shù)濃度數(shù)據(jù)，基于SVR 模型構(gòu)建古田水庫葉綠素a 參數(shù)反演模型，并進行了精度檢驗。通過已建立的反演模型，對2021 年3、6、9、12 月等4 個月份的古田水庫葉綠素a 濃度進行遙感監(jiān)測，最終獲取了古田水庫4 個季節(jié)葉綠素a 濃度的時空變化分布特征，得出結(jié)論如下：

（1）古田水庫葉綠素a 濃度與波段組合間的相關(guān)性明顯高于單波段，葉綠素a 與波段組合（B5-B4）/B3 的相關(guān)性最高。

（2）采用測試樣本對葉綠素a 參數(shù)的SVR 反演模型進行檢驗，結(jié)果表明，葉綠素a 參數(shù)的SVR 模型反演效果較好，擬合度高達0.8 以上。SVR 反演模型可以有效地應(yīng)用于小樣本的葉綠素a 濃度遙感監(jiān)測。

（3）反演結(jié)果表明，春、夏2 季古田水庫葉綠素a 濃度整體偏高，主要是受溫度和降水的影響?？臻g分布上，3 月和6 月庫區(qū)內(nèi)部葉綠素a 濃度略高于上游，呈南高北低趨勢，9 月和12月從南向北遞增，呈現(xiàn)上游明顯高于下游的特征。葉綠素a 濃度問題主要與上游鄉(xiāng)鎮(zhèn)生活排污、農(nóng)業(yè)污染、畜禽養(yǎng)殖等人類活動密切相關(guān)。