改進SVR的內陸水體COD高光譜遙感反演

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

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

引 言

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

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

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

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

1 研究區域與數據

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

1.1 實測數據

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

圖1 濰河及峽山水庫現場實驗點位圖

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

1.2 遙感數據

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

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

2 改進SVR的內陸水體COD高光譜遙感反演模型

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

2.1 光譜數據與COD的相關性分析

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

(1)

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

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

(2)

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

2.2 基于SA-PSO算法優化的SVR模型

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

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

(3)

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

由于SVR模型存在參數選取難的問題，如何確定最優參數直接影響到SVR模型的學習和泛化能力。根據鳥群覓食行為提出的粒子群算法PSO常用于對SVR模型參數進行優化選取，將只有位置和速度兩個屬性的粒子模擬為鳥，每個粒子在空間中的極值Pbest作為潛在的最優解Gbest; 所有粒子根據式(4)不斷調整各自的位置和速度，直至獲得最優解。同時，在粒子群更新過程中引入模擬退火算法SA可提高PSO算法的全局尋優能力，即利用其在尋優過程中的突跳能力對粒子群算法進行改進。SA算法本質是模擬高溫物體退火過程中尋找全局最優解的過程。初始溫度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算法優化的SVR模型，關鍵是利用SA-PSO算法優化支持向量回歸模型中判別函數的懲罰因子C和核函數中核的寬度g。建立SA-PSO優化SVR模型的主要步驟如圖2所示。

圖2 SA-PSO優化SVR流程

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

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

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

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

(5)輸出最優參數C和g。

3 實驗與結果討論

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

3.1 光譜數據與COD的相關性分析

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

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

Pearson相關性分析的輸入反射率形式包括單波段和波段比值等形式。研究表明, 波段比值組合可部分消除水表面光滑度和微波等其他環境因素的影響[15]，能在一定程度上有效提高水質參數反演的精度。對各個采樣點的COD濃度分別與每組光譜數據的波段比值組合做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反演

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

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

圖5 SA-PSO-SVR精度分析

3.3 SA-PSO-SVR精度分析

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

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

圖6 SVR, BP神經網絡，LR模型預測評估

3.4 COD濃度空間分布分析

將SA-PSO-SVR模型應用于OHS數據，得到的COD濃度空戒分布如圖7所示。可以看出，COD的濃度呈現部分區域濃度高的特點，在韓信壩、峽山水庫的東北部區域、渠河與濰河的交叉口及輝村與金口大橋之間的濃度明顯偏高。

圖7 COD濃度分布情況

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

4 結 論

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

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

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

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