利用高分一號衛(wèi)星與XGBoost模型的水體總氮和總磷監(jiān)測技術(shù)

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

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

0 引言

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 技術(shù)路線

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

圖1 技術(shù)路線

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

1)水質(zhì)參數(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，藍(lán))、B2(0.52～0.59 nm，綠)、B3(0.63～0.69 nm，紅)、B4(0.77～0.89 nm，近紅外)4個波段。本文采用8月11日GF-1 WFV準(zhǔn)同步數(shù)據(jù)。利用ENVI軟件對GF1-WFV1數(shù)據(jù)進(jìn)行預(yù)處理。DEM數(shù)據(jù)為GDEMV2 30 m。采用3次卷積插值完成正射校正，以15 m空間分辨率的Landsat-8 OLI影像為基準(zhǔn)圖像完成影像的幾何精校正，輻射定標(biāo)采用2018年GF-1 WFV1絕對輻射定標(biāo)系數(shù)，F(xiàn)LAASH模塊進(jìn)行大氣校正，部分輸入?yún)?shù)如表1所示。使用歸一化水體指數(shù)(normalized difference water index，NDWI)(式(2))對研究區(qū)水體進(jìn)行掩膜處理，得到水體影像。

(1)

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

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

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

2.1 模型輸入因子分析

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

將校正后的波段反射率信息經(jīng)過單波段、雙波段和三波段間的減法(如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相關(guān)系數(shù)，結(jié)果如圖3所示。取極顯著相關(guān)(p<0.01)且相關(guān)系數(shù)最大的波段組合作為輸入因子，則TN的輸入因子為波段組合C34、C33、C45、C31，TP為C156、C47、C174、C144，詳見表2。由表2可知，TN對B2、B3和B4的波段組合最敏感，TP對B1和B2的波段組合最敏感。Koponen等[25]的研究也表明，采用多波段組合可以在一定程度上消除水體污染物、表面粗糙度以及微波的時空變化等干擾，與單波段相比，更能真實反映水質(zhì)參數(shù)對光譜曲線的影響。本文的結(jié)論與Koponen的研究相類似。

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

表2 模型輸入因子

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

1)XGBoost。XGBoost是一種基于決策樹的監(jiān)督模型，是集成學(xué)習(xí)方法的一種，由華盛頓大學(xué)的Chen等[26]提出，其核心在于新的模型在相應(yīng)損失函數(shù)梯度方向建立，修正殘差的同時控制模型的復(fù)雜度。本研究以相關(guān)性最高的波段組合C34、C33、C45和C31作為輸入因子，以實測TN作為輸出，XGBoost模型基于樹模型進(jìn)行提升計算，采用線性回歸目標(biāo)函數(shù)，多次調(diào)參對比結(jié)果，最終設(shè)定XGBoost最大深度為3，限定孩子節(jié)點中最小的樣本權(quán)重和為3，詳細(xì)參數(shù)設(shè)定見表3。訓(xùn)練所得模型中波段組合的重要性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神經(jīng)網(wǎng)絡(luò)。BP神經(jīng)網(wǎng)絡(luò)由Webos在1974年提出，主要特點是信號向前傳遞，誤差反向傳遞，調(diào)整層間的權(quán)值、閾值以達(dá)到預(yù)設(shè)精度要求[27]。BP神經(jīng)網(wǎng)絡(luò)構(gòu)建中，隱含層神經(jīng)元節(jié)點數(shù)的選擇直接影響網(wǎng)絡(luò)對復(fù)雜問題的映射能力[28]，隱含層結(jié)點數(shù)過少，則網(wǎng)絡(luò)收斂速度減慢，精度降低，結(jié)點數(shù)過多，則導(dǎo)致計算量增加，模型過擬合，泛化能力降低。但目前尚無確定隱含層最佳結(jié)點數(shù)的直接方法[29]。因此，通常采用的方法是多次實驗對比分析。本實驗經(jīng)過60個訓(xùn)練樣本訓(xùn)練不同節(jié)點數(shù)的單隱含層BP神經(jīng)網(wǎng)絡(luò)模型，節(jié)點數(shù)3～10時，模型的R2和RMSE如表4所示，結(jié)果表現(xiàn)為9節(jié)點時R2最高，4節(jié)點次之。分別計算剩余27個測試樣本在8種不同節(jié)點數(shù)的R2和RMSE，結(jié)果如表5所示，R2表現(xiàn)為4節(jié)點最高，7節(jié)點次之。綜上，將單層4節(jié)點BP神經(jīng)網(wǎng)絡(luò)模型作為TN的BP神經(jīng)網(wǎng)絡(luò)最優(yōu)模型，此時模型的R2為0.769，RMSE為3.47×10-2mg·L-1。

3)數(shù)學(xué)統(tǒng)計模型。取訓(xùn)練樣本(n=60)，以波段組合反射率C34、C33、C45、C31為自變量，TN濃度為因變量，構(gòu)建多元線性回歸模型，模型中波段組合的權(quán)重為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神經(jīng)網(wǎng)絡(luò)的TN模型精度(n=60)

表5 BP神經(jīng)網(wǎng)絡(luò)的TN預(yù)測精度(n=27)

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

(2)

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

(3)

4)TN模型評價。以30%(27個)的采樣點數(shù)據(jù)作為模型測試樣本，將波段組合反射率分別輸入XGBoost、單層4節(jié)點BP神經(jīng)網(wǎng)絡(luò)和數(shù)學(xué)統(tǒng)計模型，輸出TN預(yù)測值，作預(yù)測值與實測值的散點圖(圖4)。由圖4可知，XGBoost模型(XG)、BP神經(jīng)網(wǎng)絡(luò)模型(BP)、數(shù)學(xué)統(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神經(jīng)網(wǎng)絡(luò)和數(shù)學(xué)統(tǒng)計模型分別高了4.19%和10.35%，RMSE比數(shù)學(xué)統(tǒng)計模型降低了6.34%。由此可見，在3種模型中，XGBoost的TN反演模型具有最高的R2和較低的RMSE，對研究區(qū)TN濃度的監(jiān)測效果明顯優(yōu)于BP神經(jīng)網(wǎng)絡(luò)和數(shù)學(xué)統(tǒng)計模型。

圖4 TN預(yù)測值與實測值比較

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

1)XGBoost。以相關(guān)性最高的波段組合C156、C47、C144和C174作為輸入因子，以實測TP作為輸出，XGBoost基于樹模型進(jìn)行提升計算，采用線性回歸目標(biāo)函數(shù)，對比多次調(diào)參結(jié)果，選定XGBoost最大深度為3，限定孩子節(jié)點中最小的樣本權(quán)重和為1，詳細(xì)參數(shù)設(shè)定如表6所示。訓(xùn)練所得模型中波段組合的重要性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神經(jīng)網(wǎng)絡(luò)。以隨機抽取的70%的數(shù)據(jù)作為訓(xùn)練樣本，計算不同節(jié)點數(shù)的單隱含層BP神經(jīng)網(wǎng)絡(luò)模型的R2和RMSE，結(jié)果如表7所示，R2表現(xiàn)為10節(jié)點>9節(jié)點>7節(jié)點>其他。以剩余30%的數(shù)據(jù)作為測試樣本，分別計算8種不同節(jié)點數(shù)模型的R2和RMSE，結(jié)果如表8所示，R2表現(xiàn)為5節(jié)點>6節(jié)點>10節(jié)點>其他。因此，將單層10節(jié)點BP神經(jīng)網(wǎng)絡(luò)模型作為TP的BP神經(jīng)網(wǎng)絡(luò)最優(yōu)模型，R2為0.911，RMSE為1.96×10-3mg·L-1。

表7 BP神經(jīng)網(wǎng)絡(luò)的TP建模精度(n=60)

表8 BP神經(jīng)網(wǎng)絡(luò)的TP預(yù)測精度(n=27)

3)統(tǒng)計學(xué)模型。以波段組合反射率C156、C47、C144和C174為自變量，TP濃度為因變量，抽取70%的數(shù)據(jù)作為模型訓(xùn)練樣本，構(gòu)建多元線性回歸模型，模型中波段組合反射率的權(quán)重為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ù)學(xué)統(tǒng)計模型包含B1和B2波段，模型為式(5)。

(5)

4)TP模型評價。將訓(xùn)練樣本(n=27個)的波段組合反射率分別輸入XGBoost、單層10節(jié)點BP神經(jīng)網(wǎng)絡(luò)和線性數(shù)學(xué)統(tǒng)計模型中，輸出TP預(yù)測值，作預(yù)測值與實測值的散點圖(圖5)。由圖5可知，XGBoost模型(XG)、BP神經(jīng)網(wǎng)絡(luò)模型(BP)、數(shù)學(xué)統(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神經(jīng)網(wǎng)絡(luò)和數(shù)學(xué)統(tǒng)計的R2相比，分別高了18.88%和20.16%。與數(shù)學(xué)統(tǒng)計模型相比,XGBoost模型的RMSE低了26.22%。由此可見，XGBoost模型能有效提高預(yù)測精度，對研究區(qū)TP濃度具有更好的反演效果，更適用于研究水體的TP監(jiān)測。

2.4 氮、磷遙感評價模型應(yīng)用

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

圖5 TP預(yù)測值與實測值比較圖

圖6 TN、TP空間分布

3 結(jié)束語

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

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

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

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