曹 暢,王勝蕾,李俊生,趙紅莉,申 維,謝 婭

(1:中國地質大學(北京)地球科學與資源學院，北京 100083)(2:中國科學院空天信息創(chuàng)新研究院，數(shù)字地球重點實驗室，北京 100094)(3:北京大學遙感與地理信息研究所，北京 100871)(4:中國科學院大學電子電氣與通信工程學院,北京 100049)(5:中國水利水電科學研究院, 水資源研究所,北京 100048)

水體富營養(yǎng)化是水體中營養(yǎng)鹽增多而導致浮游植物爆發(fā)性生長的一種現(xiàn)象，富營養(yǎng)化的湖庫水體透明度和溶解氧濃度下降，容易暴發(fā)藍藻水華，對水環(huán)境質量造成嚴重的威脅[1-3]. 因此，動態(tài)監(jiān)測和評價湖庫水體的營養(yǎng)狀態(tài)具有非常重要的意義. 水體營養(yǎng)狀態(tài)的常規(guī)監(jiān)測方法是水面采集水樣并送到實驗室內測量葉綠素a、總氮、總磷等水質參數(shù)，進而計算營養(yǎng)狀態(tài)指數(shù)，然后根據(jù)營養(yǎng)狀態(tài)指數(shù)分級可以將水體營養(yǎng)狀態(tài)分為貧營養(yǎng)、中營養(yǎng)和富營養(yǎng)[4-6]. 這種方法雖然精度較高但費時、費力、費用高，而且只能獲得個別采樣點的結果. 與之相比，基于衛(wèi)星遙感監(jiān)測水體營養(yǎng)狀態(tài)具有快速、范圍廣、成本低等優(yōu)勢，適合開展大范圍動態(tài)監(jiān)測.

基于衛(wèi)星遙感監(jiān)測水體營養(yǎng)狀態(tài)主要是基于衛(wèi)星遙感反演的葉綠素a濃度計算水體營養(yǎng)指數(shù)，進而得到水體的營養(yǎng)狀態(tài)[7-8]. Thiemann等利用IRE-1C數(shù)據(jù)反演德國梅克倫堡州湖泊水體中的葉綠素a濃度，基于卡爾森模型評價其水體的富營養(yǎng)化狀態(tài)[9]；Matthews等通過測量內陸水體葉綠素a濃度，監(jiān)測非洲南部多個湖庫水體10年間的富營養(yǎng)化及其變化情況[10]. Gilerson等利用紅色和近紅外波段算法評價美國布拉斯加州的葉綠素a濃度，觀測該州湖泊水體的營養(yǎng)狀態(tài)[11]；朱利等利用GF-1號衛(wèi)星WFV數(shù)據(jù)對太湖葉綠素a濃度、透明度、懸浮物濃度進行遙感監(jiān)測，并評價水體富營養(yǎng)狀態(tài)和水質參數(shù)空間分布[12]；Le等在三波段算法的基礎上發(fā)展四波段解析算法提高三波段算法的性能，比較兩種算法在反演太湖葉綠素a濃度的精度，并對太湖營養(yǎng)狀態(tài)進行評價[13].

然而，由于內陸水體光學特性復雜，受浮游植物、懸浮泥沙、有色可溶性有機物共同影響，而且隨區(qū)域和季節(jié)變化大，因為不同季節(jié)的水體溫度變化對于藻類生長的影響、地表徑流和風速變化對于懸浮物的影響往往存在差異，進而導致很多湖庫水體葉綠素a反演算法都具有很強的區(qū)域和季節(jié)局限性[14-16]，因此限制了基于衛(wèi)星遙感的大范圍水體營養(yǎng)狀態(tài)監(jiān)測. 與傳統(tǒng)的基于葉綠素a反演結果的營養(yǎng)狀態(tài)衛(wèi)星遙感監(jiān)測不同，Wang等[17]發(fā)展了基于水色指數(shù)的水體營養(yǎng)狀態(tài)等級評價方法，并且成功應用于2012年全球大型湖庫水體營養(yǎng)狀態(tài)評價，利用文獻和公報中的100個水體實地監(jiān)測結果進行評價得到精度為80%，這為開展大范圍湖庫水體營養(yǎng)狀態(tài)評價奠定了基礎.

隨著工農業(yè)的快速發(fā)展，我國湖庫水體富營養(yǎng)化的趨勢較為嚴重，開展全國范圍內重點湖庫營養(yǎng)狀態(tài)的動態(tài)監(jiān)測具有重要意義，有利于宏觀掌握全國湖庫營養(yǎng)狀態(tài)的時空分布，為管理決策提供數(shù)據(jù)支持. 但是，《2018年中國生態(tài)環(huán)境狀態(tài)公報》只監(jiān)測了111個湖庫的營養(yǎng)狀態(tài)，前幾年監(jiān)測的湖庫數(shù)量更少. 因此，非常有必要開展基于衛(wèi)星遙感的全國范圍內湖庫的營養(yǎng)狀態(tài)監(jiān)測，作為公報的補充，一方面可以從空間上監(jiān)測更多的湖庫，另一方面可以從時間上追溯湖庫營養(yǎng)狀態(tài)的變化. 因此本文將利用基于水色指數(shù)的營養(yǎng)狀態(tài)評價方法[17]，生產全國范圍內重點湖庫的營養(yǎng)狀態(tài)遙感監(jiān)測產品，進而分析全國重點湖庫營養(yǎng)狀態(tài)空間分布規(guī)律及其主要影響因素.

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

1.1 數(shù)據(jù)源

搭載在Terra和Aqua衛(wèi)星的中分辨率成像光譜儀(MODIS)具有較高的時間分辨率，一天內可以覆蓋全球兩次(部分赤道低緯度地區(qū)除外). MOD09數(shù)據(jù)是MODIS(Terra)的地表反射率產品，包含MODIS的1～7共7個波段. MOD09A1數(shù)據(jù)是八天合成的MOD09數(shù)據(jù)，降低了部分有云等質量不好數(shù)據(jù)的影響而且將全球范圍劃分為460個規(guī)則的網(wǎng)格分塊，便于生產大范圍產品并進行時序產品疊加分析，因而在業(yè)務化的大范圍水質監(jiān)測中具有重要優(yōu)勢. 因此本文選擇MOD09A1數(shù)據(jù)作為全國重點湖庫營養(yǎng)狀態(tài)業(yè)務化監(jiān)測的遙感數(shù)據(jù)源，利用覆蓋全國重點湖庫的18塊MOD09A1產品，網(wǎng)格分塊行列號包括：h23v04、h23v05、h24v04、h24v05、h25v03、h25v04、h25v05、h25v06、h26v04、h26v05、h26v06、h27v04、h27v05、h27v06、h28v05、h28v06、h28v07、h29v06.

1.2 研究區(qū)

MOD09A1數(shù)據(jù)的空間分辨率是500 m，為了降低岸邊混合像元或臨近像元的影響，本文主要研究面積大于5 km2的湖庫水體，也就是水體的聯(lián)通像元數(shù)量大于20個.

在全國面積大于5 km2的湖庫中選擇了144個主要湖庫，包括111個湖泊和33個水庫. 111個湖泊中包括洱海等5個飲用水源地，也包括巢湖等較渾濁和撫仙湖等較清潔的106個非飲用水源地湖泊；33個水庫包括密云水庫等19個飲用水源地，以及于橋水庫等較渾濁和新安江水庫等較清潔的14個非飲用水源地水庫. 這些湖庫的空間分布覆蓋了中國的五大湖區(qū)，其中20個位于蒙新高原湖區(qū)，34個位于青藏高原湖區(qū)，7個位于云貴高原湖區(qū)，66個位于東部平原湖區(qū)，17個位于東北山地與平原湖區(qū).

2 研究方法

本文首先對MOD09A1地表反射率產品進行離水反射率校正，然后提取湖庫水體分布范圍，進而計算水色指數(shù)(FUI)，并利用FUI生產水體營養(yǎng)狀態(tài)等級，最后分析全國重點湖庫營養(yǎng)狀態(tài)的空間分布規(guī)律.

2.1 離水反射率校正

MOD09A1是地表反射率產品，已對其7個波段進行了氣溶膠校正、卷積云校正和瑞利散射校正[18]，但是沒有對天空光及太陽耀斑的水面反射進行校正. 本文采用基于近紅外波段(NIR)和短波紅外波段(SWIR)的暗目標法進行校正[19]，將MOD09A1的地表反射率(R(λ))減去這個像元的近紅外和短波紅外的最小值，再除以π，轉換成離水反射率(Rrs(λ)).

2.2 湖庫水體提取

在對MOD09A1數(shù)據(jù)進行離水反射率校正的基礎上，還要對湖庫的水體邊界進行提取. 水體在短波紅外波段具有強吸收作用，因而水體像元在短波紅外波段反射率很低，相反陸地像元的反射率較高[20]，因此本文采用基于短波紅外波段灰度直方圖的雙峰谷值閾值自動化確定的方法提取水體邊界[21]. 為了降低接近陸地區(qū)域混合像元和光學淺水的影響，將得到的水體部分向內腐蝕一個像元. 最后，去除聯(lián)通像元數(shù)小于20的水體，只保留面積在5 km2以上的水體.

2.3 水色指數(shù)計算

利用Wang等[22]提出的從太湖MOD09數(shù)據(jù)產品中提取FUI指數(shù)的方法對全國主要湖庫的FUI指數(shù)進行求解.FUI是一種比色表，將水體顏色從深藍色到黃褐色劃分為21個級別. 利用MODIS的第1(紅)、3(藍)、4(綠)波段計算CIE顏色空間中的三刺激值[x]、[y]、[z][23]，并進一步計算CIE顏色空間色度圖中找到顏色對應的坐標(x,y)，進而找到色度角α，在FUI指數(shù)查找表中找到與色度角最接近的色度值，此值對應的FUI值就是水體的FUI值[22].

2.4 基于FUI的營養(yǎng)狀態(tài)等級評價

采用基于FUI和紅波段閾值法對湖庫營養(yǎng)狀態(tài)進行評價[22]. 貧營養(yǎng)水體一般比較清潔，藻類濃度比較低，水體一般呈現(xiàn)藍色；富營養(yǎng)水體較為渾濁，藻類濃度較高，水體顏色多為黃綠色；中營養(yǎng)水體的藻類濃度以及水體顏色一般介于貧營養(yǎng)和富營養(yǎng)之間.FUI可以表示水體顏色，因此可以使用FUI的閾值分割來評價水體的營養(yǎng)狀態(tài)；當FUI<7時，水體為貧營養(yǎng)；當7≤FUI<10時，水體為中營養(yǎng)；當FUI≥10時，水體為富營養(yǎng). 但是，也有一些中營養(yǎng)水體的FUI≥10，使用上述分割模型會被誤判為富營養(yǎng)水體. 通過分析發(fā)現(xiàn)這些誤判水體的紅波段離水反射率(Rrs(645))一般較小，因此增加一個閾值分割限制條件進行二次限定：當Rrs(645)<0.00625時，即使FUI≥10，水體仍然是中營養(yǎng)[22].

2.5 精度評價

可以利用與MODIS衛(wèi)星同步測量的水面實測數(shù)據(jù)計算水體營養(yǎng)狀態(tài)等級，對MODIS監(jiān)測結果進行精度評價. 于2018到2019年在太湖和于橋水庫分別獲取了12個和46個采樣點的葉綠素a濃度實測數(shù)據(jù)，利用這些葉綠素a濃度數(shù)據(jù)可以計算綜合營養(yǎng)狀態(tài)指數(shù)(TLI)[2]. 由于這些點大部分屬于富營養(yǎng)，為了增加貧營養(yǎng)水體的檢驗數(shù)據(jù)，從文獻[24]中獲取了2014-2016年青藏高原納木錯、塔若錯、色林錯、多爾索洞錯、赤布張錯5個湖泊的透明度數(shù)據(jù)，利用營養(yǎng)狀態(tài)指數(shù)計算公式：TLI(Chl.a)=10(2.5+1.086 ln Chl.a)、TLI(TP)=10(9.436+1.624 ln TP)、TLI(TN)=10(5.453+1.694 ln TN)、TLI(SD)=10(5.118-1.94 ln SD)、TLI(COD)=10(0.109+ 2.661 ln COD)(式中, 葉綠素a(Chl.a)單位為mg/m3；透明度(SD)單位為m; 其他指標單位均為mg/L)可以計算綜合營養(yǎng)狀態(tài)指數(shù)，通過TLI數(shù)值大小和分級標準對湖庫水體進行營養(yǎng)狀態(tài)分級[2]. 利用這63個采樣點實測數(shù)據(jù)計算得到的營養(yǎng)狀態(tài)等級與同步MODIS衛(wèi)星遙感監(jiān)測的營養(yǎng)狀態(tài)等級進行對比，計算得到基于MODIS監(jiān)測營養(yǎng)狀態(tài)等級的總體精度為88.9%. 其中，貧營養(yǎng)的判別精度比較高，少量誤判主要發(fā)生在部分富營養(yǎng)和中營養(yǎng)之間. 總體上來說，營養(yǎng)狀態(tài)評價的精度基本上可以滿足大范圍應用研究需求. Guan等[25]在2018年實地測量的葉綠素a濃度數(shù)據(jù)中有26個湖泊與本文研究區(qū)域一致，采用綜合營養(yǎng)狀態(tài)指數(shù)法計算26個湖泊的營養(yǎng)狀態(tài)，有25個湖泊與本文研究結果相同，MODIS的監(jiān)測精度達到96.0%，鄒偉等[26]在2018年7-8月份采集了長江中下游地區(qū)19個湖庫水體的葉綠素a濃度，其中與本文研究區(qū)重疊的湖庫共有13個，基于葉綠素a濃度計算的營養(yǎng)狀態(tài)等級為12個湖庫為富營養(yǎng)化，1個湖庫為中營養(yǎng). 這13個湖庫中的12個的營養(yǎng)狀態(tài)調查結果與本文基于MODIS監(jiān)測的結果一致，精度為92.3%. Song等[27]對2014-2018年夏季全國面積大于80000 m2的湖庫進行大范圍高精度水體透明度反演，利用綜合營養(yǎng)狀態(tài)指數(shù)計算得到東北山地與平原湖區(qū)和東部平原湖區(qū)水體較多呈現(xiàn)富營養(yǎng)狀態(tài)，對比圖1，與本文研究的結果一致.

2.6 全國重點湖庫營養(yǎng)狀態(tài)產品生產及空間分析

基于上述技術流程，編寫IDL程序，實現(xiàn)基于MOD09A1數(shù)據(jù)的全國重點湖庫營養(yǎng)狀態(tài)產品自動化批量生產. 由于水色指數(shù)的營養(yǎng)狀態(tài)評價方法主要是基于藻類對水體顏色的影響，因此選擇藻類生長情況最好的夏季來分析全國重點湖庫水體營養(yǎng)狀態(tài)的空間分布情況，進而分析其主要影響因素.

3 結果與討論

3.1 全國重點湖庫營養(yǎng)狀態(tài)空間分布格局

基于2018年6-8月之間的MOD09A1數(shù)據(jù)，生產了全國144個重點湖庫的營養(yǎng)狀態(tài)夏季均值產品. 為了便于統(tǒng)計分析，將每個湖庫所有像元的營養(yǎng)狀態(tài)求一個均值，得到每個湖庫平均的營養(yǎng)狀態(tài)等級，其中貧營養(yǎng)、中營養(yǎng)、富營養(yǎng)的湖庫數(shù)量分別為23、34、87個，144個重點湖庫營養(yǎng)狀態(tài)空間分布如圖1所示，可以看出我國湖庫水體的營養(yǎng)狀態(tài)分布不均勻，東部湖庫水體以中到富營養(yǎng)狀態(tài)為主，尤其是長江中下游和東北山地與平原湖區(qū)湖庫水體富營養(yǎng)化比例非常高，西部湖庫水體以貧到中營養(yǎng)狀態(tài)為主，尤其是青藏高原湖區(qū)貧營養(yǎng)比例很高. 不同水體類型營養(yǎng)狀態(tài)等級比例如表1，其中湖泊型飲用水源地的富營養(yǎng)化比例要高于水庫型飲用水源地.

表1 基于MODIS監(jiān)測的2018年夏季全國重點湖庫中不同類型水體的營養(yǎng)狀態(tài)等級比例

五大湖區(qū)中各營養(yǎng)狀態(tài)等級所占的比例如圖2所示，東北山地與平原湖區(qū)水體富營養(yǎng)狀態(tài)比例最大，東部平原湖區(qū)次之，2個湖區(qū)均無貧營養(yǎng)狀態(tài)水體. 中營養(yǎng)狀態(tài)水體在蒙新高原湖區(qū)數(shù)量最多. 青藏高原湖區(qū)貧營養(yǎng)狀態(tài)水體數(shù)量最多，其次是云貴高原湖區(qū). 在五大湖區(qū)中各選1個典型湖泊，營養(yǎng)狀態(tài)分布如圖3所示. 其中，興凱湖(東北山地與平原湖區(qū)，中俄跨界湖)、巢湖(東部平原湖區(qū))水體富營養(yǎng)化程度較高；博斯騰湖(蒙新高原湖區(qū))整體上為中營養(yǎng)狀態(tài)；色林錯(青藏高原湖區(qū))、撫仙湖(云貴高原湖區(qū))整體上為貧營養(yǎng)，水質情況較好.

圖1 基于MODIS監(jiān)測的2018年夏季全國重點湖庫水體營養(yǎng)狀態(tài)等級分布

自然環(huán)境是影響水體營養(yǎng)狀態(tài)的重要因素. 青藏高原湖區(qū)的海拔較高, 多在3500 m以上，氣溫較低[28-29]，水中藻類進行光合作用的速率受低水溫限制[30]，藻類數(shù)量一般較少，水體多為貧營養(yǎng)狀態(tài)[31]. 云貴高原湖區(qū)深水湖和淺水湖交錯分布，該湖區(qū)水體營養(yǎng)狀態(tài)差異較大[32]. 東北地區(qū)的中部為平原, 周圍三面環(huán)山，位于平原地區(qū)的湖泊面積較小，且湖水較淺[33]，降水多集中在夏季，會帶入營養(yǎng)鹽等物質，造成富營養(yǎng)狀態(tài)水體數(shù)量較多. 東部平原湖區(qū)人口密度高，面源污染帶來大量的營養(yǎng)鹽；同時長江中下游地區(qū)主要湖泊的水深多在10 m以下[34-35]，湖泊底泥在環(huán)境條件適宜時與湖水進行交互，也成為營養(yǎng)鹽的主要來源之一[36]. 蒙新高原湖區(qū)降水少，蒸發(fā)強，該湖區(qū)草場退化且土地沙漠化，造成地表徑流減少，再加上水土流失導致入湖水體中營養(yǎng)鹽濃度上升，水體營養(yǎng)狀態(tài)較高[37-38].

人類活動也是影響水體營養(yǎng)狀態(tài)的重要因素，湖泊的富營化水平與社會經(jīng)濟發(fā)展進程密切相關[39]. 東部平原湖區(qū)和東北山地湖區(qū)人口眾多，加上經(jīng)濟發(fā)展速度快，工業(yè)點源污染以及農業(yè)非點源污染都會為周邊湖庫帶來營養(yǎng)物質[40-41]，導致湖庫水體的營養(yǎng)狀態(tài)整體偏高. 青藏高原湖區(qū)海拔較高，自然環(huán)境較惡劣，人口密度小，城市化進程緩慢，水體受人為干擾因素較小，水體多為貧營養(yǎng)狀態(tài).

3.2 營養(yǎng)狀態(tài)等級監(jiān)測結果與環(huán)境公報的對比

基于FUI的營養(yǎng)狀態(tài)評價方法是利用影響水體顏色的葉綠素a和透明度信息，因此應該使用藻類生長最充分的夏季，才能對比出來不同水體的差異. 而且很多北方和青藏高原湖庫在冬季結冰，也無法使用. 綜上所述，基于FUI的營養(yǎng)狀態(tài)評價方法主要適用于夏季. 《2018中國生態(tài)環(huán)境狀態(tài)公報》監(jiān)測了111個湖庫2018年平均的營養(yǎng)狀態(tài)，其中有40個湖庫與本文研究區(qū)重疊. 經(jīng)過對比發(fā)現(xiàn)，其中，有24個湖庫(占60%)的營養(yǎng)狀態(tài)與本文監(jiān)測結果一致，另外16個湖庫(占40%)的營養(yǎng)狀態(tài)不一致，主要是公報評價為中營養(yǎng)，而本文評價為富營養(yǎng). 公報方法與本文方法不同，主要體現(xiàn)在：(1)公報是基于葉綠素a、透明度、總氮、總磷、化學需氧量5個參數(shù)計算綜合營養(yǎng)狀態(tài)指數(shù)，而本文是基于水體顏色信息評價水體營養(yǎng)狀態(tài)，其中主要是利用了對水體顏色影響較大的葉綠素a和透明度信息；(2)公報一般是基于每月一次測量的少數(shù)幾個采樣點的數(shù)據(jù)的評價結果，而本文是基于時間序列的全湖遙感圖像平均的結果. 本文的營養(yǎng)狀態(tài)評價方法可以作為公報的補充，一方面可以用于監(jiān)測公報沒有監(jiān)測的湖庫，另一方面對于公報監(jiān)測的湖庫，還可以補充其空間分布和時間變化信息.

圖2 基于MODIS監(jiān)測的2018年夏季全國五大湖區(qū)湖庫營養(yǎng)狀態(tài)等級比例

圖3 基于MODIS監(jiān)測的2018年夏季全國五大湖區(qū)典型湖泊水體營養(yǎng)狀態(tài)

4 結論

本文以MODIS八天合成的地表反射率產品MOD09A1為遙感數(shù)據(jù)源，利用基于水色指數(shù)的營養(yǎng)狀態(tài)分級方法，監(jiān)測全國144個重點湖庫營養(yǎng)狀態(tài)等級. 利用63個水面實測數(shù)據(jù)計算得到的營養(yǎng)狀態(tài)等級與同步MODIS衛(wèi)星遙感監(jiān)測結果進行對比，得到MODIS的監(jiān)測精度為88.9%，精度基本上可以滿足大范圍應用研究需求.

利用經(jīng)過檢驗的方法，基于MODIS數(shù)據(jù)生產了2018年夏季全國144個重點湖庫的營養(yǎng)狀態(tài)等級，分析了全國重點湖庫營養(yǎng)狀態(tài)的空間分布情況(圖2). 貧營養(yǎng)、中營養(yǎng)、富營養(yǎng)的湖庫比例分別為16%、24%、60%，富營養(yǎng)化的比例很高. 營養(yǎng)狀態(tài)空間上分布不均勻，東部地區(qū)湖庫(主要是東部平原和東北山地與平原湖區(qū))多為富營養(yǎng)狀態(tài)，東北山地與平原湖區(qū)的富營養(yǎng)化湖庫比例最多達到94%，東部平原湖區(qū)的富營養(yǎng)化湖庫比例達到80%. 西部湖區(qū)湖庫(特別是青藏高原湖區(qū))多為貧營養(yǎng)狀態(tài)，青藏高原湖區(qū)貧營養(yǎng)湖庫比例達到62%，富營養(yǎng)化水體僅占9%. 海拔和地表溫度等自然因素與工業(yè)點源和農業(yè)面源污染等人為因素是湖庫營養(yǎng)狀態(tài)的重要影響因素.