謝勇，李凱云，，李家國，鮑云飛，朱利

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

0 引言

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

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

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

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

1.1 研究區(qū)概況

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

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

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

表1 Sentinel-2A波段信息

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

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

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

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

2 水污染遙感識別方法

2.1 水域提取

王大釗等[15]使用Sentinel-2A影像對比分析了四種水體提取方法。結(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)測斷面1 km緩沖區(qū)(圖1中紅色部分)內(nèi)使用ArcMap軟件生成100 m×100 m的均勻樣本點。由于1 km緩沖區(qū)范圍內(nèi)含有非水域部分，因此去除非水域樣本點后，最終獲取位于水域的樣本點數(shù)量942個，包含一般水體斷面樣本點731個、污染斷面的樣本點211個。其中，樣本點的3/4用于構(gòu)建污染水體遙感識別算法，1/4樣本點用于算法精度驗證。一般水體和污染水體斷面平均反射率光譜曲線和其標準差分布如圖2、圖3所示。

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

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

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

2.3 水污染遙感識別算法

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

目前，國內(nèi)外學(xué)者提出了多種波段組合算法對湖泊、河流水質(zhì)進行遙感識別，如使用單波段閾值法(RNIR、R705 nm、Rred)對礦區(qū)水體、湖泊水體進行遙感識別[18-20]，使用波段比值法(Rred/RNIR、R0.665 μm/R0.705 μm、R0.663 μm/R0.490 μm)對一些水質(zhì)參數(shù)(葉綠素a，懸浮物等)進行污染水體監(jiān)測[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]等識別算法對污染水體進行識別。本文結(jié)合污染水體和一般水體的光譜特征差異(圖2和圖3)，分別使用單波段閾值法ρgreen、波段比值法ρred/ρre1、波段差值法ρgreen-ρblue，對珠海市污染水體進行識別，并參考光譜反射率均值和標準差確定污染水體分布閾值。識別公式如表3所示。

表3 四種遙感識別算法公式

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

(2)

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

3 結(jié)果與討論

3.1 污染水體空間分布

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

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

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

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

3.2 遙感識別精度驗證及誤差分析

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

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖14 監(jiān)測斷面識別結(jié)果統(tǒng)計

3.5 污染水體時空分布及成因分析

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

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

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

4 結(jié)束語

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

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

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

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

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