連云港主要入海河流葉綠素a及綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演

王 林,王 祥,周 超,王新新,孟慶輝,陳艷攏

國家海洋環(huán)境監(jiān)測(cè)中心,遼寧 大連 116023

引 言

我國近岸和近海海域的主要污染物80%以上來自陸源排污,每年上百億噸的工業(yè)和生活廢水將大量的污染物質(zhì)攜帶排放入海,造成近岸海域水質(zhì)普遍較差;而入海河流在陸源污染物向海洋輸送的過程中發(fā)揮著重要作用。近幾年,隨著習(xí)近平生態(tài)文明思想的深入貫徹,“河長(zhǎng)制”、“灣長(zhǎng)制”等管理制度的廣泛推行,以及污染防治攻堅(jiān)戰(zhàn)的全面打響,入海河流與近岸海域的水質(zhì)狀況得到穩(wěn)步改善,但尚存在波動(dòng)情況,污染防治形勢(shì)依然嚴(yán)峻。連云港地處江蘇東北部,海州灣西岸,東瀕黃海,入海河流主要包括灌河、古泊善后河及臨洪河等,2020年以來水質(zhì)明顯提升,15條入海河流平均水質(zhì)全部消除劣Ⅴ類,但氮磷等營養(yǎng)鹽仍為主要污染物質(zhì),部分入海河流葉綠素a(chlorophyll a,Chl a)濃度非常高,呈現(xiàn)明顯的富營養(yǎng)化狀態(tài)。

隨著各類衛(wèi)星影像的不斷涌現(xiàn),遙感技術(shù)在海洋、湖泊、水庫及河流等水生態(tài)環(huán)境監(jiān)測(cè)與評(píng)價(jià)中得到廣泛應(yīng)用,尤其針對(duì)海洋和大型湖泊生態(tài)環(huán)境要素的定量遙感研究已持續(xù)多年[1-6],技術(shù)相對(duì)較為成熟。近幾年,高分衛(wèi)星與無人機(jī)遙感技術(shù)得到快速發(fā)展,針對(duì)狹長(zhǎng)河流水體生態(tài)環(huán)境要素的定量遙感研究逐漸出現(xiàn)[7],反演參數(shù)包括濁度(turbidity,Tur)[8]、懸浮顆粒物(suspended particulate matter,SPM)[9]、總氮(total nitrogen,TN)[10-11]、總磷(total phosphorus,TP)[10-11]、葉綠素a[12-15]及水體營養(yǎng)狀態(tài)[16]等。富營養(yǎng)化是入海河流和近岸海域典型的水質(zhì)污染問題,是水體中營養(yǎng)物質(zhì)富集導(dǎo)致浮游藻類過度生長(zhǎng)的生態(tài)異?，F(xiàn)象,對(duì)區(qū)域生態(tài)環(huán)境及人類生活安全構(gòu)成嚴(yán)重威脅[17]。因此,對(duì)入海河流中葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)(trophic level index,TLI(Σ))的定量遙感研究具有重要的實(shí)際應(yīng)用價(jià)值。Shi等[14]采用Gaofen-6和Sentinel-2衛(wèi)星數(shù)據(jù),通過機(jī)器學(xué)習(xí)模型反演了北京溫榆河、北運(yùn)河等水體中葉綠素a濃度,較好反映了內(nèi)陸河流葉綠素a濃度的空間分布特征。馬馳[15]以東北松嫩平原河流、湖泊水體為研究對(duì)象,建立了基于Sentinel-2A窄近紅外與深藍(lán)波段的葉綠素a濃度一元二次模型,反演結(jié)果顯示研究區(qū)內(nèi)葉綠素a濃度較低,水質(zhì)較好。張勇等[16]基于Landsat-8衛(wèi)星影像,建立了合肥市環(huán)城河總氮、總磷及氨氮濃度遙感反演模型,并運(yùn)用評(píng)分法進(jìn)行了富營養(yǎng)化評(píng)價(jià),最終確定了各河段水體的富營養(yǎng)化狀態(tài)。入海河流下游受海洋潮汐作用影響,存在一定的咸淡水交界區(qū)域,與一般的內(nèi)陸河流水體有一定差別,目前針對(duì)入海河流水體葉綠素a濃度及綜合營養(yǎng)狀態(tài)指數(shù)的定量遙感研究還相對(duì)較少。

研究以連云港市境內(nèi)灌河、古泊善后河、薔薇河及臨洪河等入海河流為研究對(duì)象,采用Sentinel-2 MSI影像及葉綠素a、總氮、總磷等現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù),建立了該區(qū)域主要入海河流葉綠素a及綜合營養(yǎng)狀態(tài)指數(shù)的遙感定量反演模型。

1 實(shí)驗(yàn)部分

1.1 研究區(qū)域與調(diào)查站點(diǎn)分布

研究區(qū)域選擇連云港市四條主要的入海河流,包括灌河、古泊善后河、薔薇河及臨洪河。其中,灌河是蘇北地區(qū)唯一在干流上沒有建閘的黃金入海通道,受兩岸工業(yè)排污和船舶航運(yùn)等影響,水質(zhì)較差且水體常年渾濁;古泊善后河作為蘇北地區(qū)行洪排澇的骨干河道,具備防洪、排澇、供水、灌溉、航運(yùn)等綜合功能,受沿線污染匯入影響,河流水質(zhì)超標(biāo)事件時(shí)有發(fā)生[18];薔薇河沿線入河支流排口多,主要受農(nóng)村農(nóng)業(yè)面源污染、酸洗石英砂點(diǎn)源污染等因素影響,河水水質(zhì)很不穩(wěn)定[19];臨洪河是薔薇河的下游河段,經(jīng)臨洪閘分界,于臨洪河口入海,故將薔薇河與臨洪河作為一條入海河流進(jìn)行后續(xù)研究。

現(xiàn)場(chǎng)調(diào)查于2022年6月24、25、27日進(jìn)行,調(diào)查參數(shù)包括總氮、總磷及葉綠素a等水質(zhì)要素。共計(jì)23個(gè)站點(diǎn),其中薔薇河/臨洪河9個(gè),古泊善后河7個(gè),灌河7個(gè),研究區(qū)域范圍及具體站點(diǎn)布設(shè)情況如圖1所示。

圖1 研究區(qū)域外業(yè)調(diào)查站點(diǎn)Fig.1 Field survey stations in the study area

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

1.2.1 水質(zhì)參數(shù)

葉綠素a濃度測(cè)量采用熒光法;總氮濃度測(cè)量采用堿性過硫酸鉀消解紫外分光光度法;總磷濃度測(cè)量采用鉬酸銨分光光度法。

1.2.2 衛(wèi)星影像數(shù)據(jù)

衛(wèi)星影像數(shù)據(jù)從歐洲航天局哥白尼開放數(shù)據(jù)訪問中心網(wǎng)站(https://scihub.copernicus.eu/)下載。采用2022年6月25日Sentinel-2A MSI L2A影像數(shù)據(jù)與現(xiàn)場(chǎng)實(shí)測(cè)水質(zhì)數(shù)據(jù)進(jìn)行匹配,提取與實(shí)測(cè)水質(zhì)站點(diǎn)對(duì)應(yīng)的光譜反射率數(shù)據(jù);建模后應(yīng)用于2022年6月25日Sentinel-2A MSI L2A影像,得到連云港主要入海河流葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演結(jié)果,并進(jìn)行空間分布特征分析。

1.3 綜合營養(yǎng)狀態(tài)指數(shù)分析方法

采用中國環(huán)境監(jiān)測(cè)總站關(guān)于“湖泊(水庫)富營養(yǎng)化評(píng)價(jià)方法及分級(jí)技術(shù)規(guī)定”方法進(jìn)行綜合營養(yǎng)狀態(tài)指數(shù)計(jì)算,見式(1)和式(2)

(1)

(2)

式(1)和式(2)中:TLI(Σ)(無量綱)為綜合營養(yǎng)狀態(tài)指數(shù);Wj(無量綱)為第j種評(píng)價(jià)指標(biāo)營養(yǎng)狀態(tài)指數(shù)的相關(guān)權(quán)重;TLI(j)(無量綱)為第j種評(píng)價(jià)指標(biāo)的營養(yǎng)狀態(tài)指數(shù)。

綜合營養(yǎng)狀態(tài)指數(shù)的分級(jí)為:TLI(Σ)<30時(shí)為貧營養(yǎng);30≤TLI(Σ)≤50時(shí)為中營養(yǎng);5070時(shí)為重度富營養(yǎng)。

河流水體富營養(yǎng)化受多種環(huán)境因子影響,其中最為重要的兩個(gè)參數(shù)即是氮、磷,而葉綠素a則是富營養(yǎng)化最重要的表征指標(biāo)[20],對(duì)水體綜合營養(yǎng)狀態(tài)的評(píng)估均具有決定作用。因此,將TN、TP及Chl a作為綜合營養(yǎng)狀態(tài)指數(shù)的評(píng)價(jià)指標(biāo),分別計(jì)算得到各評(píng)價(jià)指標(biāo)的營養(yǎng)狀態(tài)指數(shù)TLI(j),確定其相關(guān)權(quán)重指數(shù)Wj分別為0.282 76、0.296 72、0.420 52,進(jìn)而計(jì)算得到綜合營養(yǎng)狀態(tài)指數(shù)TLI(Σ)。

1.4 模型評(píng)價(jià)

為了篩選葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)遙感反演的最佳波段,采用相關(guān)系數(shù)(r)進(jìn)行評(píng)價(jià),見式(3)

(3)

葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)反演模型的優(yōu)劣,采用決定系數(shù)(R2)、平均絕對(duì)百分比誤差(mean absolute percentage error,MAPE)和均方根誤差(root mean square error,RMSE)進(jìn)行評(píng)價(jià),見式(4)、式(5)和式(6)

(4)

(5)

(6)

2 結(jié)果與討論

2.1 入海河流營養(yǎng)狀態(tài)背景

現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)分析可知,薔薇河/臨洪河、古泊善后河、灌河水質(zhì)均較差,總氮超標(biāo)現(xiàn)象最為突出,大多站點(diǎn)未達(dá)到水質(zhì)Ⅲ類考核要求,劣五類水質(zhì)站點(diǎn)比例達(dá)65%,可能與現(xiàn)場(chǎng)調(diào)查時(shí)間為豐水期,受沿線污染大量匯入影響較大有關(guān)。進(jìn)一步對(duì)比發(fā)現(xiàn),Chl a、TP及TLI(Σ)均值以薔薇河/臨洪河最高,分別達(dá)42.83 μg·L-1、0.30 mg·L-1、67.75,古泊善后河次之,灌河最低;TN均值以灌河最高,達(dá)2.96 mg·L-1,薔薇河/臨洪河次之,古泊善后河最低。具體調(diào)查結(jié)果如表1所示。

表1 連云港主要入海河流的營養(yǎng)狀態(tài)參數(shù)統(tǒng)計(jì)Table 1 Statistics on trophic status parameters for Lianyungang major seagoing rivers

2.2 遙感反演敏感波段篩選

提取與水質(zhì)現(xiàn)場(chǎng)實(shí)測(cè)站點(diǎn)相匹配的水體光譜反射率R(λ),并分別與Chl a、TLI(Σ)進(jìn)行線性擬合,各波段的相關(guān)系數(shù)如圖2所示。整體看來,400～2 200 nm波段范圍內(nèi)R(λ)與Chl a、TLI(Σ)均存在負(fù)相關(guān)關(guān)系,即R越高水體越渾濁,Chl a、TLI(Σ)則越低,表明渾濁水體光線透射較差,不利于浮游藻類生存生長(zhǎng);且TLI(Σ)與R(λ)的相關(guān)性略高于Chl a,與TLI(Σ)計(jì)算時(shí)除Chl a之外還引入了TP、TN有關(guān);400～750 nm可見光波段R(λ)與Chl a、TLI(Σ)的相關(guān)性明顯高于其他波段,可見光波段的相關(guān)性先增加后減小,490、560、665 nm三個(gè)波段的相關(guān)性最佳,其中R(λ)與Chl a的相關(guān)系數(shù)分別為-0.697、-0.681、-0.693,R(λ)與TLI(Σ)的相關(guān)系數(shù)分別為-0.728、-0.744、-0.706。故將490、560、665 nm三個(gè)波段作為連云港主要入海河流葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演的敏感波段。

圖2 Sentinel-2 MSI主要波段反射率與葉綠素a濃度、綜合營養(yǎng)狀態(tài)指數(shù)的相關(guān)系數(shù)Fig.2 Correlation coefficients of Sentinel-2 MSI major band reflectance with Chl a concentration and TLI(Σ)

2.3 遙感反演模型建立

以490、560、665 nm為敏感波段,并選擇單波段、雙波段、三波段組合形式,在Chl a對(duì)數(shù)坐標(biāo)下采用線性、乘冪及指數(shù)模型進(jìn)行擬合分析,R2最高而MAPE、RMSE最低時(shí),則模型最佳,結(jié)果如表2所示?？砂l(fā)現(xiàn),單波段模型中,以R(665)為自變量,Chl a對(duì)數(shù)坐標(biāo)下的乘冪函數(shù)表現(xiàn)最佳,擬合R2為0.67,MAPE為47.34%,RMSE為12.89 μg·L-1;雙波段模型中,常用的波段比值模型表現(xiàn)較差,MAPE均超過60%,不適用于連云港主要入海河流水體,而二元一次函數(shù)模型略優(yōu)于波段比值模型,但整體表現(xiàn)均較差;采用三元一次函數(shù)擬合的三波段模型表現(xiàn)一般,也不適用于該區(qū)域入海河流水體。因此,將lg(Chl a)=0.267×R(665)-0.982作為連云港主要入海河流葉綠素a濃度遙感定量反演的最優(yōu)模型。

表2 連云港主要入海河流葉綠素a濃度遙感定量反演模型的精度評(píng)價(jià)Table 2 Accuracy evaluation of Chl a concentration in Lianyungang major seagoing rivers by using remote sensing inversion models

同樣以490、560、665 nm為敏感波段,選擇單波段、雙波段、三波段組合形式,并采用線性、乘冪函數(shù)進(jìn)行綜合營養(yǎng)狀態(tài)指數(shù)模型擬合分析,R2最高而MAPE、RMSE最低時(shí),則模型最佳,結(jié)果如表3所示。可發(fā)現(xiàn),單波段模型中,以R(560)為自變量,并采用乘冪函數(shù)進(jìn)行擬合的模型表現(xiàn)最佳,擬合R2為0.61,MAPE為4.36%,RMSE為3.45;雙波段模型中,波段比值模型表現(xiàn)較差,MAPE均超過6%,但二元一次函數(shù)模型明顯優(yōu)于波段比值模型;采用三元一次函數(shù)擬合的三波段模型表現(xiàn)較佳,但遜于單波段模型。因此,將TLI(Σ)=44.898×R(560)-0.241作為連云港主要入海河流綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演的最優(yōu)模型。

表3 連云港主要入海河流綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演模型的精度評(píng)價(jià)Table 3 Accuracy evaluation of TLI(Σ) in Lianyungang major seagoing rivers by using remote sensing inversion models

2.4 空間分布特征

將上節(jié)中最優(yōu)模型應(yīng)用于2022年6月25日Sentinel-2A MSI L2A影像,反演得到連云港主要入海河流Chl a和TLI(Σ),如圖3和圖4所示。整體看來,薔薇河/臨洪河Chl a最高,古泊善后河次之,灌河最低,且各河流上游Chl a普遍高于下游,可能與漲潮時(shí)海水涌入河流下游有關(guān);類似地,薔薇河/臨洪河TLI(Σ)最高,古泊善后河次之,灌河最低,且各河流上游富營養(yǎng)化程度普遍高于下游。進(jìn)一步分析可知,薔薇河/臨洪河80%以上水域Chl a大于10 μg·L-1,且大于20 μg·L-1水域占比48.68%;而灌河90%以上水域Chl a則小于10 μg·L-1,小于5 μg·L-1水域占比53.89%;古泊善后河85%以上水域Chl a在5～20 μg·L-1范圍,5～10 μg·L-1水域占比52.52%。薔薇河/臨洪河、古泊善后河及灌河基本均為富營養(yǎng)水體,其中薔薇河/臨洪河以中度富營養(yǎng)水體為主占比91.15%,重度富營養(yǎng)水體占比8.62%;古泊善后河也以中度富營養(yǎng)水體為主占比85.70%,而輕度富營養(yǎng)水體占比11.85%;灌河則以輕度富營養(yǎng)水體為主占比72.52%,中度富營養(yǎng)水體占比25.04%。具體統(tǒng)計(jì)結(jié)果如表4所示。

表4 不同區(qū)間葉綠素a和綜合營養(yǎng)狀態(tài)指數(shù)所占水域面積的比例統(tǒng)計(jì)Table 4 Interval statistics of watershed areas with different levels of Chl a and TLI(Σ)

圖3 連云港主要入海河流葉綠素a濃度(2022.06.25)Fig.3 Chl a concentration in major seagoing rivers in Lianyungang (2022.06.25)

圖4 連云港主要入海河流綜合營養(yǎng)狀態(tài)指數(shù)(2022.06.25)Fig.4 TLI(Σ) in major seagoing rivers in Lianyungang (2022.06.25)

以上模型研究發(fā)現(xiàn),分別以665和560 nm為敏感波段的單波段模型更適用于連云港主要入海河流葉綠素a濃度和綜合營養(yǎng)狀態(tài)指數(shù)的遙感定量反演,這與湖泊、海洋水體存在一定差別。在湖泊、海洋相關(guān)研究中,560和665 nm波段多用于懸浮物濃度的遙感反演,且560 nm波段還作為參比波段用于葉綠素a濃度的遙感反演。分析原因在于,與湖泊或海洋相比,河流水體的渾濁度更高,水體光學(xué)特性更為復(fù)雜,不同河流或河段光線透射進(jìn)入水體的深度存在一定差異,導(dǎo)致浮游藻類的生存生長(zhǎng)狀況呈現(xiàn)較大的不同,進(jìn)而使得葉綠素a濃度、綜合營養(yǎng)狀態(tài)指數(shù)與水體渾濁度具有負(fù)相關(guān)關(guān)系,且常用于渾濁度遙感反演的特征波段成為葉綠素a濃度、綜合營養(yǎng)狀態(tài)指數(shù)反演的最優(yōu)波段。

將反演模型應(yīng)用于衛(wèi)星影像,得到灌河、古泊善后河、薔薇河/臨洪河Chl a和TLI(Σ)空間分布結(jié)果,發(fā)現(xiàn)灌河Chl a、TLI(Σ)均最低,古泊善后河次之,薔薇河/臨洪河最高,這與河流水體的渾濁程度有直接關(guān)聯(lián)。進(jìn)一步對(duì)比分析可知,灌河各調(diào)查站點(diǎn)SPM均值為443.81 mg·L-1,古泊善后河為65.41 mg·L-1,薔薇河/臨洪河為33.34 mg·L-1,Chl a、TLI(Σ)與SPM呈負(fù)相關(guān)關(guān)系,表明SPM越高,水體越渾濁,光線透射越差,不利于浮游藻類生長(zhǎng),Chl a越低;TLI(Σ)由Chl a、TN、TP計(jì)算而得,其中Chl a是TLI(Σ)最重要的表征參數(shù),權(quán)重指數(shù)最高,對(duì)TLI(Σ)影響最大,就灌河而言,即使TN、TP相對(duì)較高,但其Chl a尤其低,綜合作用下TLI(Σ)仍最低。

3 結(jié) 論

(1)受河流沿線污染匯入影響,薔薇河/臨洪河、古泊善后河、灌河水質(zhì)均較差。其中,薔薇河/臨洪河Chl a、TP及TLI(Σ)均值最高,古泊善后河次之,灌河最低;灌河TN均值最高,薔薇河/臨洪河次之,古泊善后河最低。

(2)葉綠素a濃度、綜合營養(yǎng)狀態(tài)指數(shù)與可見光波段反射率的相關(guān)性明顯高于其他波段,最佳波段為490、560和665 nm。其中,R(λ)與Chl a的相關(guān)系數(shù)分別為-0.697、-0.681、-0.693,R(λ)與TLI(Σ)的相關(guān)系數(shù)分別為-0.728、-0.744、-0.706。

(3)以R(665)為自變量,在葉綠素a濃度對(duì)數(shù)坐標(biāo)下的乘冪模型為其遙感定量反演的最優(yōu)模型:lg(Chl a)=0.267×R(665)-0.982(R2=0.67,MAPE=47.34%,RMSE=12.89 μg·L-1);而以R(560)為自變量的乘冪模型是綜合營養(yǎng)狀態(tài)指數(shù)遙感定量反演的最優(yōu)模型:TLI(Σ)=44.898×R(560)-0.241(R2=0.61,MAPE=4.36%,RMSE=3.45)。

(4)將最優(yōu)模型應(yīng)用于2022年6月25日Sentinel-2A MSI L2A影像,得到連云港主要入海河流葉綠素a濃度和綜合營養(yǎng)狀態(tài)指數(shù)的遙感反演結(jié)果,發(fā)現(xiàn)薔薇河/臨洪河中度富營養(yǎng)水體的占比為91.15%,重度富營養(yǎng)水體占比為8.62%;古泊善后河中度富營養(yǎng)水體占比為85.70%,輕度富營養(yǎng)水體占比11.85%;而灌河則以輕度富營養(yǎng)水體為主占比72.52%,中度富營養(yǎng)水體占比25.04%。