基于BP神經(jīng)網(wǎng)絡(luò)算法的密云水庫水質(zhì)參數(shù)反演研究

馬豐魁，姜群鷗,2*，徐藜丹，梁勇，王榮臣，蘇帥

1.北京林業(yè)大學(xué)水土保持學(xué)院，100083 北京；2.北京林業(yè)大學(xué)/水土保持與荒漠化防治教育部重點實驗室，100083 北京；3.北京市密云水庫管理處，101512 北京

水是人類生存發(fā)展必需的自然資源，是社會經(jīng)濟發(fā)展的重要物質(zhì)保障，是生態(tài)環(huán)境的控制性因素。然而隨著中國經(jīng)濟社會的發(fā)展，水資源需求量日益增加，水污染整體呈現(xiàn)惡化態(tài)勢。水資源短缺和水污染已經(jīng)成為制約中國可持續(xù)發(fā)展的瓶頸因素（曾昭等，2013）。北京是一座重度資源型缺水、水污染嚴重的城市，經(jīng)濟快速發(fā)展，城市擴張以及人口持續(xù)增長導(dǎo)致北京市的水質(zhì)不斷惡化（高陽等，2018）。在這種背景下，密云水庫作為北京市最大的飲用水源供應(yīng)地，其水質(zhì)安全不僅關(guān)乎北京市兩千多萬市民的生命健康和生活用水問題，而且對于北京市發(fā)展和首都各項職能的正常運轉(zhuǎn)具有重大影響，引起了國家和政府的高度重視。保障城市水源安全、提高城鎮(zhèn)飲用水源的安全保障水平，是城市生態(tài)文明建設(shè)的重要內(nèi)容，也是生態(tài)文明建設(shè)的重要組成和基礎(chǔ)保障（馬巍等，2016）。因而探究監(jiān)測密云水庫主要水質(zhì)參數(shù)變化的有效方法，實現(xiàn)對密云水庫主要水質(zhì)指標的快速反演尤為重要。

傳統(tǒng)的水質(zhì)監(jiān)測一般以實地采樣分析為主，可以較為準確地測得多種水質(zhì)參數(shù)，但是需要耗費大量的時間和成本，且只能獲得點狀數(shù)據(jù)，所得數(shù)據(jù)在時間和空間上都不連續(xù)，無法滿足大范圍、實時的水質(zhì)監(jiān)測要求。（趙玉芹等，2009）。水質(zhì)遙感技術(shù)作為一種具有區(qū)域化監(jiān)測能力的方法，可以快速、實時的獲取水質(zhì)的時空分布，在水質(zhì)監(jiān)測中日益受到重視（曹引等，2019）。常規(guī)的遙感水質(zhì)參數(shù)反演方法主要有經(jīng)驗方法、半經(jīng)驗法和半分析法（張海威等，2017）。翟召坤等（2018）基于高分一號衛(wèi)星數(shù)據(jù)和實測數(shù)據(jù)，采用典型的經(jīng)驗方法對潘家口水庫水質(zhì)參數(shù)和敏感波段進行了線性擬合，得到水質(zhì)參數(shù)遙感反演模型。林劍遠等（2019）基于航空和水表高光譜遙感數(shù)據(jù)，利用半經(jīng)驗法對浙江省嘉興市主城區(qū)河網(wǎng)水質(zhì)進行了監(jiān)測。石亮亮（2019）基于實測數(shù)據(jù)，開發(fā)了水體IOPs及CDOM反演的半分析算法，對東海、千島湖水域進行了反演。這些方法本質(zhì)上都是通過線性回歸來實現(xiàn)水質(zhì)反演，而水環(huán)境系統(tǒng)具有強烈的非線性和不確定性特征，傳統(tǒng)的線性回歸不能完全反映其變化規(guī)律，通過線性回歸分析會導(dǎo)致結(jié)果不準確（Chen et al.，2018）。

近年來，人工神經(jīng)網(wǎng)絡(luò)理論的迅速發(fā)展為水質(zhì)評價和預(yù)測提供了一種新的有效方法。人工神經(jīng)網(wǎng)絡(luò)具有獨特的分布并行處理、非線性映射和自適應(yīng)學(xué)習(xí)等優(yōu)點，BP神經(jīng)網(wǎng)絡(luò)是其中一種較為成熟的非線性函數(shù)逼近的方法，所得結(jié)果更為準確（李輝東等，2015；黃俊等，2015；易湘生等，2012）。國內(nèi)外許多專家學(xué)者運用神經(jīng)網(wǎng)絡(luò)進行了水質(zhì)參數(shù)的研究。Gazzaz et al.（2012）通過建立神經(jīng)網(wǎng)絡(luò)模型完成了對馬來西亞金塔河的水質(zhì)參數(shù)的預(yù)測，研究結(jié)果表明，該模型預(yù)測結(jié)果可以解釋95.4%左右的實測水質(zhì)參數(shù)的變化。Wu et al.（2009）基于MODIS數(shù)據(jù)集，采用BP神經(jīng)網(wǎng)絡(luò)、多元線性回歸等方法對巢湖的水質(zhì)參數(shù)進行了反演，結(jié)果表明，BP神經(jīng)網(wǎng)絡(luò)模型優(yōu)于多元線性回歸，相對誤差低于35%。Song et al.（2011）使用Landsat TM數(shù)據(jù)和實測數(shù)據(jù)，建立了經(jīng)驗回歸和神經(jīng)網(wǎng)絡(luò)，對查干湖的水質(zhì)參數(shù)進行反演，結(jié)果表明神經(jīng)網(wǎng)絡(luò)的均方根誤差小于10%，而回歸的均方根誤差小于25%。Xiao et al.（2017）以漢江中下游為研究區(qū)，建立了VIP-BP模型估計葉綠素含量，結(jié)果表明，VIP-BP模型可以準確真實的反映葉綠素含量變化，且誤差明顯低于3波段模型。

上述研究表明，神經(jīng)網(wǎng)絡(luò)與遙感技術(shù)的結(jié)合可以很好地實現(xiàn)水質(zhì)參數(shù)的反演。在此背景下，本研究基于遙感和GIS技術(shù)，采用BP神經(jīng)網(wǎng)絡(luò)算法，結(jié)合地面監(jiān)測數(shù)據(jù)和Landsat 8遙感影像，分別建立了總磷、總氮、氨氮、COD（化學(xué)需氧量）4個主要水質(zhì)參數(shù)的神經(jīng)網(wǎng)絡(luò)模型，實現(xiàn)了密云水庫主要水質(zhì)參數(shù)的反演。旨在為密云水庫水質(zhì)監(jiān)測探索一種便捷、可靠、高效的方法，從而對密云水庫主要水質(zhì)參數(shù)變化做出科學(xué)評估，為密云水庫水質(zhì)監(jiān)測和保護工作提出重要的科學(xué)建議。

1 材料與方法

1.1 研究區(qū)概況

圖1 密云水庫區(qū)位圖Fig.1 Location of Miyun Reservoir Area

密云水庫位于北京市東北部密云縣城北13 km處（116°07′—117°30′E，40°14′—41°05′N），地處燕山群山之中（圖1）。庫區(qū)建成于1960年9月，流域面積15788 km2，最大水面面積為188 km2，最大庫容43.75 ×108m3，最大水深43.5 m，是華北地區(qū)最大的水庫。密云水庫有2大入庫河流，分別是白河和潮河，分為白河、潮河和內(nèi)湖3個庫區(qū)（張敏等，2019）。該區(qū)屬于暖溫帶半濕潤季風(fēng)氣候，年均溫 10.5 ℃，年均最低氣溫-18 ℃，最高氣溫38 ℃，降水主要集中在6—8月（李東青等，2015）。庫區(qū)主要功能是防洪、灌溉、發(fā)電等，為北京市、天津市和河北省服務(wù)。隨著北京城市化和經(jīng)濟的快速發(fā)展，水資源短缺的局面日益加重，密云水庫于1981年開始專為北京市供水，且生活用水比重逐年上升，當(dāng)前已成為北京市最大的飲用水源供應(yīng)地，被譽為“北京生命之水”。自20世紀90年代以來，水體的富營養(yǎng)化已成為密云水庫最主要的水質(zhì)污染問題（許爾琪等，2018）。因此開展對密云水庫的水質(zhì)參數(shù)的反演研究十分必要。

1.2 數(shù)據(jù)來源

1.2.1 地面監(jiān)測數(shù)據(jù)

研究團隊于2018年10月17日乘船進入密云水庫，在規(guī)劃采樣區(qū)域利用取水器距離水面30 cm深處共采集了80個樣品，同時記錄了采樣點的GPS坐標。采樣點的布設(shè)原則是盡可能均勻分布，敏感區(qū)域多布設(shè)一些取樣點。采集回樣品之后，當(dāng)天將樣品送檢相關(guān)部門，在實驗室對每個樣品進行測試：水質(zhì)總氮采用化學(xué)發(fā)光器測量基態(tài)NO2，將結(jié)果換算為總氮；水質(zhì)氨態(tài)氮采用SEALAA3連續(xù)流動分析儀測試；水質(zhì)總磷采用鉬酸銨分光光度法測定；水質(zhì)化學(xué)需氧量采用快速消解分光光度法測定。（陳穎等，2018；廖蕾等，2016；GB/T 11893—1989；HJ/T 399—2007）最終得到80采樣點的總磷、氨氮、總磷、COD（化學(xué)需氧量）4個指標的質(zhì)量濃度。

1.2.2 Landsat 8遙感影像

本研究采用的遙感數(shù)據(jù)為Landsat 8，遙感影像下載自中國科學(xué)院遙感與數(shù)字地球研究所（http://eds.ceode.ac.cn/）和地理空間數(shù)據(jù)云（http://www.gscloud.cn/sources），空間分辨率為 30 m×30 m。選取研究區(qū)的影像時間跨度為 2013—2018年，其中用于構(gòu)建神經(jīng)網(wǎng)絡(luò)模型的遙感影像時間為2018年10月17日，與水質(zhì)樣品采集時間一致。由于密云水庫從12月至次年2月有冰期，該時段遙感影像的反射率無法真實反映水質(zhì)參數(shù)的變化，再加上一些時段的遙感影像受到云層的影響，最終選取2013—2018年非冰期共36期云量低于30%的遙感影像，代表非冰期主要水質(zhì)指標的變化。原始影像含5個數(shù)據(jù)集：多光譜數(shù)據(jù)集（1—7波段）、全色波段數(shù)據(jù)（8波段）、卷云波段數(shù)據(jù)（9波段）、熱紅外數(shù)據(jù)（10、11波段）和質(zhì)量波段數(shù)據(jù)（12波段）以及含有整個圖像內(nèi)容的MTL文件，共13個文件。本研究選取多光譜數(shù)據(jù)（1—7波段）作為反演密云水庫水質(zhì)的原始數(shù)據(jù)。

1.3 數(shù)據(jù)處理

1.3.1 Landsat 8遙感影像預(yù)處理

本研究首先對Landsat 8影像進行輻射定標和大氣校正，使遙感影像表達真實的反射率（秦雁等，2011）。處理完成后將Landsat 8多光譜數(shù)據(jù)的7個波段分別導(dǎo)出為tif圖像，便于后期建模和預(yù)測。再通過歸一化水體指數(shù)NDWI，對研究區(qū)水域邊界進行提取。然后使用提取出的密云水庫的水域邊界對遙感影像進行裁剪，得到研究區(qū)各波段的遙感影像。最后將遙感影像進行轉(zhuǎn)換，得到研究區(qū)全區(qū)的反射率值。

1.3.2 采樣點反射率提取

本研究將記錄的采樣點的 GPS坐標離散到空間中，并將采樣點與Landsat 8影像疊加，用ArcGIS提取Landsat 8各波段在采樣點的反射率值，將提取的采樣點反射率數(shù)據(jù)進行整理，作為建立神經(jīng)網(wǎng)絡(luò)模型的輸入數(shù)據(jù)。

1.3.3 敏感因子提取

水體的光譜特征由其中的各種物質(zhì)對光輻射的吸收和散射性質(zhì)所決定，當(dāng)水體受到污染時，光譜特征會發(fā)生變化，所以遙感影像各波段與水質(zhì)參數(shù)之間存在一定的相關(guān)性。使用 SPSS軟件，對Landsat 8的1—7波段的反射率和4個水質(zhì)參數(shù)（總氮、氨氮、總磷、COD）的實測值進行皮爾遜相關(guān)性分析，以相關(guān)性分析為輔助，在保留最大信息量的原則下剔除干擾波段，得到水質(zhì)參數(shù)的敏感波段組合（表1），使得BP神經(jīng)網(wǎng)絡(luò)的擬合效果最優(yōu)。

1.4 研究方法

1.4.1 BP神經(jīng)網(wǎng)絡(luò)

BP（Back Propagation）神經(jīng)網(wǎng)絡(luò)是按誤差反向傳播算法訓(xùn)練的多層前饋式網(wǎng)絡(luò)。由Rumelhart和Mccelland提出，其運算原理為使用最小二乘法使網(wǎng)絡(luò)的誤差平方和達到最小。網(wǎng)絡(luò)在訓(xùn)練過程中利用錯誤反饋傳播不斷調(diào)整閾值和權(quán)值，使輸出結(jié)果達到最佳值，實現(xiàn)最佳的擬合效果（盛夏等，2013）。典型的BP神經(jīng)網(wǎng)絡(luò)拓撲結(jié)構(gòu)包括輸入層、隱含層和輸出層，每一層包含若干神經(jīng)元，上下層神經(jīng)元之間全連接，而同層神經(jīng)元之間無連接，輸入層與隱含層神經(jīng)元的連接值是網(wǎng)絡(luò)的權(quán)值，學(xué)習(xí)過程的目標是找到最優(yōu)的權(quán)值集，該權(quán)值集可以產(chǎn)生正確的輸出（李麗等，2008；張斌等，2012）。模型的基本建立過程是：給定神經(jīng)網(wǎng)絡(luò)輸入和期望值，由輸入層傳遞到隱含層，經(jīng)隱含層處理后，再傳遞到輸出層，由輸出層處理產(chǎn)生結(jié)果，此為信息前向傳播；計算實際輸出與期望輸出之間的誤差，將誤差值沿網(wǎng)絡(luò)反向傳播，并修正連接權(quán)值，此為誤差反向傳播；給定另一個輸入，重復(fù)上述過程，直到網(wǎng)絡(luò)輸出的誤差達到允許的范圍或達到設(shè)定的訓(xùn)練次數(shù)為止，得到研究所需的神經(jīng)網(wǎng)絡(luò)模型（張青等，2016；劉淑梅等，2015；劉曉芬等，2011）。

表1 水質(zhì)參數(shù)敏感波段Table 1 Water quality parameter sensitive band

1.4.2 納什系數(shù)

納什效率系數(shù)（Nash-Sutcliffeefficiency coefficient），用于驗證水文模型模擬結(jié)果的好壞，納什效率系數(shù)（E）的計算公式為：

式中：Qo指觀測值；Qm指模擬值；Qt表示第t時刻的某個值；Qo表示觀測值的平均值。E取值為負無窮至1，表示模式質(zhì)量好，模型可信度高；E接近 0，表示模擬結(jié)果接近觀測值的平均值水平，即總體結(jié)果可信，但過程模擬誤差大；E遠遠小于0，則模型不可信。

2 結(jié)果與分析

2.1 BP神經(jīng)網(wǎng)絡(luò)模型的構(gòu)建

本研究采用的是3層結(jié)構(gòu)的BP神經(jīng)網(wǎng)絡(luò)，在MATLAB中使用feedforwardnet函數(shù)，默認的訓(xùn)練函數(shù)為 trainlm（Levenberg-Marquardt算法），選擇傳遞函數(shù)tansig，輸出函數(shù)pureline，隱含層節(jié)點數(shù)通過試選法確定。以2018年10月17日的Landsat 8遙感影像提取的 80個采樣點敏感波段反射率為輸入數(shù)據(jù)，以同天對應(yīng)的地面監(jiān)測水質(zhì)參數(shù)作為目標數(shù)據(jù)，構(gòu)建一個3層的BP神經(jīng)網(wǎng)絡(luò)，包括一個輸入層、隱含層和輸出層。feedforwardnet函數(shù)將輸入數(shù)據(jù)分為3份，70%用于訓(xùn)練，15%用于檢驗，15%用于測試，防止過度擬合。feedforwardnet默認隱含層節(jié)點數(shù)為 10，學(xué)習(xí)步長為 0.05，迭代次數(shù)為50000，學(xué)習(xí)誤差為 1e-10。通過調(diào)整隱含層節(jié)點數(shù)、學(xué)習(xí)步長和迭代次數(shù)，對模型進行反復(fù)訓(xùn)練，直到BP神經(jīng)網(wǎng)絡(luò)模型的r值達到最優(yōu)（表2），即r值不再上升為止。得到4個參數(shù)的最優(yōu)神經(jīng)網(wǎng)絡(luò)模型，網(wǎng)絡(luò)結(jié)構(gòu)見（表3）。

表2 BP神經(jīng)網(wǎng)絡(luò)模型相關(guān)系數(shù)rTable 2 The correlation coefficient of BP neural network model r

表3 各水質(zhì)參數(shù)BP-ANN結(jié)構(gòu)Table 3 BP-ANN structure of each water quality parameter

2.2 BP神經(jīng)網(wǎng)絡(luò)模型的驗證

在訓(xùn)練區(qū)中隨機挑選 10個樣點，以訓(xùn)練好的BP神經(jīng)網(wǎng)絡(luò)分別對各采樣點的總氮、氨氮、總磷、COD含量進行預(yù)測，得到預(yù)測結(jié)果對BP神經(jīng)網(wǎng)絡(luò)模型做進一步驗證（表4）。

從表4可以看出，驗證點BP神經(jīng)網(wǎng)絡(luò)的預(yù)測值和實測值比較接近，特別是在平均值水平上預(yù)測值和實測值非常接近，總氮在誤差0.1的情況下可以認為預(yù)測值和實測值的平均值相等，氨氮和總磷在誤差 0.01的情況下可以認為預(yù)測值與實測值的平均值相等，COD的預(yù)測值與實測值的差距相對較大，但是 COD含量實測值的變化幅度較大，而且COD含量的國家水質(zhì)標準分類區(qū)間大，反演結(jié)果雖然與預(yù)測值有差距，仍可以在一定程度上反映密云水庫的COD含量變化情況。

為了更直觀的表現(xiàn)預(yù)測和實測值的變化情況，將實測值與預(yù)測值繪制成折線圖（圖2）。從圖2可以看出，總氮、氨氮、總磷含量的預(yù)測值與實測值比較接近，變化趨勢基本一致，預(yù)測值相對于實測值的浮動較小，預(yù)測效果較好。COD含量的預(yù)測值與實測值的變化趨勢大致相同，但是預(yù)測值相對于實測值的偏差相對較大，COD的預(yù)測效果相對較差。因此本研究所建立的4個BP神經(jīng)網(wǎng)絡(luò)模型中，總氮、氨氮、總磷具有較好的反演能力，COD的反演能力相對較差。

為進一步驗證BP神經(jīng)網(wǎng)絡(luò)模型質(zhì)量，采用納什效率系數(shù)和均方根誤差，以 10個驗證點的觀測值和預(yù)測值進行計算，結(jié)果見（表5）。

從計算結(jié)果來看，氨氮、總磷和 COD的納什系數(shù)接近1，說明這幾個水質(zhì)參數(shù)的模型質(zhì)量較好，模型可信度高，其中是氨氮和總磷納什系數(shù)較高，而總氮的納什系數(shù)更接近于 0，說明總氮的神經(jīng)網(wǎng)絡(luò)模型其模擬結(jié)果接近觀測值的平均值水平，即總氮模擬結(jié)果總體結(jié)果可信，但模擬過程誤差大。均方根誤差顯示，總磷、總氮和氨氮的均方根誤差均小于0.03，其中最小的為總磷，其次是氨氮，再次是總氮，說明總磷、總氮和氨氮的預(yù)測值和實測值較為接近。COD的均方根誤差大于1，說明 COD的真實值和預(yù)測值有偏差。綜合來看，本研究建立的氨氮和總磷的BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測效果較好，其次是氨氮，COD的神經(jīng)網(wǎng)絡(luò)模型預(yù)測效果較差。

表4 BP神經(jīng)網(wǎng)絡(luò)預(yù)測值與實測值的比較Table 4 Comparison between measured values and predicted values of BP neural network mg·L-1

圖2 BP神經(jīng)網(wǎng)絡(luò)模型結(jié)果驗證Fig.2 Verification of BP neural network Model

表5 納什效率系數(shù)、均方根誤差Table 5 NSE (Nash-Sutcliffe efficiency coefficient),RMSE(Root Mean Square Error)

2.3 密云水庫主要水質(zhì)參數(shù)的變化分析

2.3.1 主要水質(zhì)指標的年際變化

本研究基于構(gòu)建的 BP神經(jīng)網(wǎng)絡(luò)模型反演了2013—2018年密云水庫非結(jié)冰期每月的主要水質(zhì)參數(shù)，估算了每年主要水質(zhì)參數(shù)總氮、氨氮、總磷、COD的年平均值（圖3），分析2013—2018年密云水庫水質(zhì)的年際變化情況。

結(jié)果顯示，2013—2018年密云水庫的主要水質(zhì)指標總氮、氨氮、總磷、COD整體均呈下降趨勢?？偟|(zhì)量濃度在2013年最高，達到0.932 mg·L-1，2017 年最低，為0.868 mg·L-1，2015—2018 年年際波動較小，含量基本穩(wěn)定；氨氮質(zhì)量濃度在2013—2018年變化比較平緩，2015達到最高值0.018 mg·L-1，2015—2017年連續(xù)下降，2017年為最低值 0.0099 mg·L-1；總磷質(zhì)量濃度變化為先呈上升趨勢，2015年達到最高值0.027 mg·L-1，2015—2018年總磷含量表現(xiàn)為逐年小幅下降且趨勢變緩；COD變化趨勢與氨氮相似，2013—2015年呈上升的趨勢，2015年達到最高值 12.974 mg·L-1，2015—2017年快速下降，2017年為最低值6.17 mg·L-1，年際波動比較大。

從 2013—2018年間密云水庫年均水質(zhì)指標變化來看，氨氮和 COD的含量一直處于Ⅰ類水質(zhì)標準，且在Ⅰ類類別中含量較小，指標優(yōu)良；而總磷近年為Ⅱ類水質(zhì)標準，總氮為Ⅲ類水質(zhì)。總體來看，密云水庫水質(zhì)較好，且年際變化顯示近年來密云水庫各項水質(zhì)指標均呈下降趨勢，水質(zhì)整體呈上升趨勢。同時由年際變化可以看出，2018年總氮、氨氮、COD比上年有所增長。因此，密云水庫的水質(zhì)風(fēng)險仍然存在，保護與監(jiān)測工作不容忽視。

2.3.2 主要水質(zhì)參數(shù)的季節(jié)變化分析

本研究選取的是 2013—2018非結(jié)冰期遙感影像，因冬季密云水庫區(qū)域人為活動較少，生物生長緩慢，庫區(qū)水質(zhì)比較穩(wěn)定，故只研究春、夏、秋 3個季節(jié)主要水質(zhì)指標的變化狀況。為了剔除偶然年份的影響，得到研究時段密云水庫水質(zhì)參數(shù)的季節(jié)性變化特征，本研究將反演結(jié)果按季節(jié)求得多年平均值，獲取了4個水質(zhì)指標不同季節(jié)的變化分布圖（圖4）。

圖3 主要水質(zhì)參數(shù)的年際變化Fig.3 Interannual variation of major water quality parameters

由圖4可以看出，密云水庫主要水質(zhì)參數(shù)總氮、氨氮、總磷、COD的季節(jié)變化趨勢基本一致，均表現(xiàn)為春季濃度最高，秋季次之，夏季濃度最低，即密云水庫夏季水質(zhì)最優(yōu)，后因水資源的補充不足，水質(zhì)參數(shù)含量升高，導(dǎo)致秋春季水質(zhì)較差，春季達到最差。

2.3.3 密云水庫主要水質(zhì)參數(shù)的空間變化特征分析

本研究在空間柵格尺度上估算了 2013—2018年密云水庫水質(zhì)參數(shù)總氮、氨氮、總磷、COD的含量，揭示了主要水質(zhì)參數(shù)的空間分異與變化特征。

2.3.3.1 水質(zhì)參數(shù)總氮的空間變化特征

密云水庫總氮含量按照國家水質(zhì)分類標準總體以0.7 5mg·L-1和1 mg·L-1為分界點劃分為三類區(qū)域。從反演結(jié)果來看，總氮質(zhì)量濃度比較集中分布于0.75—1 mg·L-1這個區(qū)段。為了更加顯著突出總氮含量的空間變化，將這一區(qū)段再以 0.83 mg·L-1和0.91mg·L-1為分界點細分，共分為五級（圖5）。

從分級結(jié)果來看，北部區(qū)域的總氮含量高于其他區(qū)域。其中，2013年水庫的總氮質(zhì)量濃度主要集中在0.91—1 mg·L-1區(qū)域，中部和北部地區(qū)含量高于其他區(qū)域。2014年總氮含量明顯下降，集中區(qū)域轉(zhuǎn)變?yōu)?.83—0.91 mg·L-1區(qū)域，西部和北部岸邊含量顯著高于其他區(qū)域。2015年水域面積縮小，總氮含量明顯降低，東北部含量高于其他區(qū)域。到2016年水域面積擴大，雖然總氮含量仍然集中在0.83—0.91 mg·L-1區(qū)段，但是相比2014年整體降低了，北部水域總氮含量高于其他區(qū)域。2017年水域面積進一步擴大，總氮含量在0.75—0.83 mg·L-1區(qū)段分布面積擴大。2018年水域面積達到近年最大，總氮含量相比 2017年有所上升，主要分布區(qū)域位于東南部地區(qū)，這些區(qū)域遙感影像的云量偏高，這會造成估算值比實際值偏高。總體來看，總氮含量較高的水域面積呈減小趨勢，且集中分布于密云水庫北部和東南部。

2.3.3.2 水質(zhì)參數(shù)氨氮的空間變化特征

2013—2018年 密云水庫的氨氮質(zhì)量濃度小于0.03 mg·L-1。為了顯示密云水庫氨氮含量的空間差異性，以 0.006、0.012、0.018、0.024 mg·L-1為分界點將氨氮含量分為五級（圖6）。2013—2015年密云水庫的氨氮含量以0.012—0.018 mg·L-1為主，北部和東北部含量較高。2014年西部水域氨氮含量明顯降低，2015年東北部0.018—0.024 mg·L-1的氨氮水域面積不斷擴大。2016年之后密云水庫氨氮含量明顯下降，以0.012—0.018 mg·L-1為主；2017年氨氮含量進一步下降，以0.006—0.012 mg·L-1為主；2018年東部地區(qū)氨氮含量較上年有所上升。從密云水庫氨氮的空間分布上看，2015年較氨氮含量較高的區(qū)域分布面積最大，2015年后氨氮含量呈下降趨勢，2017年氨氮含量較低的區(qū)域分布面積最大。從水質(zhì)參數(shù)氨氮標準來看，2013—2018年密云水庫的水質(zhì)良好，氨氮含量處于Ⅰ類水質(zhì)標準，且遠低于Ⅰ類閾值。

圖4 主要水質(zhì)參數(shù)季節(jié)變化Fig.4 Seasonal changes in major water quality parameters

圖5 2013—2018密云水庫水質(zhì)參數(shù)總氮空間分布圖Fig.5 Spatial distribution of total nitrogen in Miyun Reservoir during the period of 2013-2018

2.3.3.3 水質(zhì)參數(shù)總磷的空間變化特征

圖6 2013—2018密云水庫水質(zhì)參數(shù)氨氮空間分布圖Fig.6 Spatial distribution of ammonia nitrogen in Miyun Reservoir during the period of 2013-2018

圖7 2013—2018密云水庫水質(zhì)參數(shù)總磷空間分布圖Fig.7 Spatial distribution of total phosphorus in Miyun Reservoir during the period of 2013-2018

總磷按國家水質(zhì)分類標準標分為四級（圖7）。2013年密云水庫的總磷質(zhì)量濃度集中在 0.01—0.025 mg·L-1區(qū)域，整體含量較低。2014年總磷含量略有升高，集中在 0.025—0.05 mg·L-1區(qū)域，北部和南部區(qū)域含量較高。2015年總磷含量較高區(qū)域主要分布于東北部。2016年密云水庫北部0.025—0.05 mg·L-1總磷分布面積擴大。而到2017年之后北部區(qū)域總磷含量明顯降低；2018年南部地區(qū)總磷含量升高?？傮w來看，密云水庫的總磷含量處于Ⅱ類和Ⅲ類水質(zhì)標準。因此，密云水庫仍需加強上中游地區(qū)污染物總磷的排放和監(jiān)管。

2.3.3.4 水質(zhì)參數(shù)COD的空間變化特征

為了體現(xiàn)COD質(zhì)量濃度的空間變化，以7.5、15、20、30 mg·L-1為分界點將COD含量劃分為五級（圖8）。2013年密云水庫的COD質(zhì)量濃度以7.5—15 mg·L-1為主，東北部區(qū)域COD含量較高。東北部COD含量較高的區(qū)域隨著時間推移而不斷擴大，2015年COD含量較高的區(qū)域擴到最大。2016年東北部COD含量最高等級的區(qū)域迅速縮小，而20—30 mg·L-1的COD面積不斷擴大。到2017年COD含量全面下降，COD含量主要集中在0—7.5 mg·L-1區(qū)段。2018年東部COD含量上升為 7.5—15 mg·L-1。總體而言，密云水庫的COD含量處于Ⅰ類和Ⅱ類水質(zhì)標準，符合飲用水源要求。

圖8 2013—2018密云水庫水質(zhì)參數(shù)COD空間分布特征Fig.8 Spatial distribution of COD in Miyun Reservoir during the period of 2013-2018

3 討論

分析結(jié)果表明，密云水庫的主要水質(zhì)參數(shù)總氮、氨氮、總磷、COD含量雖然年際之間有波動，但整體呈下降趨勢，表明密云水庫的水質(zhì)在不斷改善，水質(zhì)保護工作取得了顯著成果。從各項水質(zhì)參數(shù)的季節(jié)平均值來看，夏季含量最低，秋季次之，春季最高，可能是由于夏季降雨較多，密云水庫的水源得到補給，水體得到更新，導(dǎo)致水質(zhì)相比其他季節(jié)更佳。密云水庫主要水質(zhì)指標的空間分布顯示，庫區(qū)北部、和東部水質(zhì)指標偏高，這種分布可能與北部和東部村莊密集，以及密云水庫兩大入庫河流：潮河和白河，從北部入庫有關(guān)。秦麗歡等（2017）研究表明，深水區(qū)和朝河及白河的過度區(qū)域總磷含量較高，這與本研究所得空間分布吻合。此外，班靜雅等（2019）研究表明南水北調(diào)工程使得庫區(qū)水位大幅上升，庫區(qū)周邊大面積土壤被永久淹沒，產(chǎn)生了土壤污物溶出，模擬結(jié)果表明對密云水庫水質(zhì)影響范圍小、程度低。筆者認為，這也可能是導(dǎo)致北部和東部水域指標含量較高的另一個原因?？傮w來看，密云水庫水質(zhì)較好，氨氮和COD兩個參數(shù)為Ⅰ類水質(zhì)，總磷為Ⅱ類，總氮為Ⅲ類，各項水質(zhì)參數(shù)指標滿足飲用水源要求。然而密云水庫水質(zhì)監(jiān)測和保護工作不能放松，總磷和總氮在個別年份也存在小幅超標的現(xiàn)象。

本研究開展的基于BP神經(jīng)網(wǎng)絡(luò)算法反演密云水庫水質(zhì)參數(shù)的研究達到了一定精度，但也存在一些不足。由于Landsat 8數(shù)據(jù)有限，再加上密云地區(qū)特殊的氣象條件，一些時段的遙感影像受到霧霾和云層的影響無法使用，導(dǎo)致遙感影像無法被充分利用，以后的研究需要進一步探索針對水質(zhì)反演的大氣校正模型和遙感影像的去云技術(shù)，從而減小云量對于輻射值的影響，使遙感影像得到充分利用。另外，本研究采用的多光譜遙感影像分辨率為 30 m×30 m，而記錄采樣點坐標的手持GPS的精度為15 m，因此兩者的配準會產(chǎn)生一定的誤差，期待以后的研究通過技術(shù)手段的進步，提高遙感影像與采樣點的配準程度，使數(shù)據(jù)能夠更加準確的反映水庫的水質(zhì)指標的變化。

4 結(jié)論

本研究基于Landsa 8遙感數(shù)據(jù)與地面監(jiān)測水質(zhì)數(shù)據(jù)，采用遙感和GIS技術(shù)，通過構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型，實現(xiàn)了密云水庫2013—2018總磷、總氮、氨氮和COD（化學(xué)需氧量）4個水質(zhì)參數(shù)的反演，并對水質(zhì)參數(shù)的時空變化做了分析，得到如下結(jié)論：

（1）主要水質(zhì)參數(shù)的Landsat 8敏感波段為：總氮為1、4波段，氨氮為1—7波段，總磷為1、3—7波段，COD為2—5波段。對BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測結(jié)果驗證顯示，密云水庫水質(zhì)指標的預(yù)測值與實測值在平均值水平上非常接近，納什系數(shù)和均方根誤差評價結(jié)果表明，研究構(gòu)建的總磷和氨氮的BP神經(jīng)網(wǎng)絡(luò)模型質(zhì)量較好，可以很好實現(xiàn)水質(zhì)反演，而氨氮和COD模型質(zhì)量不佳，但在平均值水平上仍然具有準確性。

（2）密云水庫主要水質(zhì)參數(shù)總氮、氨氮、總磷、COD的含量在2013—2018年總體呈下降趨勢，年均值顯示，氨氮和 COD兩個參數(shù)處于Ⅰ類水質(zhì)，總磷為Ⅱ類，總氮為Ⅲ類，主要各項水質(zhì)指標滿足飲用水源要求，需注意加強總磷和總氮的控制。

（3）從季節(jié)變化特征分析顯示，4個水質(zhì)均表現(xiàn)為春季最高、秋季次之、夏季最低。

（4）密云水庫主要水質(zhì)參數(shù)的空間差異性較大，總體來看，水庫北部和東部水域的4個水質(zhì)參數(shù)含量相對較高，這與北部和東部村莊分布，以及密云水庫的兩大入庫河流：潮河和白河，從北部入庫有關(guān)。