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

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

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

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

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

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

1 研究方法

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

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

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

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

(1)

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

1.1.2無(wú)人機(jī)載高光譜監(jiān)測(cè)

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

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

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

(2)

(3)

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

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

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

圖1 4項(xiàng)和21項(xiàng)水質(zhì)監(jiān)測(cè)指標(biāo)CWQI值的關(guān)系

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

1.3 PLSR反演模型建立與評(píng)價(jià)

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

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

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

2 結(jié)果與分析

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

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

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

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

為分析光譜分辨率對(duì)水體PLSR反演模型估算效果的影響，需找出水體CWQI高光譜遙感監(jiān)測(cè)最佳光譜分辨率。采用分段三次Hermite插值多項(xiàng)式法對(duì)原始數(shù)據(jù)進(jìn)行重采樣，采樣間隔依次為約 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)的建模與驗(yàn)證精度評(píng)價(jià)結(jié)果研究模型估算精度隨光譜分辨率的變化。結(jié)果(表2)表明，在各光譜分辨率的數(shù)據(jù)集中，原始光譜數(shù)據(jù)集的PLSR反演模型的估算效果最差；對(duì)原始數(shù)據(jù)進(jìn)行不同光譜分辨率的重采樣后，均有助于提升PLSR反演模型的估算效果。PLSR反演模型的驗(yàn)證集RMSE隨著光譜分辨率的降低呈現(xiàn)先減小后增大的趨勢(shì)。光譜分辨率為50 nm、波段數(shù)為11個(gè)時(shí)，驗(yàn)證集RMSE最小達(dá)0.790，MRE為18.63%，模型估算效果最優(yōu)。但總體上重采樣后各驗(yàn)證集的PLSR反演模型估算效果均處于接近的水平。該現(xiàn)象表明，在CWQI的遙感監(jiān)測(cè)中，過(guò)高的光譜分辨率可能會(huì)造成數(shù)據(jù)冗余和引入環(huán)境噪聲，造成模型復(fù)雜化并影響模型穩(wěn)定性。反之，在較低的光譜分辨率下，數(shù)量較少的波段可在基本涵蓋水體主要污染物特征譜段的同時(shí)，增強(qiáng)光譜信噪比，從而提高估算精度。

表2 不同波段數(shù)的PLSR建模效果評(píng)價(jià)

2.3 提取的主成分?jǐn)?shù)對(duì)PLSR反演模型的影響

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

圖3 不同主成分?jǐn)?shù)與驗(yàn)證集RMSE的關(guān)系

因此，確定建模最佳主成分?jǐ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為不同波長(zhǎng)處的水體光譜反射率。

2.4 PLSR反演模型估算效果評(píng)價(jià)

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

圖4 驗(yàn)證集樣本CWQI實(shí)測(cè)值與反演模型計(jì)算值比較

2.5 反演結(jié)果

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

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

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

3 結(jié) 語(yǔ)

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