基于PLSR的珠江口城市河流水質(zhì)高光譜反演

黃 華，李茂億，陳吟暉，陳 耿，劉海龍，邢前國，蔡建楠

(1.廣東環(huán)境保護工程職業(yè)學(xué)院，廣東 佛山 528216；2.中山市環(huán)境監(jiān)測站，廣東 中山 528403； 3.廣東省環(huán)境監(jiān)測中心，廣東 廣州 510308；4.中國科學(xué)院煙臺海岸帶研究所，山東 煙臺 264003)

珠江口地區(qū)河網(wǎng)密布，河流水質(zhì)狀況直接影響城市的整體風(fēng)貌和居民生活質(zhì)量[1]。近年來隨著國家《水污染防治行動計劃》的深入推進，城市河流的綜合整治和水質(zhì)改善已成為珠江口各市的重要工作任務(wù)。為實現(xiàn)對河流水質(zhì)狀況監(jiān)控及河流整治效果的評估，往往需要配套開展密集、長期的水質(zhì)監(jiān)測工作，而采用常規(guī)的水樣現(xiàn)場采集和實驗室分析的方法存在耗時長、花費高、采樣點位有限等局限，無法滿足新形勢下水環(huán)境管理工作的需要[2]。

水質(zhì)遙感技術(shù)基于水體水質(zhì)與光譜特征之間的響應(yīng)關(guān)系實現(xiàn)對水質(zhì)指標的反演[3]，獲得的水質(zhì)信息具有時間和空間上的相對連續(xù)性，可以作為傳統(tǒng)水質(zhì)監(jiān)測方法的有效補充，近年來得到迅速發(fā)展。其中，高光譜水質(zhì)遙感技術(shù)由于具有納米級的光譜分辨率和連續(xù)的光譜譜段，極大地增加了水體光譜信息量，有助于精準建立各種水質(zhì)參數(shù)與光譜響應(yīng)關(guān)系以提高反演模型的精度，成為水質(zhì)遙感研究領(lǐng)域的熱點。目前，相關(guān)研究主要面向的監(jiān)測對象為水庫和海洋水體，而珠江口河網(wǎng)水體受陸源物質(zhì)大量輸入和海洋潮汐因素的共同作用，有著與內(nèi)陸二類水體和開闊大洋水體不同的水體光學(xué)特征，開展珠江口水體水質(zhì)光學(xué)遙感研究具有重要的意義[4]。目前城市河流水質(zhì)遙感有關(guān)研究中的目標水質(zhì)參數(shù)主要為水體濁度、懸浮物、葉綠素a、有色可溶性有機物(CDOM)等光學(xué)活性參數(shù)[5-7]，對于河流水質(zhì)綜合指標的遙感反演研究較少。

本文以珠江口城市中山市為研究區(qū)，通過開展河流高光譜監(jiān)測和同步水質(zhì)分析，研究受測水體高光譜特征與城市河流水質(zhì)指數(shù)(city river water quality index，CWQI)的關(guān)系，并基于偏最小二乘回歸(partial least squares regression， PLSR)建立高光譜數(shù)據(jù)與河流水質(zhì)指數(shù)的反演模型，同時探索水體高光譜監(jiān)測的最佳光譜分辨率，以期為城市河流水質(zhì)光譜監(jiān)測及有關(guān)高光譜傳感器的設(shè)計和應(yīng)用提供參考。

1 研究方法

中山市地處珠江口西岸，北接廣州，與深圳隔海相望，緊鄰港澳，處在粵港澳大灣區(qū)的地理中心位置。中山市境內(nèi)河流密布，不同規(guī)模的內(nèi)河涌、排洪渠高達1 000余條。改革開放以來，中山市經(jīng)濟迅猛發(fā)展，但隨著人口和產(chǎn)業(yè)聚集程度不斷提高，河流水體受到了一定程度的污染。

1.1 河流光譜監(jiān)測和數(shù)據(jù)處理

1.1.1地面高光譜監(jiān)測

于2019年7—9月選擇中山市境內(nèi)具有代表性的河流水體進行水質(zhì)高光譜監(jiān)測，共獲得151個點位的現(xiàn)場實測數(shù)據(jù)。高光譜監(jiān)測使用美國Ocean Optics公司USB4000光譜儀，光譜分辨率約為 0.20 nm，光譜范圍為345.34～1 046.12 nm，視場角為15°。每次測量前均使用標準白板進行儀器定標。測量時將光譜儀探頭垂直置于離河流水面 0.5 m 處測量河流水體上行輻射亮度Lu、天空輻射亮度Ls和太陽下行輻照度Ed，每個點位重復(fù)進行3次測量后取平均值作為測定結(jié)果。測點的光譜反射率Rrs計算公式[8]為

(1)

式中γ為天空光在水體表面的菲涅爾反射率，通常取經(jīng)驗值0.028?？紤]到400～900 nm波段范圍以外的光譜反射率更易受外界條件及儀器設(shè)計誤差等因素的影響，可能對監(jiān)測結(jié)果帶來信號噪聲[9]，因此僅選擇400～900 nm范圍內(nèi)共2 534個波段的原始光譜反射率數(shù)據(jù)用于構(gòu)建反演模型。

1.1.2無人機載高光譜監(jiān)測

無人機遙感監(jiān)測技術(shù)具有靈活性大、影像分辨率高、時效性強、成本低等優(yōu)點，是水環(huán)境遙感研究的一個重要發(fā)展方向[10]。在利用地面實測高光譜構(gòu)建CQWI反演模型后，將其推廣應(yīng)用于無人機高光譜影像中，以評價模型應(yīng)用于無人機載高光譜監(jiān)測的適用性。無人機載高光譜監(jiān)測采用北京智科遠達數(shù)據(jù)技術(shù)有限公司ZK-VNIR-FPG480型高光譜成像掃描儀，其有效譜段范圍為400～1 000 nm，光譜通道數(shù)大于或等于270，光譜分辨率小于或等于 2.8 nm，空間分辨率為9 cm(100 m飛行高度)，視場角為26°。高光譜成像掃描儀搭載在大疆M600 PRO無人機平臺上，在晴朗天氣條件下，根據(jù)設(shè)定的航飛路線采集研究水體的光譜，并同步采集靶標的光譜數(shù)據(jù)，對水體數(shù)據(jù)進行定標運算得到水體反射率數(shù)據(jù)。無人機航拍區(qū)域主要為內(nèi)河與外江交匯處、主要水工設(shè)施周邊、碼頭作業(yè)區(qū)等潛在的水質(zhì)突變點區(qū)域，以盡可能同時獲取不同水質(zhì)水體的光譜信息。

1.2 CWQI構(gòu)建和水質(zhì)監(jiān)測

目前，國內(nèi)外對于地表水環(huán)境質(zhì)量評價方法主要有單因子水質(zhì)類別評價法、水質(zhì)指數(shù)法和不確定性方法三大類[11]，其中水質(zhì)指數(shù)法可綜合多項水質(zhì)指標的監(jiān)測結(jié)果，計算過程簡單，評價結(jié)果便于比較。本文采用原環(huán)境保護部《城市地表水環(huán)境質(zhì)量排名技術(shù)規(guī)定(試行)》中規(guī)定的河流水質(zhì)指數(shù)作為反映河流水質(zhì)狀況的綜合指標，單項指標(除溶解氧外)和溶解氧的CWQI計算公式分別為

(2)

(3)

式中：Ii、IDO分別為水質(zhì)指標i和溶解氧的CWQI；ρi、ρ(DO)分別為水質(zhì)指標i和溶解氧的質(zhì)量濃度測定值；ρsi、ρs(DO)分別為水質(zhì)指標i和溶解氧Ⅲ類地表水對應(yīng)的質(zhì)量濃度標準限值。所有水質(zhì)指標的CWQI值相加即為該水樣的CWQI。

《城市地表水環(huán)境質(zhì)量排名技術(shù)規(guī)定(試行)》推薦CWQI計算的基本水質(zhì)指標為GB 3838—2002《地表水環(huán)境質(zhì)量標準》表1中除水溫、糞大腸菌群和總氮以外的21項指標。為優(yōu)化監(jiān)測項目，對中山市各河流監(jiān)測斷面2018—2019年的數(shù)據(jù)(共計水質(zhì)監(jiān)測樣本2 669個)進行了分析，結(jié)果顯示中山市各河流的主要超標因子為化學(xué)需氧量(COD)、總磷(TP)、氨氮(NH3-N)和溶解氧(DO)，這4項水質(zhì)監(jiān)測指標計算得到的CWQI值(I4)與21項水質(zhì)監(jiān)測指標計算得到的CWQI值(I21)之比的均值和標準差分別為53.7%和0.08%，相關(guān)系數(shù)r=0.94(P<0.05)，兩者具有顯著的相關(guān)性(圖1)，表明I4可較好地反映河流水質(zhì)狀況，故本文中選取上述4項指標作為與光譜測量同步進行的水質(zhì)監(jiān)測指標，并計算I4值作為河流水質(zhì)綜合評價結(jié)果。

表1 樣品水質(zhì)特征

圖1 4項和21項水質(zhì)監(jiān)測指標CWQI值的關(guān)系

各河流監(jiān)測點的水樣采集與地面高光譜測量同步進行，水質(zhì)監(jiān)測指標為COD、NH3-N、TP和DO。其中DO的質(zhì)量濃度使用美國YSI公司Pro 2030多參數(shù)水質(zhì)分析儀現(xiàn)場測定，其他指標取河流表層 0～20 cm處混合樣，裝入潔凈聚乙烯塑料瓶中并加入相應(yīng)保存劑后置于低溫保存箱及時送至實驗室分析。COD質(zhì)量濃度采用重鉻酸鉀法測定，NH3-N質(zhì)量濃度采用納氏試劑光度法測定，TP質(zhì)量濃度采用鉬銻抗分光光度法測定[12]，并使用空白樣、平行樣、標準樣品和加標回收等措施對樣品監(jiān)測全過程進行了質(zhì)量控制。水質(zhì)樣品各特征指標見表1，可見多數(shù)樣品水質(zhì)屬于Ⅳ或Ⅴ類，主要定類水質(zhì)指標為TP和NH3-N。

1.3 PLSR反演模型建立與評價

由于水體反射率數(shù)據(jù)波段數(shù)眾多，相鄰波段間隔緊密，而不同波段光譜反射率數(shù)據(jù)間最大方差膨脹因子大于80 000，表明不同波段光譜反射率存在極嚴重的多重相關(guān)性[13]。為提高反演模型精度，采用PLSR建立水體不同波段光譜反射率與CWQI的反演模型。PLSR基于高維投影思想，通過提取主成分的方法對系統(tǒng)中的數(shù)據(jù)信息進行分解和篩選，實現(xiàn)對系統(tǒng)解釋性最強的綜合變量的有效抽取[13-15]；PLSR常用于樣本量小、變量間存在多重相關(guān)性的數(shù)據(jù)處理，適用于中高光譜遙感數(shù)據(jù)反演模型的建立[16-17]。

PLSR建模的基本步驟[18]為：①標準化處理水體光譜反射率數(shù)據(jù)矩陣X和水體CWQI數(shù)據(jù)矩陣Y，得到的標準化矩陣分別記為E0和F0；②分別提取E0和F0第一對成分T1和U1，T1和U1分別為E0和F0的線性組合并使之相關(guān)性最大；③分別求E0和F0在T1上的回歸，并得到回歸殘差矩陣E1和F1；④用殘差矩陣E1和F1代替E0和F0重復(fù)上述步驟；⑤檢查收斂性以確定提取的主成分數(shù)。

將151個樣本原始數(shù)據(jù)隨機劃分為建模集(120個樣本)和驗證集(31個樣本)兩部分，其中建模集用建立PLSR反演模型；驗證集用于反演模型對CWQI反演效果的驗證。模型驗證效果分別采用相關(guān)系數(shù)(r)、均方根誤差(RMSE)、平均相對誤差(MRE)進行評價。

2 結(jié)果與分析

2.1 光譜反射率與CWQI的關(guān)系

根據(jù)CWQI計算結(jié)果，151個水質(zhì)樣品的CWQI值范圍為1.36～10.75，平均值為3.60，標準差為1.91。將不同CWQI值范圍的水樣進行分類并計算其高光譜反射率均值，結(jié)果(圖2)顯示，在400～700 nm波段范圍內(nèi)光譜反射率對不同CWQI值的水樣有較好的區(qū)分度，水體光譜反射率均值隨CWQI的增大依次遞減，其中I4≤2的水體反射率均值在整個波段范圍內(nèi)均高于其他類別的水體。但在 720～900 nm波段范圍內(nèi)，出現(xiàn)I4>5的較差水質(zhì)水體的光譜分辨率高于2

圖2 不同CWQI范圍水樣對應(yīng)的高光譜反射率均值

2.2 不同光譜分辨率對PLSR反演模型的影響

為分析光譜分辨率對水體PLSR反演模型估算效果的影響，需找出水體CWQI高光譜遙感監(jiān)測最佳光譜分辨率。采用分段三次Hermite插值多項式法對原始數(shù)據(jù)進行重采樣，采樣間隔依次為約 0.20 nm、1 nm、10 nm、20 nm、30 nm、40 nm、50 nm、60 nm、70 nm和80 nm，連同兩端波段的總波段數(shù)分別為2 534、900、51、26、18、14、11、10、9和8。在此基礎(chǔ)上建立不同光譜分辨率樣本集的PLSR反演模型，并根據(jù)相應(yīng)的建模與驗證精度評價結(jié)果研究模型估算精度隨光譜分辨率的變化。結(jié)果(表2)表明，在各光譜分辨率的數(shù)據(jù)集中，原始光譜數(shù)據(jù)集的PLSR反演模型的估算效果最差；對原始數(shù)據(jù)進行不同光譜分辨率的重采樣后，均有助于提升PLSR反演模型的估算效果。PLSR反演模型的驗證集RMSE隨著光譜分辨率的降低呈現(xiàn)先減小后增大的趨勢。光譜分辨率為50 nm、波段數(shù)為11個時，驗證集RMSE最小達0.790，MRE為18.63%，模型估算效果最優(yōu)。但總體上重采樣后各驗證集的PLSR反演模型估算效果均處于接近的水平。該現(xiàn)象表明，在CWQI的遙感監(jiān)測中，過高的光譜分辨率可能會造成數(shù)據(jù)冗余和引入環(huán)境噪聲，造成模型復(fù)雜化并影響模型穩(wěn)定性。反之，在較低的光譜分辨率下，數(shù)量較少的波段可在基本涵蓋水體主要污染物特征譜段的同時，增強光譜信噪比，從而提高估算精度。

表2 不同波段數(shù)的PLSR建模效果評價

2.3 提取的主成分數(shù)對PLSR反演模型的影響

確定提取的主成分數(shù)是PLSR反演模型建立和優(yōu)化的另一重要環(huán)節(jié)。目前PLSR建模確定主成分數(shù)的方法主要有舍一交叉驗證方法、分批交叉驗證方法、分裂樣本交叉驗證方法、隨機樣本交叉驗證方法等[17]，其中舍一交叉驗證方法的應(yīng)用較為廣泛，但該方法得到的建模主成分數(shù)可能不是最優(yōu)值，主成分數(shù)過多或過少可能會導(dǎo)致欠擬合或過擬合，影響模型的估算準確度[23]。本文采用50 nm分辨率的水體光譜反射率數(shù)據(jù)，分別提取2～10個主成分進行PLSR建模，以驗證集RMSE最小時對應(yīng)的主成分數(shù)為最佳的主成分數(shù)。圖3表明，模型的RMSE隨主成分數(shù)的增加呈先下降后增加的趨勢，當(dāng)主成分數(shù)為8時，模型的RMSE達到最低值。

圖3 不同主成分數(shù)與驗證集RMSE的關(guān)系

因此，確定建模最佳主成分數(shù)為8，并建立PLSR反演模型如下：

I=4.84-23.33R400+22.62R450+53.53R500+

4.69R550-24.19R600-102.32R650+

75.96R700+25.38R750-30.93R800+

46.15R850-19.99R900

(4)

式中Rj為不同波長處的水體光譜反射率。

2.4 PLSR反演模型估算效果評價

以31個驗證集樣本的CWQI實測值為縱坐標、CWQI高光譜反演模型(式(4))計算值為橫坐標，繪制其實測值與計算值的對比散點圖(圖4)。結(jié)果表明，該反演模型的估算精度良好，RMSE、MRE分別為0.768和18.1%，說明該反演模型在有效縮減了高光譜輸入?yún)?shù)的前提下，仍然保持了良好的泛化能力和反演精度。

圖4 驗證集樣本CWQI實測值與反演模型計算值比較

2.5 反演結(jié)果

由于無人機航拍覆蓋范圍較大，本文未進行同步的水質(zhì)分析，主要對基于無人機高光譜影像的CWQI反演效果做初步定性分析，并重點關(guān)注河流水質(zhì)突變區(qū)域CWQI反演結(jié)果的合理性。從應(yīng)用結(jié)果(圖5)上看，反演得到的CWQI值可較合理地反映出河流水質(zhì)的空間差異。如圖5(a)反演結(jié)果突出了地理位置接近，但因水閘控制導(dǎo)致連通性差的水體水質(zhì)的顯著差異；圖5(b)反演結(jié)果反映了河流交匯處水體混合過程中的局部水質(zhì)差異；圖5(c)反演結(jié)果突出了碼頭船舶作業(yè)引起的水體擾動和水質(zhì)變化。該結(jié)果表明，基于地面高光譜測點數(shù)據(jù)建立的CWQI反演模型應(yīng)用于空間影像數(shù)據(jù)反演CWQI具有較好的可行性。

(a) 某水閘兩側(cè)水體

受限于無人機航程的限制，每次作業(yè)時間約在20 min左右，對于河流水質(zhì)大范圍監(jiān)測的能力較為不足，僅能重點開展局部區(qū)域的水質(zhì)高光譜反演。同時，現(xiàn)場航飛采用多航帶作業(yè)方式，但對于水體流動性強，漲落潮劇烈的水域，水體流動導(dǎo)致水質(zhì)的時空變化不能被無人機同步采集；這種準同步的水體光譜數(shù)據(jù)與實際河流的水質(zhì)分布會有一定差異。但整體上反演模型應(yīng)用于無人機高光譜影像的CWQI反演可取得良好效果，表明該方法具有良好的業(yè)務(wù)化運行潛力。

3 結(jié) 語

2019年7—9月在珠三角地區(qū)中山市選擇代表性河流進行了水體高光譜反射率監(jiān)測和同步水質(zhì)監(jiān)測，并以基于COD、TP、NH3-N和DO 4項水質(zhì)指標計算得到的CWQI來表征河流水質(zhì)狀況。CWQI與水體高光譜反射率的相關(guān)分析表明，不同波段高光譜反射率對水質(zhì)狀況有較好的光學(xué)敏感性；基于PLSR建立高光譜數(shù)據(jù)與CWQI的最優(yōu)反演模型，其驗證數(shù)據(jù)集RMSE和MRE分別為0.768和18.1%，具有良好的反演精度。將該模型應(yīng)用于無人機高光譜影像的CWQI反演，可以較合理地反映出河流水質(zhì)的空間差異。盡管無人機航程、水流變化復(fù)雜等因素可能對無人機高光譜水質(zhì)反演產(chǎn)生影響，但水質(zhì)高光譜反演模型和無人機結(jié)合，可為河流水質(zhì)綜合評價提供一種大范圍覆蓋、可視化表達的技術(shù)手段，具有良好的業(yè)務(wù)化運行潛力。