哨兵2 衛星綜合水質指標的河流水質遙感監測方法

王歆暉，田 華，季鐵梅，鞏彩蘭，胡 勇，李 瀾，何志杰

（1.中國科學院大學 物理科學學院，北京 100049；2.中國科學院 上海技術物理研究所，上海 200083；3.上海市水文總站，上海 200232；4.中國科學院 紅外成像光譜重點實驗室，上海 200083）

0 引言

內陸河流作為重要的資源和載體，是生態系統的綠色生命線，關系到人類的生存和發展。然而，河流水環境極易受到污染，特別是城市河道，隨著現代城市的飛速發展與擴張，人口密度增加，給城市河流水環境帶來了挑戰。從2017 年起，上海市就推行“河長制”，目的在于加強對上海市市管河道的監管，對水污染進行防治，保護水資源。因此，保障居民的正常生產生活及推進可持續發展，對城市河流水質快速、準確的監測工作提出了高要求。目前，河流水質監測主要采用人工實地采樣與實驗室分析的方法，該方法結果準確，但只能獲取河流截面采樣點數據［1-2］。與其相比，內陸河流遙感監測具有以下優勢：1）遙感技術具有快速、大范圍和周期性的特點，可以彌補常規定時、定點監測的不足，節省大量人力、物力和財力［3］；2）隨著對水體光譜特性的深入研究以及水質參數反演模型的改進，遙感圖像能夠較為準確地獲取河流水質信息，與地面實測水文參數、地理位置信息等結合，可有效發現污染水體以及水質變化趨勢，作為水質管理的參考。

隨著新的衛星傳感器不斷升空，時間分辨率和空間分辨率不斷提高，水質遙感監測從定性發展到定量，可實現遙感反演的水質指標也由葉綠素a、懸浮物，發展到總氮、總磷、高錳酸鹽指數等［4-7］。當前，利用衛星遙感技術對水質監測的研究主要集中在區域跨度較大的研究對象上，例如大型湖泊水庫和河口海岸的監測，對小尺度和微型態環境問題研究相對較少。MATTEWS［8］利用MERIS 數據發展了改進的最大峰值算法，對內陸和近岸水體葉綠素a 濃度進行探測，研究了非洲南部50 多個湖庫水體從2002 到2012 年的富營養化及其變化狀況。HOU等［9］利用MODIS 地表反射率數據645 nm 波段反演了長江中下游102 個大型湖泊的水體總懸浮濃度，并分析了其變化趨勢。國內學者臧友華［10］將遺傳算法和支持向量機組合，提出一種新的反演模型GA-SVM 模型，采用ETM+數據對渭河水質參數葉綠素a、透明度、總磷進行了反演。鞏彩蘭等［11］對黃浦江進行了水體光譜測量和同步水質采樣，分析了各水質參數與反射率及相關反射率特征的相關性，建立了常規水質參數與反射率之間的關系模型。以上研究主要是針對部分具有光學特性的遙感水質參數或者是國家《地表水環境質量標準》中規定的單一水質參數，雖然簡單直觀地反映了單因子污染狀況，但是缺乏對河流水質的綜合研究。因此，本文利用5 項主要水質參數建立河流綜合水質遙感模型，以綜合水質指標來表征河流水質狀況。

本文采用Sentinel-2A 衛星遙感影像數據和同步地面河流斷面實測水質參數數據，對5 項水質參數溶解氧、高錳酸鹽指數、五日生化需氧量、氨氮以及總磷進行因子分析，在此基礎上建立水質遙感反演的模型，得到綜合水質指標用于確定水質類別。

1 實驗數據及預處理

1.1 實驗數據

本文采用歐洲航空局發射的Sentinel-2A 衛星搭載的多光譜成像儀獲取的遙感影像，影像成像日期為2018 年4 月9 日，區域為上海市，具體研究選擇上海市青浦區和松江區的部分河流流域，研究區域如圖1所示。影像覆蓋13個光譜波段，幅寬達290 km，重訪周期10 d。從可見光和近紅外到短波紅外，具有不同的空間分辨率。

圖1 研究區域衛星影像band 8、4、3 假彩色合成圖Fig.1 False color composite images of bands 8，4，and 3 in the study area

考慮空間分辨率，本研究用到的波段為band 2、band 3、band 4、band 8（見表1）。

表1 Sentinel-2A 衛星部分波段信息Tab.1 Partial band information of Sentinel-2A satellite

水質參數采樣數據獲取的是衛星成像同一天的上海市河流斷面實地采樣分析數據，按照《地表水環境質量標準（GB 3838—2002）》對上海市主要河流94 個斷面采樣點溶解氧（DO）、高錳酸鹽指數（CODMn）、五日生化需氧量（BOD5）、氨氮（NH3-N）、總磷（TP）共5 項指標進行評價，如圖2 所示。總體上，溶解氧、高錳酸鹽指數、總磷評價等級較好，基本達到3 類水的標準。超標的主要水質參數是氨氮，通常來源于飼料、水生動物的排泄物、肥料及動植物尸體分解，未處理或處理過的城市生活和工業廢水、各種浸濾液和地表徑流是氨氮含量升高的主要原因［12］。

圖2 各水質參數所屬水體類型分布情況圖Fig.2 Water grade summary of each water quality parameter

1.2 實驗數據預處理

對于經過輻射定標等預處理后的反射率影像，運用歸一化差異水體指數（NDWI）對2018年4月9 日的實驗影像研究區進行河流水體提取。對NDWI 圖像設置掩模模板提取水體，閾值往往設為0，即NDWI 圖像大于0 的像素被認為是水體［13］：

式中：bandgreen為第3 波段（綠光）反射率；bandNIR為第8 波段（近紅外）反射率。

將影像與水體提取掩模模板疊加對比，掩模結果如圖3 所示，其中白色部分為識別的水體，而黑色部分為非水體。目視解譯結果表明，歸一化差異水體指數NDWI 能夠較好地識別出河流。

圖3 2018 年4 月研究區水體NDWI 提取掩模結果Fig.3 Mask results of water extraction with NDWI in the study area in April 2018

2 基于因子分析的水質遙感反演方法

2.1 數據波段篩選

在正式進行建模前對波段進行了篩選，結合單波段與各主要水質參數的相關性以及OIF 指數，來選擇對水質參數變化較為敏感的特征波段或者波段組合。本文分析了各單波段與各水質參數的相關性，整體相關性較低，Pearson 系數整體不超過0.4，所以Sentinel-2 單波段并不適合進行水質遙感反演，需要考慮其他波段組合。研究區主要超標的水質參數為氨氮，各波段與氨氮相關系數按照從大到小順序排列，波段順序為band 3、band 4、band 2、band 8。OIF 指數［14］是進行波段篩選的依據之一，通常波段數據的標準差越大，包含的信息量越多；而波段之間的相關系數越小，則波段之間獨立性越高，數學表達式為

式中：Si為第i個波段的標準差；Rij為第i、j兩波段的相關系數。

不同三波段組合的OIF 指數見表2，按照從大到小的順序排列來選擇最優組合方案，理想組合順序為：1）band 2、band 4、band 8；2）band 3、band 4、band 8；3）band 2、band 3、band 8。因組合1 和組合2的OIF 指數值差異并不明顯，結合研究區主要超標的水質參數氨氮與波段的相關性，最終確定選用組合2 波段3、4、8 組合。

2.2 綜合因子分析理論

因子分析是一種處理多元變量數據的統計方法，研究眾多變量之間的內在關系，對可觀測變量進行概括和抽象，用較少的因子變量來最大程度上解釋原始的觀測變量［15-16］，其主要表達形式為

式中：X1，X2，…，Xp為p個存在相關關系的可觀測變量；F1，F2，…，Fm為m個互相獨立的因子，又稱為公共因子，是不可觀測的變量；ε1，ε2，…，εp為特殊因子，是無法包含在公共因子中的部分；aij(1≤i≤p，1≤j≤m)為第i個變量與第j個因子的相關系數。

基于因子分析建立遙感反演模型，首先計算原始變量的相關系數矩陣，然后利用相關系數矩陣結合反演模型和先驗知識提取公共因子，通常只取方差大于1 或者特征值大于1 的因子，或者按照因子的累計方差貢獻率來確定，一般要求至少達到70%。因子Fj的方差貢獻是式（3）中對應列相關系數的平方和，根據如下公式計算方差貢獻率，方差貢獻率進行累加即為累計方差貢獻率：

接著進行因子旋轉，其目的是通過坐標變換使因子實際意義更易解釋，最后根據如下公式計算綜合因子得分：

式中：ωj(1≤j≤n)為對應因子Fj的方差貢獻率。

2.3 基于綜合水質因子的遙感反演模型

由于水質參數本身數量級的不一致，為了避免對模型造成影響，對數據進行最大最小化標準化，將其映射到-1 至1 范圍內。采用因子分析法對研究區水質參數進行分析，首先KMO 值為0.763，大于閾值0.5，說明了變量之間是存在相關性的，符合要求；然后是Bartlett 球形檢驗的結果，Sig 檢驗概率值為0.000，小于閾值0.05，說明結果顯著，表明數據可以進行因子分析。

河流水質參數因子分析的主要因子成分如圖4所示。主要因子1 的特征值為2.8，方差貢獻率為56.74%；主要因子2 的特征值為0.95，方差貢獻率為18.93%；前3 個主要因子的方差貢獻率之和為87.5%。因此，水質參數監測數據原始5 個變量的信息能夠用變換后的3 個主要因子來表達。

圖4 監測數據的主要因子成分Fig.4 Major factor components of the monitoring data

由表3 可知，河流水質參數因子分析的旋轉載荷矩陣和貢獻率，主要抽取了方差貢獻率較大的前3 個主要因子成分，因子成分F1主要受高錳酸鹽指數和五日生化需氧量控制，系數分別為0.928 和0.832，這兩項指標越大，反映水體受有機物的污染越嚴重。因子成分F2主要由總磷和氨氮控制，系數分別為0.921 和0.649。通常是城市生活洗滌污水、垃圾等包含的無機鹽類，與水體的富營養化有關。因子成分F3則主要由溶解氧控制，系數值最大為0.970。

表3 因子旋轉載荷矩陣與方差貢獻率Tab.3 Factor loading matrix and variance contribution rate

5 項主要水質參數作為式（1）中各可觀測變量，經過因子旋轉后借助表3 的因子旋轉載荷矩陣，可以得到如下3 個主要因子成分的表達式：

將各個公共因子的方差貢獻率等參數代入式（5），得到水質綜合因子得分模型為

為了得到最終水質遙感反演模型，需要建立各個水質參數和遙感影像波段3、4、8 之間的關系。以影像波段反射率數據作為自變量，水質參數數據作為因變量，采用最小二乘法確定多元線性回歸模型的待定參數，當實測值與回歸模型得到的擬合值的差值平方和達到最小時的模型參數即為所求待定參數。根據研究數據得到的回歸模型結果見表4。由回歸分析結果可見，溶解氧、高錳酸鹽與波段之間的關系相對較密切。

表4 不同水質參數的線性回歸模型Tab.4 Regression models of different water quality parameters

根據得出的綜合因子得分模型以及水質參數和波段的回歸分析結果，計算得到最終水質遙感反演模型為

3 水質遙感反演結果與精度分析

對預處理［17］后的2019 年4 月9 日的研究區影像，應用綜合水質因子反演模型進行反演，得到研究區的綜合水質圖，如圖5 所示。一般情況下，綜合水質F值越小，代表水體越潔凈，否則水體則存在水質參數超標，或者受到污染的可能性。從結果圖來看，主干河流的值基本低于支流，是符合實際水質情況的。

圖5 2019 年4 月研究區綜合水質圖Fig.5 Comprehensive water quality map of the study area in April 2019

在綜合水質圖的基礎上，據不同水質等級的情況分別設定閾值，最終確定不同情況下綜合水質F對應的水質等級，繪制綜合水質等級圖，如圖6 所示。因衛星遙感影像分辨率的限制，圖中小河流存在斷裂、破碎等狀況。從結果圖看Ⅲ類和Ⅳ類水主要分布在主干河流，相對而言，Ⅴ類水體在支流較多，整體水質狀況良好，劣5 類水體只出現在局部區域。

隨機選取30 個樣本數據，將反演的綜合水質等級結果與實際測量采樣的結果進行對比驗證，如圖7 所示，水質類別2、3、4、5、6 分別對應標準水質類別Ⅱ類、Ⅲ類、Ⅳ類、Ⅴ類以及劣Ⅴ類。

圖7 實測水質與反演綜合水質類別對比圖Fig.7 Comparison of the measured water quality and the inversion comprehensive water quality

在大河道區域準確度能夠達到71.2%，但是細小的支流準確度只能夠達到42.8%，誤差的主要來源如下：1）地面實測采樣水質數據與衛星影像時間和空間的不一致性；2）僅以5 項主要水質參數表征河流水質狀況，仍存在其他因素的干擾；3）衛星數據空間分辨率的限制會帶來的小河流混合像元問題。

4 結束語

利用2018 年4 月Sentinel-2A 遙感影像以及同步獲取的河流斷面實測水質參數數據，結合5 項主要水質參數溶解氧、高錳酸鹽指數、五日生化需氧量、氨氮以及總磷，在因子分析的基礎上建立水質遙感反演的模型，得到綜合水質指標用于確定水質類別。實驗結果表明：1）上海市主要超標的水質參數指標為氨氮，水污染治理應著重處理難降解的有機物以及含氮的可溶性無機物等；2）綜合水質指標值越小，內陸河流水質狀況越良好，水體更清潔；3）本文方法對于大河流的水質等級判定的正確率達到了71.2%。基于因子分析的水質遙感反演模型能夠應用于內陸河流水質遙感監測，為水環境部門提供參考信息。但是本文僅考慮了5 項水質參數，對于總體水質狀況的表達還略有不足，需要考慮其他水質遙感監測指標，輔助綜合水質類別判定；另外，該模型僅適用于一個季節，獲取長時間序列的影像以及同步水質參數數據，能夠對模型進行改進，并且利用多時相數據對內陸河流水質狀況進行動態變化研究。