趙力，盧修元，譚海，馬天浩

(1.四川農業(yè)大學 水利水電學院，四川 雅安 625014；2.自然資源部國土衛(wèi)星遙感應用中心，北京 100048；3.遼寧工程技術大學 測繪與地理科學學院，遼寧 阜新 123000)

0 引言

長期以來農業(yè)面源污染形勢嚴峻，而化肥施用量高居不下是主因。2014—2018年，全國單位耕地化肥施用量由0.83 t·ha-1增長為0.93 t·ha-1[1]，化肥中的氮、磷營養(yǎng)鹽進入水體，導致水體富營養(yǎng)化，影響水生環(huán)境，危害人體健康及動植物的生存，阻礙當?shù)厣鐣洕沙掷m(xù)發(fā)展。因此，快速監(jiān)測總氮(TN)、總磷(TP)的含量，客觀地反映水體富營養(yǎng)化程度，為污染水體治理提供依據(jù)是十分必要的[2]。傳統(tǒng)的水體氮、磷濃度監(jiān)測以點位監(jiān)測為主，由人工到現(xiàn)場采集水樣并進行化學分析，該方法耗時費力，且樣點稀疏，只能反映取樣點位的水質情況，對未取樣位置，無法準確判定其水質。而衛(wèi)星遙感技術通過經驗、半經驗模型或是分析模型建立衛(wèi)星遙感觀測信息與地面點位實測數(shù)據(jù)間的關系[3]，可以高效監(jiān)測水質。潘邦龍等[4]利用環(huán)境一號衛(wèi)星(HJ-1)超光譜遙感數(shù)據(jù)，建立了TN的三波段回歸克里格模型；王小平等[5]采用Landsat-8 OLI數(shù)據(jù)，建立了艾比湖流域TP的4波段反演模型；錢彬杰[6]利用Landsat-8 OLI衛(wèi)星遙感影像，建立了氨氮和TP波段反射率比值的回歸模型；王磊等[7]基于Landsat-8 OLI影像，對丹江口水庫的TP濃度進行了反演；林劍遠等[8]采用芬蘭機載高光譜成像系統(tǒng)獲取的遙感數(shù)據(jù)，運用波段比值對嘉興市水體的TP和TN進行反演，為內陸城市河網水質監(jiān)測提供了參考。上述學者對水體氮、磷鹽的遙感監(jiān)測提供了有益探索。然而，由于不同建模方式的原理及反映復雜映射關系的能力不同，導致模型預測的效果差異較大，實測值與預測值的決定系數(shù)R2低至0.50，高至0.99[9]。此外，目前國內外遙感水質監(jiān)測研究主要利用中低分辨率遙感數(shù)據(jù)，多數(shù)研究僅針對較大的內陸湖泊或海洋等水體[10]。

利用遙感進行水質監(jiān)測中的建模方法，以數(shù)學統(tǒng)計法最為常見[11]，但因其難以反映復雜的水體光學性質，導致所得模型穩(wěn)定性較差[12]。近年來，以BP神經網絡為代表的機器學習因為具有自適應、自學習、高效率和容錯性等優(yōu)點，在遙感領域也得到了廣泛使用[13-15]。而XGBoost(extreme gradient boosting)模型是提升算法(boosting)的最新研究成果之一，不僅具有較強的自我學習能力和映射能力，同時還彌補了BP神經網絡無法反映輸入因子的重要性的問題，已被廣泛用于圖像處理[16]和大數(shù)據(jù)處理[17]等領域。就數(shù)據(jù)源而言，目前常用的MODIS、MERIS、Landsat、HJ-1等中低分辨率數(shù)據(jù)限制了研究尺度。而2013年我國發(fā)射升空的高分一號衛(wèi)星(GF-1)搭載有高空間分辨率的寬幅多光譜相機(WFV)，為小尺度水體的遙感監(jiān)測提供了基礎條件[18-19]。盡管已有基于GF-1數(shù)據(jù)并采用機器學習方法反演水質的研究，但將GF-1數(shù)據(jù)與XGBoost模型聯(lián)合使用進行水質監(jiān)測的研究還鮮有報道。

本研究以長江上游農業(yè)生產水匯集的白水湖水庫為研究區(qū)，以水體TN、TP為研究對象，分別采用XGBoost、BP神經網絡2種機器學習方法，構建水質遙感反演模型，并與數(shù)學統(tǒng)計模型比較，旨在選出適用于研究區(qū)域的水質監(jiān)測模型，為長期、便捷、大范圍的水質監(jiān)測提供技術支撐。

1 材料與方法

1.1 研究區(qū)概況

白水湖水庫位于長江上游，川西平原北部(31°30′55″N～31°32′10″N，104°14′45″～104°16′07″E)，屬中亞熱帶濕潤季風氣候，冬干春旱，夏季旱澇交錯，秋多陰雨，水域面積2 km2，庫容1 672萬m3，平均水深10 m，湖區(qū)四面環(huán)田，在雨水沖刷和地表徑流作用下，農田里過量營養(yǎng)鹽直接進入水體，造成水環(huán)境污染。

1.2 技術路線

采用ENVI 5.3軟件校正GF-1 WFV1圖像，提取研究水域4個波段的反射率?；趯崪yTN和TP數(shù)據(jù)，通過計算水質參數(shù)與波段及波段組合的Pearson相關系數(shù)，識別各個水質參數(shù)的敏感波段或波段組合，將其作為模型的輸入因子。隨機選取70%的樣本作為訓練集，剩余30%樣本作為測試集，利用Python 3.7軟件構建XGBoost模型，利用MATLAB 軟件構建數(shù)學統(tǒng)計模型和BP神經網絡模型，以R2和均方根誤差(RMSE)為依據(jù)，比選模型。反演TN、TP，運用ArcGIS 10.2軟件制作水質專題圖，并分析TN和TP空間分布特征。具體技術路線如圖1所示。

圖1 技術路線

1.3 數(shù)據(jù)采集及預處理

1)水質參數(shù)測量。2019年8月20日，在白水湖采集表層水樣87個，樣點分布如圖2所示。TN采用過硫酸鉀氧化紫外分光光度法測定(HJ 636-2012)，TP采用鉬酸銨分光光度法(GB 11893-89)測定。

圖2 研究水域樣點分布

2)衛(wèi)星遙感影像數(shù)據(jù)獲取與處理。GF-1 WFV影像空間分辨率為16 m，重訪周期為4 d，具有B1(0.45～0.52 nm，藍)、B2(0.52～0.59 nm，綠)、B3(0.63～0.69 nm，紅)、B4(0.77～0.89 nm，近紅外)4個波段。本文采用8月11日GF-1 WFV準同步數(shù)據(jù)。利用ENVI軟件對GF1-WFV1數(shù)據(jù)進行預處理。DEM數(shù)據(jù)為GDEMV2 30 m。采用3次卷積插值完成正射校正，以15 m空間分辨率的Landsat-8 OLI影像為基準圖像完成影像的幾何精校正，輻射定標采用2018年GF-1 WFV1絕對輻射定標系數(shù)，F(xiàn)LAASH模塊進行大氣校正，部分輸入參數(shù)如表1所示。使用歸一化水體指數(shù)(normalized difference water index，NDWI)(式(2))對研究區(qū)水體進行掩膜處理，得到水體影像。

(1)

式中：Green為GF-1數(shù)據(jù)綠光波段反射率；NIR為近紅外波段反射率。

表1 FLASSH 大氣校正模塊輸入值

2 數(shù)據(jù)分析

2.1 模型輸入因子分析

衛(wèi)星傳感器上記錄的與水體信息有關的總輻射亮度主要包括水面散射輻射、水中散射輻射以及底部反射輻射3個部分。當水體污染物、水面粗糙度或微波發(fā)生變化時，水體散射和反射狀態(tài)等也會變化，總輻射亮度隨之改變，導致總輻射亮度與水質參數(shù)改變。正由于水體成分和水體狀態(tài)會影響水體的光譜特征[20]，可以通過波段組合消除或減少與水質無關的信息，保留或放大顯著相關的信息，提高反演模型精度。如通過剔除對有機污染物濃度不敏感的近紅外波段，利用敏感度高的可見光波段建立模型，可以提高水體有機污染物反演精度[21]。此外，已有研究也表明，波段組合不僅可以消除部分大氣影響[22]，還可以減少水面粗糙度變化的影響[23]，同時也能消除微波變化導致的吸收系數(shù)和后向散射系數(shù)干擾，減少水體反射率二向反射等問題[24]。所以，可以通過波段組合，減少污染物、粗糙度、微波變化等的干擾，從而把目標水質參數(shù)的特性表現(xiàn)出來。

將校正后的波段反射率信息經過單波段、雙波段和三波段間的減法(如B1-B2、B3-B4等)、比值法(如B1/B2、(B1-B3)/(B2-B4)等)、對數(shù)法(ln(B1)、ln(B1/B2) 等)，組合成254個波段反射率信息，編號為Cn(n=1～254)，利用SPSS 21軟件分別計算各波段組合信息與TN、TP的Person相關系數(shù)，結果如圖3所示。取極顯著相關(p<0.01)且相關系數(shù)最大的波段組合作為輸入因子，則TN的輸入因子為波段組合C34、C33、C45、C31，TP為C156、C47、C174、C144，詳見表2。由表2可知，TN對B2、B3和B4的波段組合最敏感，TP對B1和B2的波段組合最敏感。Koponen等[25]的研究也表明，采用多波段組合可以在一定程度上消除水體污染物、表面粗糙度以及微波的時空變化等干擾，與單波段相比，更能真實反映水質參數(shù)對光譜曲線的影響。本文的結論與Koponen的研究相類似。

圖3 波段組合信息與TN、TP 的Person相關系數(shù)

表2 模型輸入因子

2.2 TN反演模型的建立與評價

1)XGBoost。XGBoost是一種基于決策樹的監(jiān)督模型，是集成學習方法的一種，由華盛頓大學的Chen等[26]提出，其核心在于新的模型在相應損失函數(shù)梯度方向建立，修正殘差的同時控制模型的復雜度。本研究以相關性最高的波段組合C34、C33、C45和C31作為輸入因子，以實測TN作為輸出，XGBoost模型基于樹模型進行提升計算，采用線性回歸目標函數(shù)，多次調參對比結果，最終設定XGBoost最大深度為3，限定孩子節(jié)點中最小的樣本權重和為3，詳細參數(shù)設定見表3。訓練所得模型中波段組合的重要性C34>C45>C33>C31分別為32、31、23和0，包含波段B2、B3和B4，模型的決定系數(shù)R2為0.970，RMSE為1.24×10-2mg·L-1。

表3 XGBoost的TN反演模型參數(shù)

2)BP神經網絡。BP神經網絡由Webos在1974年提出，主要特點是信號向前傳遞，誤差反向傳遞，調整層間的權值、閾值以達到預設精度要求[27]。BP神經網絡構建中，隱含層神經元節(jié)點數(shù)的選擇直接影響網絡對復雜問題的映射能力[28]，隱含層結點數(shù)過少，則網絡收斂速度減慢，精度降低，結點數(shù)過多，則導致計算量增加，模型過擬合，泛化能力降低。但目前尚無確定隱含層最佳結點數(shù)的直接方法[29]。因此，通常采用的方法是多次實驗對比分析。本實驗經過60個訓練樣本訓練不同節(jié)點數(shù)的單隱含層BP神經網絡模型，節(jié)點數(shù)3～10時，模型的R2和RMSE如表4所示，結果表現(xiàn)為9節(jié)點時R2最高，4節(jié)點次之。分別計算剩余27個測試樣本在8種不同節(jié)點數(shù)的R2和RMSE，結果如表5所示，R2表現(xiàn)為4節(jié)點最高，7節(jié)點次之。綜上，將單層4節(jié)點BP神經網絡模型作為TN的BP神經網絡最優(yōu)模型，此時模型的R2為0.769，RMSE為3.47×10-2mg·L-1。

3)數(shù)學統(tǒng)計模型。取訓練樣本(n=60)，以波段組合反射率C34、C33、C45、C31為自變量，TN濃度為因變量，構建多元線性回歸模型，模型中波段組合的權重為C31>C45>C34>C33，分別為0.518 12、0.144 13、0.138 03和0，模型的R2和RMSE分別為0.755和3.50×10-2mg·L-1，詳見式(2)。

表4 BP神經網絡的TN模型精度(n=60)

表5 BP神經網絡的TN預測精度(n=27)

TNST=-0.138 03×C34-0.144 13×C45-0.518 12×C31+0.765 847

(2)

代入表2中因子，則反演TN的數(shù)學統(tǒng)計模型包含4個波段，如式(3)所示。

(3)

4)TN模型評價。以30%(27個)的采樣點數(shù)據(jù)作為模型測試樣本，將波段組合反射率分別輸入XGBoost、單層4節(jié)點BP神經網絡和數(shù)學統(tǒng)計模型，輸出TN預測值，作預測值與實測值的散點圖(圖4)。由圖4可知，XGBoost模型(XG)、BP神經網絡模型(BP)、數(shù)學統(tǒng)計模型(ST)的R2分別為0.835 0、0.801 4、0.756 7；各模型的RMSE值分別為3.25×10-2、2.97×10-2、3.47×10-2mg·L-1。XGBoost模型的R2比BP神經網絡和數(shù)學統(tǒng)計模型分別高了4.19%和10.35%，RMSE比數(shù)學統(tǒng)計模型降低了6.34%。由此可見，在3種模型中，XGBoost的TN反演模型具有最高的R2和較低的RMSE，對研究區(qū)TN濃度的監(jiān)測效果明顯優(yōu)于BP神經網絡和數(shù)學統(tǒng)計模型。

圖4 TN預測值與實測值比較

2.3 TP反演模型的建立與評價

1)XGBoost。以相關性最高的波段組合C156、C47、C144和C174作為輸入因子，以實測TP作為輸出，XGBoost基于樹模型進行提升計算，采用線性回歸目標函數(shù)，對比多次調參結果，選定XGBoost最大深度為3，限定孩子節(jié)點中最小的樣本權重和為1，詳細參數(shù)設定如表6所示。訓練所得模型中波段組合的重要性C47>C156>C174>C144，分別為96、65、45和8，包含波段B1和B2，模型的R2為0.896 6，RMSE為1.97×10-3mg·L-1。

表6 XGBoost的TP反演模型參數(shù)

2)BP神經網絡。以隨機抽取的70%的數(shù)據(jù)作為訓練樣本，計算不同節(jié)點數(shù)的單隱含層BP神經網絡模型的R2和RMSE，結果如表7所示，R2表現(xiàn)為10節(jié)點>9節(jié)點>7節(jié)點>其他。以剩余30%的數(shù)據(jù)作為測試樣本，分別計算8種不同節(jié)點數(shù)模型的R2和RMSE，結果如表8所示，R2表現(xiàn)為5節(jié)點>6節(jié)點>10節(jié)點>其他。因此，將單層10節(jié)點BP神經網絡模型作為TP的BP神經網絡最優(yōu)模型，R2為0.911，RMSE為1.96×10-3mg·L-1。

表7 BP神經網絡的TP建模精度(n=60)

表8 BP神經網絡的TP預測精度(n=27)

3)統(tǒng)計學模型。以波段組合反射率C156、C47、C144和C174為自變量，TP濃度為因變量，抽取70%的數(shù)據(jù)作為模型訓練樣本，構建多元線性回歸模型，模型中波段組合反射率的權重為C156>C47>C174>C144，分別為0.354 57、0.146 55、0.026 76、0.022 17，與XGBoost模型的重要性存在差異，模型的R2和RMSE分別為0.696和3.49×10-3mg·L-1，模型如式(4)所示。

TPST=-0.355 7×C156-0.146 55×
C47-0.022 17×C144+0.026 76×
C174-0.084 13

(4)

代入表2中因子，則反演TP的數(shù)學統(tǒng)計模型包含B1和B2波段，模型為式(5)。

(5)

4)TP模型評價。將訓練樣本(n=27個)的波段組合反射率分別輸入XGBoost、單層10節(jié)點BP神經網絡和線性數(shù)學統(tǒng)計模型中，輸出TP預測值，作預測值與實測值的散點圖(圖5)。由圖5可知，XGBoost模型(XG)、BP神經網絡模型(BP)、數(shù)學統(tǒng)計模型(ST)的R2分別為：0.896 6、0.754 2、0.746 2；各模型的RMSE值分別為1.97×10-3、0.974×10-3、2.67×10-3mg·L-1。XGBoost模型與BP神經網絡和數(shù)學統(tǒng)計的R2相比，分別高了18.88%和20.16%。與數(shù)學統(tǒng)計模型相比,XGBoost模型的RMSE低了26.22%。由此可見，XGBoost模型能有效提高預測精度，對研究區(qū)TP濃度具有更好的反演效果，更適用于研究水體的TP監(jiān)測。

2.4 氮、磷遙感評價模型應用

以水質反演效果最好的XGBoost模型為基礎，將對應反射率數(shù)據(jù)輸入模型，反演得到TN、TP值。利用ArcGIS分別制作TN、TP的水質專題圖，結果見圖6。由圖6(a)可知,TN濃度進水區(qū)域高于出水區(qū)域，近岸高于湖心，符合反硝化作用下的庫區(qū)水體TN分布特征。由圖6(b)可知，TP濃度分布與TN相反，呈現(xiàn)出水區(qū)域高于進水區(qū)域，湖心高于近岸的趨勢，原因可能是水體中磷酸鹽被底泥吸附，吸附量隨含氧量降低而降低，出水口和湖心底泥含氧量較低，磷酸鹽釋放，導致水體TP增加。

圖5 TP預測值與實測值比較圖

圖6 TN、TP空間分布

3 結束語

本研究以農業(yè)生產水匯集的白水湖水庫為研究區(qū)，利用校正后的GF-1 WFV影像，采用XGBoost模型和BP神經網絡模型2種機器學習方法，反演水體TN和TP，并與數(shù)學統(tǒng)計模型比較，得到以下結論。

1)3種模型中，XGBoost模型更適合預測TN、TP，其預測值與實際值的R2分別可達到0.835、0.897。理論上而言，提高布點均勻度還可以進一步提高XGBoost模型反演TN、TP的精度。

2)同一水質參數(shù)，采用不同建模方法建模時，其輸入因子的重要性不同。TN差異較大時，XGBoost模型的C31最不重要，而數(shù)學統(tǒng)計模型中C31最重要；TP差異較小時，XGBoost模型中C47最重要，C156次之，而數(shù)學統(tǒng)計模型中C156最重要，C47次之，其他不變。對不同水質指標，敏感波段不同，TN的敏感波段為綠波、紅波和近紅外波段，TP的敏感波段為藍波和綠波。前人研究表明，水深影響綠波段反射率，懸浮物濃度增加會產生“紅移”現(xiàn)象，因此，可考慮結合水深和懸浮物濃度建立TN、TP反演模型，以進一步提高精度。

本研究結果表明，采用GF-1 WFV信息，結合XGBoost模型能準確反演水體TN、TP情況。本研究為農業(yè)生產水匯集區(qū)域的水體監(jiān)測提供了方法，也為其他水體利用遙感技術進行水質的總氮、總磷監(jiān)測提供了思路。