謝勇，李凱云，，李家國(guó)，鮑云飛，朱利

(1.南京信息工程大學(xué) 地理科學(xué)學(xué)院，南京 210044；2.中國(guó)科學(xué)院空天信息創(chuàng)新研究院，北京 100101；3.北京空間機(jī)電研究所，北京 100094；4.生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心，北京 100094)

0 引言

河流與水庫(kù)是最重要的淡水資源之一，在調(diào)整支流、農(nóng)業(yè)灌溉、防洪排水以及調(diào)節(jié)全球生態(tài)平衡中起著重要作用，是全球經(jīng)濟(jì)可持續(xù)發(fā)展與人類生存和發(fā)展的重要保證。未經(jīng)處理而直接排放的工業(yè)廢水和生活污水，會(huì)導(dǎo)致水生態(tài)系統(tǒng)被破壞、河流湖泊污染加劇、部分水域富營(yíng)養(yǎng)化，給國(guó)家?guī)?lái)了嚴(yán)重的經(jīng)濟(jì)損失[1]。因此，污染水體的監(jiān)測(cè)在陸地水資源保護(hù)與污染治理中十分重要。

傳統(tǒng)的污染水體監(jiān)測(cè)主要采用定點(diǎn)定剖面采樣進(jìn)行水質(zhì)的化學(xué)因子分析，面對(duì)大尺度的水域時(shí)，傳統(tǒng)的監(jiān)測(cè)方法則表現(xiàn)出工作量大、投入大、周期長(zhǎng)的缺點(diǎn)。而遙感影像具有時(shí)效性與宏觀性強(qiáng)的特點(diǎn)，能清楚地反映出區(qū)域甚至整個(gè)流域污染現(xiàn)狀和空間分布特征[2]。目前，污染水體的遙感識(shí)別算法常采用經(jīng)驗(yàn)/半經(jīng)驗(yàn)?zāi)Ｐ?。其中，?jīng)驗(yàn)?zāi)Ｐ偷慕⑹菍?shí)測(cè)數(shù)據(jù)和遙感數(shù)據(jù)進(jìn)行統(tǒng)計(jì)相關(guān)分析，得出二者之間的定量關(guān)系。例如文獻(xiàn)[3-6]通過(guò)波段比值、波段差值等波段組合方式分別對(duì)懸浮物濃度、透明度、葉綠素濃度、濁度等水質(zhì)參數(shù)進(jìn)行反演和預(yù)測(cè)，結(jié)果表明實(shí)測(cè)數(shù)據(jù)和估算的水質(zhì)參數(shù)值具有良好的一致性，可有效監(jiān)測(cè)水質(zhì)參數(shù)的變化，為水質(zhì)快速監(jiān)測(cè)提供了有效的技術(shù)支持。半經(jīng)驗(yàn)方法則是以輻射傳輸模型為基礎(chǔ)，通過(guò)離水輻射率的定量關(guān)系反推水體各個(gè)成分的吸收系數(shù)和散射系數(shù)，同時(shí)依靠實(shí)測(cè)數(shù)據(jù)和水體固有光學(xué)參數(shù)的轉(zhuǎn)換關(guān)系，定量計(jì)算水質(zhì)參數(shù)，如順日布等[7]檢驗(yàn)了Nechad算法[8]、多波段準(zhǔn)分析算法(quasi-analytical algorithm，QAA)[9]、最優(yōu)化模型[10]等在黃河口水域懸浮物反演的適用性，結(jié)果表明，Nechad算法、QAA算法反演效果穩(wěn)定可靠，且具有相似的空間分布特征。除常用的經(jīng)驗(yàn)/半經(jīng)驗(yàn)算法模型外，深度學(xué)習(xí)方法也逐漸成為水質(zhì)監(jiān)測(cè)的技術(shù)手段，如Kupssinskü等[11]使用線性回歸模型、支持向量機(jī)模型、隨機(jī)森林模型、人工神經(jīng)網(wǎng)絡(luò)模型等方法對(duì)葉綠素和懸浮物濃度進(jìn)行反演，結(jié)果表明，隨機(jī)森林模型適用于葉綠素a和懸浮物的估算。根據(jù)以上研究現(xiàn)狀可知：以往多數(shù)學(xué)者研究地區(qū)為河口、湖泊等限定地區(qū)(如黃河口、珠江口等)，針對(duì)內(nèi)陸河流的污染水體反演研究鮮有報(bào)道；相對(duì)于以往所涉及的中低分辨率遙感數(shù)據(jù)，如MODIS、MERIS、Landsat等，Sentinel-2A具有獨(dú)特優(yōu)勢(shì)。Sentinel-2A的優(yōu)勢(shì)具體表現(xiàn)在：空間分辨率方面，Sentinel-2A數(shù)據(jù)優(yōu)于MODIS、MERIS數(shù)據(jù)；時(shí)間分辨率方面，Sentinel-2A數(shù)據(jù)優(yōu)于Landsat數(shù)據(jù)。另外，MODIS和MERIS傳感器已接近和超過(guò)衛(wèi)星設(shè)計(jì)壽命。

基于此，本研究利用具有較高時(shí)間和空間分辨率的Sentinel-2A影像，通過(guò)分析污染水體和一般水體在遙感影像上光譜特征差異，構(gòu)建污染水體識(shí)別算法，評(píng)價(jià)水體污染監(jiān)測(cè)成效，進(jìn)一步根據(jù)水體污染空間分布監(jiān)測(cè)結(jié)果探索珠海市內(nèi)陸河流污染來(lái)源，以期為珠海市及周邊地區(qū)水環(huán)境質(zhì)量安全控制提供依據(jù)。

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

1.1 研究區(qū)概況

珠海市位于廣東省中南部，珠江出?？谖靼?，是珠江三角洲中心城市之一，地處21°48′N(xiāo)～22°27′N(xiāo)，113°03′E～114°19′E之間。珠海市河流主要為西江的出海水道，有磨刀門(mén)水道、雞啼門(mén)水道、虎跳門(mén)水道和前山水道。珠海市屬典型的南亞熱帶季風(fēng)海洋性氣候，終年氣溫較高，氣候濕潤(rùn)，雨量充沛。近年來(lái)，隨著珠海市工業(yè)化、城市化進(jìn)程的加快，來(lái)自工業(yè)發(fā)展超標(biāo)排放的工業(yè)廢水和大量未經(jīng)處理的生活污水匯入河流，導(dǎo)致河流、湖泊等水污染現(xiàn)狀進(jìn)一步加重。水資源質(zhì)量不斷下降，水環(huán)境持續(xù)惡化，由水污染引起的缺水問(wèn)題和相關(guān)事故不斷出現(xiàn)，造成不良的社會(huì)影響和較大的經(jīng)濟(jì)損失，嚴(yán)重威脅著居民的生存和社會(huì)的可持續(xù)發(fā)展。

1.2 數(shù)據(jù)獲取

1)衛(wèi)星影像。Sentinel-2A多光譜成像衛(wèi)星于2015年6月23日發(fā)射，攜帶一枚多光譜成像儀MSI(譜段參數(shù)見(jiàn)表1)，可用于陸地監(jiān)測(cè)，提供植被、土壤和水覆蓋、海岸區(qū)域等圖像[12]，還可用于緊急救援服務(wù)。由于Sentinel-2A影像易獲取，且影像寬幅較大、覆蓋面積廣、重訪周期短(10 d)、在可見(jiàn)光波段及近紅外波段空間分辨率較高，因此適用于城市大、中型河流和湖泊污染水體監(jiān)測(cè)。

表1 Sentinel-2A波段信息

由于Sentnel-2A Level-1C級(jí)影像數(shù)據(jù)為經(jīng)正射校正和亞像元級(jí)幾何精校正后的大氣表觀反射率產(chǎn)品[13]，因此本文對(duì)Sentnel-2A數(shù)據(jù)進(jìn)行了大氣校正、重采樣和裁剪等預(yù)處理操作。大氣校正的目的是為了消除大氣和光照等因素對(duì)地物反射的影響，獲取地物真實(shí)反射率數(shù)據(jù)。同時(shí)，為了保持所有通道空間分辨率大小一致，對(duì)影像進(jìn)行了10 m重采樣操作。本文對(duì)Sentinel-2A影像的大氣校正和重采樣操作均借助于歐空局發(fā)布的Sen2cor工具集和SNAP軟件完成。

經(jīng)圖像清晰度目視判別，并結(jié)合獲取時(shí)間、太陽(yáng)光照條件、地表云量占比、影像區(qū)域覆蓋范圍等因素，最終篩選出覆蓋珠海市全區(qū)的2019年1月25日、2019年9月22日和2020年1月30日三期影像，進(jìn)行污染水體遙感識(shí)別算法構(gòu)建與評(píng)估以及污染水體空間分布識(shí)別和特征分析。

2)河流斷面類型。通過(guò)珠海政府網(wǎng)站(http：//www.zhuhai.gov.cn/)發(fā)布的2018年、2019年和2020年主要江河水質(zhì)月報(bào)確定污染水體斷面。監(jiān)測(cè)斷面共七個(gè)，即前山碼頭、石角咀水閘、尖峰大橋、雞啼門(mén)大橋、布洲、珠海大橋、虎跳門(mén)水道河口。判斷污染水體斷面的依據(jù)為《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》(GB 3838—2002)中Ⅳ類及以下，一般水體斷面為Ⅲ類及以上[14]。由于珠海政府網(wǎng)站發(fā)布的月報(bào)為月統(tǒng)計(jì)結(jié)果，為了更加準(zhǔn)確判別污染水體和一般水體斷面類別，利用影像獲取前后各一個(gè)月的月報(bào)指標(biāo)對(duì)水體斷面進(jìn)行判別，結(jié)果如表2所示，其中污染水體斷面4個(gè)，一般水體斷面15個(gè)。

表2 珠海市河流斷面水體判別結(jié)果

2 水污染遙感識(shí)別方法

2.1 水域提取

王大釗等[15]使用Sentinel-2A影像對(duì)比分析了四種水體提取方法。結(jié)果表明，改正后的歸一化差分水體指數(shù)(modify normalized difference water index，MNDWI)[16]提取精度為96.74%，效果較好。因此，本文使用MNDWI提取珠海市水域信息，如式(1)所示。

(1)

式中：ρgreen為Sentinel-2A影像的第3波段反射率；ρSWIR為Sentinel-2A影像的第11波段反射率。

2.2 水體樣本提取

基于表2，本文在監(jiān)測(cè)斷面1 km緩沖區(qū)(圖1中紅色部分)內(nèi)使用ArcMap軟件生成100 m×100 m的均勻樣本點(diǎn)。由于1 km緩沖區(qū)范圍內(nèi)含有非水域部分，因此去除非水域樣本點(diǎn)后，最終獲取位于水域的樣本點(diǎn)數(shù)量942個(gè)，包含一般水體斷面樣本點(diǎn)731個(gè)、污染斷面的樣本點(diǎn)211個(gè)。其中，樣本點(diǎn)的3/4用于構(gòu)建污染水體遙感識(shí)別算法，1/4樣本點(diǎn)用于算法精度驗(yàn)證。一般水體和污染水體斷面平均反射率光譜曲線和其標(biāo)準(zhǔn)差分布如圖2、圖3所示。

圖1 監(jiān)測(cè)斷面分布(圖中①～⑦依次為：布洲、虎跳門(mén)水道河口、前山水道、石角咀水閘、珠海大橋、尖峰大橋和雞啼門(mén)大橋監(jiān)測(cè)斷面緩沖區(qū))

圖2 一般水體和污染水體斷面平均光譜曲線

圖3 一般水體(總)和污染水體(總)斷面平均光譜反射率標(biāo)準(zhǔn)差

2.3 水污染遙感識(shí)別算法

通過(guò)衛(wèi)星手段監(jiān)測(cè)污染水體，主要是利用水體自身的遙感光譜信息建立其與水質(zhì)參數(shù)的關(guān)聯(lián)模型[17]。污染水體遙感識(shí)別算法的構(gòu)建原則是通過(guò)污染水體和一般水體之間的光譜差異特征，進(jìn)行波段組合運(yùn)算來(lái)增大不同水體之間的差異?；诖?，參考當(dāng)前其他水體異常提取方法，構(gòu)建水體識(shí)別算法對(duì)珠海市全行政區(qū)進(jìn)行污染水體識(shí)別提取。

目前，國(guó)內(nèi)外學(xué)者提出了多種波段組合算法對(duì)湖泊、河流水質(zhì)進(jìn)行遙感識(shí)別，如使用單波段閾值法(RNIR、R705 nm、Rred)對(duì)礦區(qū)水體、湖泊水體進(jìn)行遙感識(shí)別[18-20]，使用波段比值法(Rred/RNIR、R0.665 μm/R0.705 μm、R0.663 μm/R0.490 μm)對(duì)一些水質(zhì)參數(shù)(葉綠素a，懸浮物等)進(jìn)行污染水體監(jiān)測(cè)[21-23]，使用波段差值法(Rgreen-Rblue、R0.705 μm-R0.754 μm)[24-25]、比值植被指數(shù)(ratio vegetation index，RVI)[26]、歸一化差分葉綠素指數(shù)(normalized difference chlorophyll index，NDCI)[27]、清潔水體指數(shù)(water cleanliness index，WCI)[28]等識(shí)別算法對(duì)污染水體進(jìn)行識(shí)別。本文結(jié)合污染水體和一般水體的光譜特征差異(圖2和圖3)，分別使用單波段閾值法ρgreen、波段比值法ρred/ρre1、波段差值法ρgreen-ρblue，對(duì)珠海市污染水體進(jìn)行識(shí)別，并參考光譜反射率均值和標(biāo)準(zhǔn)差確定污染水體分布閾值。識(shí)別公式如表3所示。

表3 四種遙感識(shí)別算法公式

由圖2和圖3污染水體和一般水體的平均光譜曲線可以得出，污染水體水表反射率在波長(zhǎng)0.490～0.560 μm之間的差值比一般水體在此波長(zhǎng)范圍內(nèi)明顯偏高，且在0.490 μm和0.560 μm處標(biāo)準(zhǔn)差較小，分別為0.014和0.119，同時(shí)污染水體水表反射率在0.705 μm與0.842 μm處的差值比一般水體要高。由此可知，Sentinel-2A對(duì)應(yīng)此光譜范圍內(nèi)的四個(gè)波段能夠很好地體現(xiàn)污染水體的光譜特征。為了進(jìn)一步增加污染水體與一般水體的光譜特征差異，使用這兩組差值的和來(lái)增大兩種不同水體類型的差異。因此，選擇綠光波段、藍(lán)光波段、紅邊1波段和近紅外波段四個(gè)波段，構(gòu)建河流污染指數(shù)法(river pollution index，RPI)來(lái)識(shí)別珠海市污染水體，如式(2)所示。

(2)

式中：ρgreen、ρblue、ρre1、ρNIR分別對(duì)應(yīng)Sentinel-2A影像的第3波段、第2波段、第5波段、第8波段反射率。

3 結(jié)果與討論

3.1 污染水體空間分布

根據(jù)構(gòu)建的四種遙感識(shí)別算法分別對(duì)2019年1月25日、2019年9月22日和2020年1月30日三期影像進(jìn)行污染水體識(shí)別，結(jié)果如圖4至圖6所示。從識(shí)別結(jié)果可以看出，珠海市主要污染水體分布在香洲區(qū)南部和金灣區(qū)東部，而位于西北部斗門(mén)區(qū)污染水體面積分布則較少。其中，波段比值法對(duì)于污染水體的識(shí)別面積明顯比單波段閾值法、波段差值法和河流污染指數(shù)法偏大，主要體現(xiàn)在將部分河流和湖泊等一般水體識(shí)別為污染水體。

圖4 2019年1月25日水污染遙感識(shí)別結(jié)果

圖5 2019年9月22日水污染遙感識(shí)別結(jié)果

圖6 2020年1月30日水污染遙感識(shí)別結(jié)果

3.2 遙感識(shí)別精度驗(yàn)證及誤差分析

為了驗(yàn)證單波段閾值法、波段比值法、波段差值法、河流污染指數(shù)法四種遙感識(shí)別算法對(duì)珠海市污染水體的識(shí)別效果，選用精度驗(yàn)證公式(式(3)))進(jìn)行驗(yàn)證，驗(yàn)證結(jié)果如圖7所示。

(3)

式中：η為驗(yàn)證水體類型的正確率；n為驗(yàn)證水體類型正確樣本數(shù)；N5為驗(yàn)證水體類型的總樣本數(shù)。

根據(jù)四種污染水體遙感識(shí)別精度驗(yàn)證結(jié)果(圖7)，構(gòu)建的河流污染指數(shù)的識(shí)別準(zhǔn)確率最高，其中污染水體的識(shí)別準(zhǔn)確率高達(dá)92.59%，單波段閾值識(shí)別準(zhǔn)確率次之，為86.79%，波段比值法對(duì)污染水體識(shí)別效果較差，識(shí)別準(zhǔn)確率僅為64.81%。從整體上來(lái)看，河流污染指數(shù)法對(duì)污染水體和一般水體的識(shí)別準(zhǔn)確率較高，適用于珠海市河流污染水體的遙感識(shí)別和提取。

圖7 四種水污染遙感識(shí)別算法精度驗(yàn)證結(jié)果

3.3 提取結(jié)果空間精細(xì)分析

本文選取楊寮水庫(kù)和前山水道部分水域遙感識(shí)別結(jié)果(圖8至圖13)進(jìn)行誤差精細(xì)分析。其中，楊寮水庫(kù)為集中式生活飲用水水源，長(zhǎng)期處于達(dá)標(biāo)狀態(tài)，判定為一般水體；前山水道水質(zhì)情況根據(jù)影像月報(bào)統(tǒng)計(jì)結(jié)果，判定2019年1月25日和2019年9月22日處于污染狀態(tài)，2020年1月30日部分水體為污染狀態(tài)。

圖9 2019年1月25日前山水道水污染識(shí)別結(jié)果

圖10 2019年9月22日楊寮水庫(kù)水污染識(shí)別結(jié)果

圖11 2019年9月22日前山水道水污染識(shí)別結(jié)果

圖12 2020年1月30日楊寮水庫(kù)水污染識(shí)別結(jié)果

圖13 2020年1月30日前山水道水污染識(shí)別結(jié)果

1)單波段閾值法。從遙感識(shí)別結(jié)果可以看出，單波段閾值法對(duì)于污染水體和一般水體的區(qū)分較為明顯，由圖8(b)和圖11(b)識(shí)別結(jié)果可知，該識(shí)別算法將前山水道部分污染區(qū)域誤判為一般水體，造成污染水體識(shí)別面積偏低。

2)波段比值法。從圖8(c)和圖12(c)楊寮水庫(kù)識(shí)別結(jié)果可以看出，波段比值法將水庫(kù)、湖泊等一般水體誤判為污染水體，將圖9(c)前山水道部分污染水體誤判為一般水體。結(jié)合圖4至圖6波段比值法對(duì)三期影像的識(shí)別結(jié)果，該方法對(duì)于污染水體的識(shí)別面積偏高，造成了識(shí)別準(zhǔn)確率下降。

3)波段差值法。由波段差值法對(duì)于楊寮水庫(kù)的識(shí)別結(jié)果可以看出，該方法對(duì)于一般水體的識(shí)別效果較好；而從前山水道圖11(d)的識(shí)別結(jié)果可知，該污染區(qū)域被誤判為一般水體，造成了污染水體識(shí)別面積偏小，識(shí)別準(zhǔn)確偏低。

4)河流污染指數(shù)法。相波段比值法對(duì)于湖泊、水庫(kù)等一般水體的誤判，單波段閾值法和波段差值法對(duì)于污染水體的識(shí)別不穩(wěn)定。河流污染指數(shù)法對(duì)于一般水體識(shí)別效果較好，對(duì)前山水道長(zhǎng)期為污染水體的識(shí)別符合實(shí)際情況?？傮w而言，河流污染指數(shù)法對(duì)于兩種水體類型的識(shí)別精度較高，適合珠海市大、中型河流和水庫(kù)污染水體的識(shí)別。

3.4 斷面識(shí)別結(jié)果分析

根據(jù)精度驗(yàn)證和分析的結(jié)果，使用河流污染指數(shù)識(shí)別算法對(duì)監(jiān)測(cè)斷面進(jìn)行識(shí)別結(jié)果統(tǒng)計(jì)(圖14)。從圖14可以得出，2019年1月至2020年1月期間，前山碼頭、石角咀水閘和虎跳門(mén)水道河口污染程度均呈下降趨勢(shì)。布洲、雞啼門(mén)大橋、尖峰大橋和珠海大橋四個(gè)斷面河流污染指數(shù)雖在2019年9月份呈現(xiàn)上升趨勢(shì)，但仍均處于一般水體狀態(tài)。

圖14 監(jiān)測(cè)斷面識(shí)別結(jié)果統(tǒng)計(jì)

3.5 污染水體時(shí)空分布及成因分析

使用河流污染指數(shù)識(shí)別算法對(duì)2019年1月25日、2019年9月22日和2020年1月30日影像的污染水體和一般水體進(jìn)行面積統(tǒng)計(jì)(圖15)。從污染水體的空間分布上來(lái)看，主要分布在香洲區(qū)和金灣區(qū)，而斗門(mén)區(qū)污染水體面積相對(duì)于香洲區(qū)和金灣區(qū)較小。其中，香洲區(qū)污染區(qū)域主要分布在前山水道、馬騮洲水道。由于前山水道水域周?chē)ㄖ锩芏容^大，居民區(qū)和商業(yè)區(qū)產(chǎn)生的生活、生產(chǎn)廢水未經(jīng)專業(yè)處理，直接排入河流，導(dǎo)致前山水道水體污染。馬騮洲水道污染分布區(qū)域相對(duì)較集中，其中，裕安圍附近九洲港貨運(yùn)碼頭、洪灣港和西域碼頭地理位置相對(duì)集中分布，使得附近水域的懸浮泥沙含量較大，船舶排出的有機(jī)物直接對(duì)水體造成了污染。金灣區(qū)的主要污染區(qū)域分布在洪鶴大橋和珠三角環(huán)線高速(金灣段)附近，主要是由于工程施工對(duì)附近的水域影響較大，造成了水體大面積污染的情況。

從2019年至2020年污染面積變化結(jié)果可以看出(圖15)，珠海市污染水體面積占比逐漸減小，其中，2019年1月和2019年9月污染水體面積分布較大，污染水體占比分別為25.11%和21.97%，主要原因?yàn)楹辁Q大橋、珠三角環(huán)線高速、金海大橋工程建設(shè)造成的大面積水體污染。2020年1月，污染水體分布面積逐漸減小，僅為2.59%，一方面由于部分工程的推進(jìn)減少了污染水體的產(chǎn)生以及居民環(huán)保意識(shí)的增強(qiáng)，另一方面也證實(shí)了珠海市水環(huán)境治理工作卓有成效。

圖15 不同時(shí)期各行政區(qū)水體污染水體和一般水體面積占比統(tǒng)計(jì)

4 結(jié)束語(yǔ)

通過(guò)分析污染水體和一般水體的光譜特征差異特征以及構(gòu)建污染水體遙感識(shí)別算法，對(duì)珠海市內(nèi)陸河流、水庫(kù)進(jìn)行遙感識(shí)別，得出以下結(jié)論。

1)在構(gòu)建的四種污染水體識(shí)別算法中，河流污染指數(shù)識(shí)別精度最高，污染水體的識(shí)別準(zhǔn)確率為92.59%，適用于珠海市河流、水庫(kù)的污染水體識(shí)別。

2)2019年至2020年期間，珠海市污染水體面積在逐漸減少，主要污染區(qū)域在前山河道和珠三角環(huán)線高速(金灣段)。

3)由于珠海市政府網(wǎng)站公布的月報(bào)為月統(tǒng)計(jì)結(jié)果，與衛(wèi)星影像瞬時(shí)過(guò)境時(shí)間存在一定的誤差，使得遙感識(shí)別算法閾值不穩(wěn)定以及識(shí)別精度下降。本文通過(guò)選用影像前后各一個(gè)月的月報(bào)結(jié)果對(duì)斷面進(jìn)行判別來(lái)減小誤差。對(duì)于構(gòu)建的河流污染指數(shù)是否適用于其他地區(qū)，還需進(jìn)行深入探究和驗(yàn)證。

致謝：感謝十三五民用航天技術(shù)預(yù)先研究項(xiàng)目(國(guó)產(chǎn)衛(wèi)星信息智能分發(fā)技術(shù))為本論文提供的技術(shù)支持。