連云港主要入海河流葉綠素a及綜合營養狀態指數遙感定量反演

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

國家海洋環境監測中心,遼寧 大連 116023

引 言

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

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

研究以連云港市境內灌河、古泊善后河、薔薇河及臨洪河等入海河流為研究對象,采用Sentinel-2 MSI影像及葉綠素a、總氮、總磷等現場實測數據,建立了該區域主要入海河流葉綠素a及綜合營養狀態指數的遙感定量反演模型。

1 實驗部分

1.1 研究區域與調查站點分布

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

現場調查于2022年6月24、25、27日進行,調查參數包括總氮、總磷及葉綠素a等水質要素。共計23個站點,其中薔薇河/臨洪河9個,古泊善后河7個,灌河7個,研究區域范圍及具體站點布設情況如圖1所示。

圖1 研究區域外業調查站點Fig.1 Field survey stations in the study area

1.2 數據獲取

1.2.1 水質參數

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

1.2.2 衛星影像數據

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

1.3 綜合營養狀態指數分析方法

采用中國環境監測總站關于“湖泊(水庫)富營養化評價方法及分級技術規定”方法進行綜合營養狀態指數計算,見式(1)和式(2)

(1)

(2)

式(1)和式(2)中:TLI(Σ)(無量綱)為綜合營養狀態指數;Wj(無量綱)為第j種評價指標營養狀態指數的相關權重;TLI(j)(無量綱)為第j種評價指標的營養狀態指數。

綜合營養狀態指數的分級為:TLI(Σ)<30時為貧營養;30≤TLI(Σ)≤50時為中營養;5070時為重度富營養。

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

1.4 模型評價

為了篩選葉綠素a和綜合營養狀態指數遙感反演的最佳波段,采用相關系數(r)進行評價,見式(3)

(3)

葉綠素a和綜合營養狀態指數反演模型的優劣,采用決定系數(R2)、平均絕對百分比誤差(mean absolute percentage error,MAPE)和均方根誤差(root mean square error,RMSE)進行評價,見式(4)、式(5)和式(6)

(4)

(5)

(6)

2 結果與討論

2.1 入海河流營養狀態背景

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

表1 連云港主要入海河流的營養狀態參數統計Table 1 Statistics on trophic status parameters for Lianyungang major seagoing rivers

2.2 遙感反演敏感波段篩選

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

圖2 Sentinel-2 MSI主要波段反射率與葉綠素a濃度、綜合營養狀態指數的相關系數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對數坐標下采用線性、乘冪及指數模型進行擬合分析,R2最高而MAPE、RMSE最低時,則模型最佳,結果如表2所示?？砂l現,單波段模型中,以R(665)為自變量,Chl a對數坐標下的乘冪函數表現最佳,擬合R2為0.67,MAPE為47.34%,RMSE為12.89 μg·L-1;雙波段模型中,常用的波段比值模型表現較差,MAPE均超過60%,不適用于連云港主要入海河流水體,而二元一次函數模型略優于波段比值模型,但整體表現均較差;采用三元一次函數擬合的三波段模型表現一般,也不適用于該區域入海河流水體。因此,將lg(Chl a)=0.267×R(665)-0.982作為連云港主要入海河流葉綠素a濃度遙感定量反演的最優模型。

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

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

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

2.4 空間分布特征

將上節中最優模型應用于2022年6月25日Sentinel-2A MSI L2A影像,反演得到連云港主要入海河流Chl a和TLI(Σ),如圖3和圖4所示。整體看來,薔薇河/臨洪河Chl a最高,古泊善后河次之,灌河最低,且各河流上游Chl a普遍高于下游,可能與漲潮時海水涌入河流下游有關;類似地,薔薇河/臨洪河TLI(Σ)最高,古泊善后河次之,灌河最低,且各河流上游富營養化程度普遍高于下游。進一步分析可知,薔薇河/臨洪河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%。薔薇河/臨洪河、古泊善后河及灌河基本均為富營養水體,其中薔薇河/臨洪河以中度富營養水體為主占比91.15%,重度富營養水體占比8.62%;古泊善后河也以中度富營養水體為主占比85.70%,而輕度富營養水體占比11.85%;灌河則以輕度富營養水體為主占比72.52%,中度富營養水體占比25.04%。具體統計結果如表4所示。

表4 不同區間葉綠素a和綜合營養狀態指數所占水域面積的比例統計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 連云港主要入海河流綜合營養狀態指數(2022.06.25)Fig.4 TLI(Σ) in major seagoing rivers in Lianyungang (2022.06.25)

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

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

3 結 論

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

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

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

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