盛 輝，池海旭，許明明*，劉善偉，萬劍華，王錦錦

1.中國石油大學(華東)，海洋與空間信息學院，山東 青島 266580 2.珠海歐比特宇航科技股份有限公司，廣東 珠海 519080

引 言

河流及水庫作為內(nèi)陸水資源的重要組成部分，對涵養(yǎng)水源、保護生態(tài)具有重要的作用。但近幾十年來，自然因素和人為因素的雙重壓力使水環(huán)境在不同程度上受到了污染，水體富營養(yǎng)化的趨勢愈加嚴峻，對水質(zhì)進行動態(tài)監(jiān)測和精細化的管理迫在眉睫。

傳統(tǒng)的水質(zhì)監(jiān)測工作主要是獲得監(jiān)測區(qū)域內(nèi)某些點的水質(zhì)濃度信息，這些有限點的信息并不能反映出整個區(qū)域的情況，而遙感水質(zhì)監(jiān)測具有實時高效、監(jiān)測范圍廣和適合于長期動態(tài)監(jiān)測等特點[1]?；瘜W需氧量(chemical oxygen demand，COD)可以用來表征水中有機物的含量，是遙感水質(zhì)監(jiān)測的重要指標。COD遙感監(jiān)測主要是研究水體反射的光譜特征和COD濃度之間的關系，從而建立反演算法。水質(zhì)參數(shù)反演的方法主要有兩類：半經(jīng)驗方法和半分析方法[2]。一些學者利用半經(jīng)驗的分析方法研究COD和遙感參數(shù)之間的關系，通過構建的線性及生物光學等模型反演得到了內(nèi)陸湖泊COD含量的分布情況，并進一步分析得到了COD的來源及影響[3]。但半經(jīng)驗模型的擬合精度范圍大多在0.6～0.8之間，因此模型的準確性還有待提高。國內(nèi)外的部分學者研究評估機器學習進行水質(zhì)遙感的可行性[4-5]，表明在充足樣本的情況下，人工神經(jīng)網(wǎng)絡能夠對水質(zhì)參數(shù)進行有效的反演分析。由于某些機器學習算法參數(shù)選擇的不確定性以及COD光譜特征的微弱性，通常需要多次重復實驗來獲得最優(yōu)參數(shù)從而提高模型的準確性。另外，研究人員在野外實驗中需要花費大量的時間和精力才能夠收集到足夠的水質(zhì)參數(shù)數(shù)據(jù)。因此，構建結構簡單、適合小樣本的水質(zhì)參數(shù)反演模型尤為重要。研究表明，經(jīng)過優(yōu)化的支持向量回歸機(support vector regression，SVR)具有結構簡單、全局最優(yōu)的特點，已經(jīng)廣泛應用到葉綠素、懸浮物、溶解氧和水體富營養(yǎng)化的水質(zhì)研究中[6-7]，但是針對COD的反演研究較少。因此，本工作采用改進的SVR方法進行COD反演。

用于COD反演的數(shù)據(jù)源主要集中于Landsat，Modis和Sentinel-2等多光譜數(shù)據(jù)[8]。高光譜數(shù)據(jù)可以捕獲由不同水質(zhì)濃度引起的光譜變化，并在內(nèi)陸水質(zhì)監(jiān)測中顯示出巨大的潛力[9-10]。但由于水體光譜影響機理的復雜性和其他水質(zhì)參數(shù)的影響[11]，狹小區(qū)域的COD反演和變化監(jiān)測需要兼有高空間分辨率和高光譜分辨率的數(shù)據(jù)提供支持。實測的水表面光譜反射率具有受大氣影響較小的特征[12]，將其與高光譜圖像結合可以更準確地估算COD濃度[13]。林建遠等通過航空高光譜數(shù)據(jù)和實測的水表面反射率反演得到了杭嘉湖平原河道的COD含量，但是航空高光譜數(shù)據(jù)的獲取難度大、成本高。珠海一號高光譜衛(wèi)星(orbita hyper spectral, OHS)數(shù)據(jù)具有高空間分辨率、高時間分辨率和易于獲取的特點。因此，采用現(xiàn)場實測光譜反射率數(shù)據(jù)結合珠海一號高光譜衛(wèi)星數(shù)據(jù)用于COD濃度的反演。

針對COD光譜特征的微弱性以及SVR參數(shù)選取困難、易陷入局部極值的問題，基于高光譜數(shù)據(jù)，提出了一種模擬退火—粒子群算法(simulated annealing-particle swarm algorithm，SA-PSO)優(yōu)化的SVR方法進行COD濃度反演。首先，基于實測的光譜數(shù)據(jù)建立濰河—峽山水庫區(qū)域的COD高光譜反演模型，減少大氣校正對反演精度的影響。然后應用于OHS高光譜數(shù)據(jù)對該區(qū)域的COD含量進行遙感估算，以期為濰河流域的水質(zhì)監(jiān)測與綜合管理和OHS衛(wèi)星數(shù)據(jù)在內(nèi)陸水質(zhì)的反演應用提供基礎數(shù)據(jù)與科學參考。

1 研究區(qū)域與數(shù)據(jù)

濰河及峽山水庫位于山東半島，是膠東地區(qū)的戰(zhàn)略水源地之一。濰河是流經(jīng)峽山水庫的主要河流之一，發(fā)源于沂水縣和莒縣，最終注入渤海，河長200多公里，流域面積近萬平方公里；峽山水庫的注入河流有濰河、渠河和浯河等河流，入庫量為8×108m3左右。峽山水庫周圍及濰河上下游區(qū)域分布著眾多村莊和工廠，該河域為周邊地區(qū)的居民飲水，農(nóng)業(yè)灌溉和工業(yè)發(fā)展做出了重要貢獻。同時，該地區(qū)的百姓生活與經(jīng)濟發(fā)展也對水質(zhì)造成了一定的影響，采用遙感技術對該地區(qū)的水質(zhì)進行監(jiān)測顯得尤為必要。

1.1 實測數(shù)據(jù)

2019年10月26日—10月28日，分別在濰河流域的古縣大橋、峽山水庫、輝村和金口大橋附近水域進行水體取樣和光譜采集，采集光譜數(shù)據(jù)12組，28日有效測得COD濃度數(shù)據(jù) 22組。將采水器采集到的表面水樣裝在琥珀瓶中，在瓶體貼標簽記錄好采樣地點、時間和當時的天氣情況并拍照保存，最終在實驗室化驗得到各采樣點COD的濃度值。

圖1 濰河及峽山水庫現(xiàn)場實驗點位圖

光譜測量采用YW-TRIOS-AWRMMS水面移動測量系統(tǒng)(光譜范圍320～943.7 nm, 共190個波段)，該儀器內(nèi)置GPS定位系統(tǒng)，攜帶方便可用于各種情況的光譜測量。測量過程中，保持伸縮桿方向和太陽入射平面135°的夾角可使3根傳感器分別有效測量太陽輻照度、水面輻照度和天空輻照度三種參數(shù)，并可在顯示屏上查看光譜曲線。在每個采樣點測量10組光譜以最小化不確定性，通過光譜數(shù)據(jù)的處理程序及Mobley(1999)的方法對光譜數(shù)據(jù)進行處理，即可得到離水輻亮度和遙感反射率等光譜數(shù)據(jù)信息參數(shù)。

1.2 遙感數(shù)據(jù)

OHS衛(wèi)星于2018年4月26日成功發(fā)射升空，彌補了我國在高光譜數(shù)據(jù)上的不足，開啟了商業(yè)航天遙感的新時代。OHS衛(wèi)星空間分辨率為10 m, 4顆高光譜衛(wèi)星的重訪周期為2.5 d，單次成像范圍為150 km×320 km，光譜分辨率為3～8 nm，光譜范圍為466～940 nm, 共32個波段。由于其具有幅寬大、高空間分辨率、高光譜分辨率和高時間分辨率的特點，非常有利于反演光學特性復雜多變的內(nèi)陸水體、小型水庫和河流的水質(zhì)參數(shù)等相關研究。OHS數(shù)據(jù)是當前能夠進行小區(qū)域水質(zhì)參數(shù)反演的最佳數(shù)據(jù)源。

采用珠海歐比特宇航科技股份有限公司提供的2景濰河流域2019年10月28日的OHS高光譜數(shù)據(jù)。首先對影像數(shù)據(jù)進行輻射定標、大氣校正、影像拼接與陸地掩膜等處理，最終得到僅有濰河流域區(qū)域在內(nèi)的影像產(chǎn)品。大氣校正是為了消除大氣分子和氣溶膠的影響，是準確獲得水質(zhì)遙感信息的前提。FLAASH模型用于遙感影像的大氣糾正具有很好的應用效果；故使用FLAASH模型進行大氣校正，然后根據(jù)采樣點的坐標位置，在影像上提取各點對應的反射率信息。

2 改進SVR的內(nèi)陸水體COD高光譜遙感反演模型

建立COD的反演模型，首先要對光譜數(shù)據(jù)進行歸一化和相關性分析等預處理，然后基于對COD敏感的波段結合SA-PSO算法建立SVR反演模型。

2.1 光譜數(shù)據(jù)與COD的相關性分析

為了削弱野外環(huán)境給光譜測量造成的影響，同時也使光譜數(shù)據(jù)和OHS影像光譜范圍相對應，采用歸一化的方法對466～940 nm的光譜數(shù)據(jù)進行預處理。歸一化的公式如式(1)

(1)

式(1)中：RN(λi)為水體遙感反射率的歸一化結果，R(λi)為原始的遙感反射率，λi為i處波長，n為466～940 nm處的波段數(shù)。

Pearson相關性分析是從統(tǒng)計學的角度研究兩個或多個隨機變量間關聯(lián)度強弱的方法。相關系數(shù)的大小可以描述變量間的密切程度，兩個變量的相關系數(shù)表達式如式(2)

(2)

式(2)中：yi是各個采樣點水質(zhì)參數(shù)的濃度，xi是各個采樣點的地表反射率。

2.2 基于SA-PSO算法優(yōu)化的SVR模型

SVR是基于機器學習的一種算法，通過核函數(shù)將非線性數(shù)據(jù)映射到高維空間構造決策函數(shù)進行線性回歸，在解決小樣本、非線性和高維模式識別問題方面具有其獨特的優(yōu)勢，經(jīng)常用于小樣本水質(zhì)參數(shù)的反演[14]。SVR模型的數(shù)學關系為式(3)

f(xi)=ωTφ(xi)+B

(3)

式(3)中，xi為樣本數(shù)據(jù)，ω為待辨識的權重，φ(xi)為非線性映射，f(xi)為特征空間中的線性函數(shù)，B為常數(shù)項。

由于SVR模型存在參數(shù)選取難的問題，如何確定最優(yōu)參數(shù)直接影響到SVR模型的學習和泛化能力。根據(jù)鳥群覓食行為提出的粒子群算法PSO常用于對SVR模型參數(shù)進行優(yōu)化選取，將只有位置和速度兩個屬性的粒子模擬為鳥，每個粒子在空間中的極值Pbest作為潛在的最優(yōu)解Gbest; 所有粒子根據(jù)式(4)不斷調(diào)整各自的位置和速度，直至獲得最優(yōu)解。同時，在粒子群更新過程中引入模擬退火算法SA可提高PSO算法的全局尋優(yōu)能力，即利用其在尋優(yōu)過程中的突跳能力對粒子群算法進行改進。SA算法本質(zhì)是模擬高溫物體退火過程中尋找全局最優(yōu)解的過程。初始溫度T的確定采用適應度和和接受概率的方法，由式(5)決定

Vi+1=ω×Vi+c1×rand×(Pbest-Xi)+

c2×rand×(Gbest-Xi)

Xi+1=Xi+Vi+1

(4)

T=(fmax-fmin)/lnp=-|Δf|/lnp

(5)

式中：fmax和fmin和Δf為初始粒子群最大、最小適應度值及其差值，p為初始接受概率，p=0.8。

建立基于SA-PSO算法優(yōu)化的SVR模型，關鍵是利用SA-PSO算法優(yōu)化支持向量回歸模型中判別函數(shù)的懲罰因子C和核函數(shù)中核的寬度g。建立SA-PSO優(yōu)化SVR模型的主要步驟如圖2所示。

圖2 SA-PSO優(yōu)化SVR流程

(1)粒子群的初始化。對粒子速度，位置，懲罰因子C和RBF核函數(shù)寬度等參數(shù)進行初始化。設置初始退火溫度T，并將最大迭代次數(shù)設置為200。

(2)計算群體中每個粒子的適應度。如果新的適應性值好于原始值，則接受新位置。否則，將保留舊位置。

(3)根據(jù)式(5)更新粒子的位置和速度。

(4)收斂性的判斷。如果滿足終止條件，則算法停止。否則，執(zhí)行退火操作，轉到(2)。

(5)輸出最優(yōu)參數(shù)C和g。

3 實驗與結果討論

進行COD的高光譜反演，首先要對光譜進行歸一化和相關性分析來確定敏感因子，然后根據(jù)敏感因子和COD濃度建立反演模型。最后通過使用與采樣點相對應的OHS反射率數(shù)據(jù)來分析模型的準確性，并將SA-PSO-SVR應用于濰河流域的OHS數(shù)據(jù)進行COD的空間分析。

3.1 光譜數(shù)據(jù)與COD的相關性分析

實測的原始光譜曲線如圖3(a)所示，根據(jù)等式(1)進行歸一化的結果如圖3(b)所示。與原始光譜曲線相比，歸一化光譜曲線的形狀發(fā)生了變化，反射谷和反射峰更加明顯。從560，680和710 nm的反射率值可以發(fā)現(xiàn)，當濃度增加時，反射峰具有向短波長方向移動而反射谷向長波方向移動的趨勢。

圖3 原始光譜(a)和歸一化后的光譜(b)

Pearson相關性分析的輸入反射率形式包括單波段和波段比值等形式。研究表明, 波段比值組合可部分消除水表面光滑度和微波等其他環(huán)境因素的影響[15]，能在一定程度上有效提高水質(zhì)參數(shù)反演的精度。對各個采樣點的COD濃度分別與每組光譜數(shù)據(jù)的波段比值組合做Pearson相關性分析，得到COD濃度與波段比值組合的相關性分析結果如圖4所示，最佳反演因子是518 nm/940.4 nm，663.6 nm/636.8 nm，729.2 nm/890.9 nm和752.3 nm/857.9 nm四種波段比值組合。

圖4 466～940 nm的實測光譜比值與COD濃度的Pearson相關性分析結果

3.2 COD反演

根據(jù)選取的最佳反演因子，用12個采樣點的光譜數(shù)據(jù)來訓練SA-PSO-SVR模型。根據(jù)相關研究及反復實驗，設置種群規(guī)模為50，進化次數(shù)為200，選取參數(shù)c1和c2均為0.5，SVR模型的C和g的取值范圍為[0.1, 1 000]和[0.1, 10]。對訓練數(shù)據(jù)進行歸一化處理后，經(jīng)過SA-PSO優(yōu)化獲得的懲罰因子C與核參數(shù)g分別為151.09和0.36。

利用28日測得的22個采樣點數(shù)據(jù)作為檢查點，在OHS影像上提取與采樣點相對應的反射率，導入SA-PSO-SVR模型得到各個檢查點的預測值。圖5為SA-PSO-SVR精度分析結果。紅色點為利用光譜數(shù)據(jù)建立模型得到的實測值和預測值的分布情況，雖然某些采樣點的預測值與實測值具有一定差異，但依然可以看出SA-PSO-SVR模型的反演結果與實測值呈現(xiàn)較好的相關一致性，其模型決定系數(shù)為0.86。綠色點為檢查點數(shù)據(jù)的預測值和實測值的分布情況，平均相對誤差(MRE)和均方誤差(RMSE)分別為9.04%和3.64 mg·L-1，說明該模型可以實現(xiàn)對COD濃度的有效反演。

圖5 SA-PSO-SVR精度分析

3.3 SA-PSO-SVR精度分析

從圖5可以看出，OHS數(shù)據(jù)部分反演結果與實測值差距略大，推測原因為數(shù)據(jù)質(zhì)量參差不齊。由于TRIOS光譜測量儀和OHS衛(wèi)星是兩種截然不同的傳感器，傳感器的信噪比等參數(shù)并不一致，因此兩種儀器獲得的數(shù)據(jù)質(zhì)量并不一致。研究表明，COD對水體的光學特性響應很微弱，為了減少河道兩邊其他非水體像元的影響，采樣點的像元位置盡可能的選在距離監(jiān)測點近的河道中央，但由于河道中央流速較快，可能導致河中央同采樣點的COD濃度存在一定差異，進而造成部分反演結果與實測值的差異較大。

另外，基于光譜數(shù)據(jù)分別建立SVR和BP神經(jīng)網(wǎng)絡和線性回歸(LR)模型，將SVR模型的懲罰因子C與核參數(shù)g設置為109.85和0.000 1；根據(jù)相關研究及反復實驗，設置BP神經(jīng)網(wǎng)絡模型的隱含層節(jié)點為12，最大訓練次數(shù)為2 000，期望誤差為0.000 1，學習速率為0.01。SVR和BP神經(jīng)網(wǎng)絡和LR模型的決定系數(shù)分別為0.712，0.598和0.494。從圖6(a)可以看出，SVR模型的整體反演效果較好，但是OHS數(shù)據(jù)的部分反演結果低于實際測量值。從圖6(b)可以看出BP神經(jīng)網(wǎng)絡模型的RMSE和MRE分別為1.93 mg·L-1和5.18%，與其他兩種模型的結果相差不大；但是從OHS數(shù)據(jù)的反演情況看，BP神經(jīng)網(wǎng)絡模型的結果陷入了局部最優(yōu)，由此也證明BP神經(jīng)網(wǎng)絡模型確實存在參數(shù)選取困難的問題。另外，從圖6(c)可以看出，光譜數(shù)據(jù)的RMSE和MRE分別為3.00 mg·L-1和8.53%，與其他模型結果相近，但OHS影像的反演結果要普遍高于實測值。綜合來看，SA-PSO-SVR模型的反演效果優(yōu)于其他模型。

圖6 SVR, BP神經(jīng)網(wǎng)絡，LR模型預測評估

3.4 COD濃度空間分布分析

將SA-PSO-SVR模型應用于OHS數(shù)據(jù)，得到的COD濃度空戒分布如圖7所示?？梢钥闯?，COD的濃度呈現(xiàn)部分區(qū)域濃度高的特點，在韓信壩、峽山水庫的東北部區(qū)域、渠河與濰河的交叉口及輝村與金口大橋之間的濃度明顯偏高。

圖7 COD濃度分布情況

結合濰河流域地理情況并查閱相關資料發(fā)現(xiàn)，濰河流域附近分布著眾多工廠，大部分企業(yè)將工業(yè)污水排入河中，這些河流(濰河、渠河等)攜帶污染物最終匯入峽山水庫，是導致峽山水庫東部沿岸區(qū)域的COD濃度明顯高于其他區(qū)域的原因之一。另外，庫區(qū)及濰河等其他河流兩岸分布著眾多村莊，百姓的農(nóng)業(yè)生產(chǎn)活動依賴于該區(qū)域的水資源，農(nóng)業(yè)生產(chǎn)帶來的化肥、農(nóng)藥等污染物也是造成該局部區(qū)域COD濃度過高的原因之一。同時，研究表明農(nóng)作物腐爛死亡降解產(chǎn)生的有機物也會造成局部區(qū)域的COD濃度增高。另外，在韓信壩和山陽村附近區(qū)域有多處攔河大壩，部分區(qū)域河流徑流量減少，遙感影像獲取的COD反射信息不充分，造成該區(qū)域的COD反演結果偏高。

4 結 論

針對SVR模型參數(shù)選取困難和易陷入局部極值的情況，本文基于水面以上的光譜數(shù)據(jù)建立了SA-PSO-SVR模型，通過珠海一號衛(wèi)星數(shù)據(jù)反演得到了濰河—峽山水庫流域COD的分布情況，得到以下結論:

(1)通過對濰河—峽山水庫區(qū)域的實測光譜進行分析，結果表明該區(qū)域實測的水面光譜具有典型內(nèi)陸湖泊水域復雜渾濁水體的光譜特征，560和710 nm附近的光譜曲線形狀呈現(xiàn)明顯的雙峰特征，反射率幅值較大。當濃度增加時，反射峰具有向短波長方向移動而反射谷向長波長方向移動的趨勢。

(2)引入SA-PSO算法對SVR模型的參數(shù)進行優(yōu)化選取，解決了局部最優(yōu)的問題，其反演效果明顯也好于其他模型。將基于實測光譜建立的SA-PSO-SVR模型應用在珠海一號衛(wèi)星數(shù)據(jù)上可以進行COD遙感估算，說明了基于實測光譜和高光譜衛(wèi)星影像的水質(zhì)參數(shù)反演方法具有良好的應用前景和推廣價值。同時，反演得到的COD濃度分布情況可以為珠海一號衛(wèi)星數(shù)據(jù)的內(nèi)陸水質(zhì)參數(shù)反演和濰河流域的綜合管理提供科學依據(jù)。

珠海一號數(shù)據(jù)的驗證結果說明了COD反演的有效性，但由于星載傳感器和地面光譜儀信噪比等儀器參數(shù)的不同，并受到其他復雜的大氣和光學水體組分的影響，該方法還需要有足夠的論據(jù)對此進行進一步的詳細分析。此外，雖然基于實測光譜建立的SA-PSO-SVR模型在濰河流域表現(xiàn)出了相對較好的性能，但在其他復雜的水環(huán)境的應用性能還有待研究。因此，接下來的工作將多次收集濰河流域的光譜、COD濃度等其他水體成分的信息，并根據(jù)長期觀測資料進一步探討COD的時空變化；同時，在包含多個內(nèi)陸水域的綜合數(shù)據(jù)集的基礎上，對SA-PSO-SVR方法進行更全面的評估。