結合機器學習與高光譜遙感的城市內河水質反演

曾江超, 黃風華, 李秉政, 高榮剛

(1. 福州大學數字中國研究院(福建), 福建 福州 350108; 2. 福建省空間信息感知與智能處理重點實驗室(陽光學院), 福建 福州 350015; 3. 空間數據挖掘與應用福建省高校工程研究中心(陽光學院), 福建 福州 350015)

0 引言

水環(huán)境治理的前提是準確獲取水體的健康狀況, 進而及時調整治理策略, 保障水體質量健康穩(wěn)定[1]. 總磷(TP)、 總氮(TN)是評價水體富營養(yǎng)化的重要指標, 主要來源是化肥和洗滌廢水排放, 研究這些指標對于監(jiān)測和治理河流水體富營養(yǎng)化有著重要意義[2].

由于內陸水體成分復雜多樣, 多光譜遙感探測水體的光譜信息較弱[3]. 而高光譜遙感帶來了豐富的光譜信息, 能夠提高內陸水體水質參數的反演精度[4]. 再者城市內河寬度較窄、 水網密集, 因此空間分辨率更高的無人機影像成為內河水環(huán)境監(jiān)測的首選方案[5].

目前, 利用衛(wèi)星和無人機高光譜影像的水質反演已有大量研究[6]. Wang等[7]基于地貌學的非點源污染模型分析了長江干流TN濃度及其對水質的影響, 但是此類半經驗模型的泛用性不高. Cao等[8]采用極端梯度增強樹對中國東部湖泊進行水質反演, 證明了機器學習算法在水體監(jiān)測的潛力. Chen等[9]采用隨機森林回歸(RF)對墨西哥灣進行TP、 TN的反演, 結果表明RF算法能卓越處理此類非線性問題. 盛輝等[10]采用改進的支持向量回歸(SVR)對濰河水域取得了較好的反演效果, 并在山東半島峽山水庫測試了該模型的泛化能力. 以上研究均表明, 采用機器學習算法進行水質反演具有較高的可行性.

本研究基于高光譜相機采集的影像, 利用機器學習算法, 結合福建省福州市晉安河的實測光譜數據, 進行TP和TN模型構建和水質反演, 從而既能豐富城市內河水質遙感反演的理論研究, 又可對內河的水體富營養(yǎng)化監(jiān)測提供重要參考價值.

1 材料與方法

1.1 研究區(qū)概況

福建省福州市晉安河歷史悠久, 貫穿晉安區(qū)、 鼓樓區(qū)和臺江區(qū), 其河流水質的好壞直接影響福州市的整體市容市貌和居民生活. 因此本研究選擇晉安河為研究區(qū), 由于研究區(qū)域長度將近7 km, 整條河流水質反演圖效果不直觀, 所以僅展示圖1晉安河東門村河段(26.09°～26.10°N, 119.310°～119.312°E)約800 m的水質參數反演結果.

圖1 研究區(qū)與采樣點示意圖Fig.1 Schematic diagram of study area and sampling points資料來源: 水經微圖(圖幅號G50H092170, G50H092171, G50H093171)http://www.rivermap.cn/

1.2 水質參數樣本采集

表1 晉安河取樣點實測水質數據統(tǒng)計

1.3 水面實測光譜數據采集

由于此次無人機航飛晉安河采用的高光譜相機波段范圍為400～1 000 nm, 波長總數為176個, 故使用Field Spec4手持光譜儀在波段為400～1 000 nm內按照3.4 nm間隔采樣. 在標準取水器采水之前, 使用手持光譜儀采用45度傾斜法采集采樣點處水面光譜信息, 單個采樣點重復采集5次, 均值作為水體反射值, 每次測定前需要采集自然光和灰板數據用于數據校正. 圖2為實測數據繪制的光譜反射率曲線, 其中波長用λ表示, 反射率物理量用R表示.

圖2 樣本點處實測光譜信息 Fig.2 Measured spectral information at sample points

1.4 無人機高光譜影像采集與處理

根據研究區(qū)域的地形制定對應航線, 通過DJI M600 Pro六旋翼無人機搭載雙利合譜GaiaSky-mini2高光譜相機獲取晉安河高光譜影像, 光譜范圍為400～1 000 nm, 光譜分辨率和空間分辨率為8 nm和0.32 m. 由于較強的水汽吸收, 波長大于900 nm時信噪比低, 反射率波動較大, 故反演建模和可視化時選擇400～900 nm波段. 而且無人機高光譜影像反射率振幅大, 本研究采用Symlet小波變換對無人機高光譜影像各像素點反射率進行處理[13]. 其中選用sym8作為小波基, 并采用5層小波分解進行去噪. 去噪前后結果如圖3所示. 小波變換公式為

(1)

圖3 影像去噪前后結果Fig.3 Results before and after denoising of image

1.5 研究方法

1.5.1波長篩選和處理算法

1) 連續(xù)投影算法(SPA)是一種降維算法, 從光譜變量中找到冗余信息最少的變量組合. 該算法通過將波長投影到其他波長上, 比較投影向量大小, 以投影向量最大的波長為待選波長, 然后基于矯正模型選擇最終的特征波長[14].

2) 皮爾遜相關系數(Pearson)可以反映兩個變量的線性相關程度, 本研究主要是篩選出能用于各水質參數反演波長和波長組合.

采用方差膨脹因子(VIF)主要是為了避免使用相似作用的因子, 減少模型的誤差. 當該值大于10時, 則說明多重共線性較強, 需要舍棄.

3) 大量研究表明, 對兩種及以上波長進行組合計算, 可以降低部分外界干擾, 增強有效的光譜信息, 從而增大模型的擬合度. 因此, 本研究構建了雙波長差異指數(ID)、 比率指數(IR)和歸一化差異指數(IND)模型對波長進行處理. 即

(2)

式中:i=1, 2, …, 176;j=1, 2, …, 176;Bi和Bj表示無人機遙感影像的波長值.

1.5.2機器學習算法和評價指標

1) 采用機器學習算法如隨機森林回歸(RF)、 支持向量回歸(SVR)、 KNN回歸(KNN)、 Adaboost回歸算法建模. 機器學習算法部分均采用Sklearn庫, 其中SVM采用的是徑向基核函數, 并采用GridSearchCV庫進行調參. AdaBoost的基本思想是通過多次迭代, 將若干個弱分類器組合成一個強分類器[15].

2) 評價指標為決定系數(R2)、 平均相對誤差(EMR)和均方根誤差(ERMS). 其中:R2值在0～1之間, 相當于分類問題的精度, 其值越高越好;EMR用來衡量估計值與真值的相對誤差百分比;ERMS用來衡量估計值與真值之間的偏差. 即

(3)

2 實驗結果與分析

2.1 相關性分析

1) 各水質間的相關性. 由于TP、 TN水質參數對光譜不敏感, 因此參照李瀾博士的處理方式, 計算TP、 TN與當天所采的其他9個水質參數的相關系數, 并將相關度較大的納入反演模型自變量中[16]. 相關系數熱力圖如圖4所示, 其中相關性大小用r表示. 由圖中可知, TP、 TN與氨氮的相關性為0.70、 0.61, 說明氨氮 指標能影響TP、 TN質量濃度, 因此本研究嘗試將影響氨氮的特征波長納入TP、 TN水質反演中[16].

2) TP、 TN波長組合相關性分析. 選取晉安河86個樣本點176波長, 逐波長遍歷構建ID、IR和IND指數, 并計算其與TP、 TN濃度的皮爾遜相關系數. 水質參數與雙波長指數相關性等值圖如圖5所示.

圖5 相關性等值線分布圖Fig.5 Correlation contour map

從圖5可以看出, 采用ID指數對兩兩波長進行處理, 發(fā)現其與TP相關性最高接近0.8, 波長組合主要分布在600、 750 nm附近; TN與雙波長ID指數相關性最高接近0.8, 相關性高的區(qū)域主要分布在800～900 nm之間, 相關性滿足河流水質模型反演要求.

2.2 自變量篩選

由于高光譜數據原始波長數據量大, 自變量個數較多, 其中無關變量也多, 從而導致模型計算量大, 穩(wěn)定性低. 因此采用合適的自變量是水質反演模型好壞的重中之重.

2.2.1SPA特征波長篩選

通過多次采用SPA算法設置不同閾值進行實驗, 實驗表明在樣本數為86時, SPA迭代波長次數最高設置為20次. 同時將各波長反射率歸一化(R′),R′范圍為[-1,1], 模型效果較優(yōu). 如圖6所示, 對總磷的波長進行篩選, SPA算法迭代到波長數為14時,RMSE處于最低值。最終總磷SPA選擇888、 885、 404、 763、 899、 866、 401、 813、 827、 870、 670、 820、 475、 795 nm這14個波長; 總氮SPA選擇763、 866、 728、 462、 813、 707、 401、 881、 578 nm這9個波長.

圖6 SPA對總磷的特征波長選擇Fig.6 SPA’s selection of characteristic bands for TP

2.2.2Pearson特征波長組合篩選

采用經驗法選用Pearson系數較高的5個指數值(包括Pearson系數最高的ID、IR、IND各1個, 其余2個為所有Pearson系數最高的雙波長指數)作為待反演水質模型的自變量, 同時采用VIF對這5個自變量進一步篩選, 剔除共線性強(VIF>10)的自變量,選出最終參與模型運算的雙波長組合. 采用如上方法對TP、 TN和氨氮進行計算, 得到TP的雙波長組合為ID(874, 867)、IR(750, 760)、ID(757, 760); TN雙波長組合為ID(853, 860)、IND(670, 663); 氨氮雙波長組合為ID(874, 867)、ID(767, 760)、IND(544, 600).

2.2.3最終采用的自變量

結合2.2.1與2.2.2, 晉安河TP、 TN反演最終采用的模型自變量如表2所示.

表2 各水質參數篩選的自變量

2.3 模型構建與精度評價

基于2.2.3所挑選的自變量作為模型訓練的輸入, 隨機挑選70%樣本點進行模型訓練, 30%樣本點進行測試. 采用10折交叉驗證建立回歸模型, 以決定系數(R2)、 平均相對誤差 (EMR)、 均方根誤差 (ERMS) 作為評價指標. 測試集對模型各指標檢驗結果如表3所示.

表3 各水質參數反演模型

由表3可知, 對于TP、 TN來說, 采用SPA+Pearson比SPA篩選算法模型的百分比偏差均有所降低,R2均有所提升, 這表明波長篩選算法結合的方式確實比單純采用單波長建模效果更優(yōu). 并且將影響氨氮波長組合加入TP、 TN反演時, TP模型R2從0.86提升到0.92,ERMS從0.008 mg·L-1降低至0.005 mg·L-1, TN模型R2從0.87提升到0.90,ERMS從0.086 mg·L-1降低至0.082 mg·L-1, 模型的效果有所提升.

測試集最優(yōu)模型驗證圖如圖7所示. 其中TP、 TN反演模型R2均為0.9以上, 實驗結果表明, 將SPA與Pearson波長篩選結合的方式, 同時考慮其他水質參數的影響, 可以取得較好的反演效果. 這對復雜內陸水體光學活性較差的水質參數反演有著借鑒意義.

圖7 晉安河最優(yōu)模型與相應水質參數驗證圖Fig.7 Verification diagram of Jin’an River optimal model and corresponding water quality parameters

2.4 模型的泛化能力

采用福州市長樂區(qū)下洞江, 對SPA+Pearson+機器學習算法模型構建的TP、 TN模型泛化能力進行驗證. 下洞江樣本點數量為87, 數據分布如表4所示, 取其中70%作為訓練集, 30%作為測試集. 實驗結果如圖8所示, 可以看出, 采用SPA+Pearson+氨氮的隨機森林算法構建的TP模型決定系數R2為0.86,ERMS為0.032 mg·L-1, 采用SPA+Pearson+氨氮的SVR算法構建的TN模型決定系數R2為0.74,ERMS為0.450 mg·L-1. 結合晉安河實例, 該方法對下洞江的水質反演效果沒有晉安河好, 其原因是下洞江水閘處于開放狀態(tài), 相較晉安河河流流速較快, 下洞江水質樣本數據相對較差. 但盡管如此, 其水質反演效果仍舊良好, 這說明利用本研究提出的算法模型泛化能力不錯, 能應用于其他河流水質反演之中.

表4 下洞江取樣點實測水質數據統(tǒng)計

圖8 下洞江最優(yōu)模型與相應水質真實值和預測值驗證圖Fig.8 Verification diagram of Xiadong River optimal model and real-predicted value

2.5 反演結果

通過上述模型之間的對比, 選擇最優(yōu)模型對東門村河段無人機高光譜影像進行TP、 TN參數的反演可視化. 可視化結果如圖9所示, 可以看出該河段TP質量濃度范圍為0.09～0.14 mg·L-1, TN質量濃度范圍為1.80～2.60 mg·L-1. 經現場調查發(fā)現, 河流兩側為居民地, 存在面源污染的情況, 其中左上角和左下角有兩個居民排污口, 因此河流左側相對橋下的TP、 TN濃度偏高, 從圖中也可以明顯看出這一趨勢, 研究結果符合實地調研情況.

圖9 東門村河段總磷和總氮反演圖Fig.9 Inversion map of TP and TN in Dongmen Village reach

3 結語

本研究基于晉安河實測水質數據和光譜數據, 通過SPA算法對光譜波長降維提取, 并采用經VIF調整后的雙波長Pearson相關性分析挑選出相關性強的波長組合, 之后再將與TP、 TN強相關的氨氮水質參數加入機器學習模型自變量中進行建模. 經測試集實測數據與模型預測結果的驗證分析, TP、 TN兩者最高模型R2均可達0.9,ERMS也相對較低, 同時在福州市下洞江對該方法模型進行驗證, 實驗表明模型擬合效果和泛化能力都較為出眾. 最后采用最佳模型對晉安河東門村河段進行了TP、 TN的水質反演可視化, 結果表明, 各模型對無人機高光譜影像水質參數濃度反演的空間分布與東門村河段的實地調研情況基本符合.

目前采用無人機遙感進行水環(huán)境治理前景很大, 但季節(jié)性差異和區(qū)域性差異對河流水質的影響較大, 未來可以研究不同季節(jié)之間的模型差異, 從而進一步增強模型的泛化能力.