大氣校正算法對高光譜反演水體葉綠素a 濃度的影響

孟祥亮，馮建飛，付萍杰，?，張家威，張雨煊，孟飛

(1.山東建筑大學測繪地理信息學院，山東 濟南 250101；2.山東省生態環境監測中心，山東 濟南 250101；3.山東科技大學測繪與空間信息學院，山東 青島 266000)

0 引言

水體葉綠素a(Chla)是水體富營養化程度的重要指標之一，對湖庫葉綠素a 濃度的廣域定量監測和實時定點分析是監測湖庫水體營養狀態的基本需求。水體葉綠素a 濃度的監測難點之一在于其空間異質性，地面監測雖然精度較高，但局限于某一特定時刻，難以滿足動態變化的監測。 遙感技術具有監測范圍廣、監測時序長的優勢[1]，但大氣中水汽分子、氣溶膠等物質會消耗掉大部分的離水輻射信號(90%～95%)，且內陸水體空間布局分散，地理狀況多變，其光學特性表現為復雜性和區域性[2－4]。 大氣校正能夠消除部分大氣影響，還原目標地物的光譜信號，高精度的大氣校正可以降低不同空間和不同時相影像大氣校正誤差對模型構建和應用帶來的不確定性，是提升水質遙感監測模型時空移植性的前提，對于提高水色反演模型的準確性和普適性具有重要意義。

學者們針對水質反演的大氣校正方法開展了多項研究。 杜挺等[5]采用HJ－1B 多光譜影像，對太湖流域進行的反演中發現光譜超立方體的快速視線大氣分析(Fast Line－of－Fight Atmospheric Analysis of Spectral Hypercubes， FLAASH)和對太陽光譜中衛星信號的二次分析(Second Simulation of the Satellite Signal in the Solar Spectrum， 6S)算法效果較好，但快速大氣校正( Quick Atmospheric Correction，QUAC)算法下的典型地物光譜出現失真現象。 商子健等[6]以歸一化植被指數(Normalized Difference Vegetation Index，NDVI)為指標，驗證了6S 大氣校正模型在GF－1 多光譜影像水體反演方面的適用性，其效果優于FLAASH 大氣校正。 曾群等[7]將不同的波段組合因子運用到大氣校正模型中，發現對于渾濁、高動態特性的水體，FLAASH 優于QUAC。 潘岑岑等[8]利用Hyspex 高光譜數據的研究指出，基于統計學的QUAC 和經驗線性法兩種大氣校正方法效果不穩定，導致其系數變化較大。 崔成嶺等[9]指出6S 大氣校正查找表精度較低的現狀，通過實時創建查找表的方式對資源1 號02D(ZY－1 02D)高光譜影像做大氣校正處理，總體精度與FLAASH 相近，且發現6S 對藍光波段的校正不完全。 劉瑤等[10]測試了資源1 號02D 高光譜影像在內陸水體葉綠素a 濃度反演方面的適用性，其波段比值模型取得了最好的效果，并指出針對于ZY－1 02D 水體圖像的降噪和大氣校正方法是未來的發展方向。 高光譜影像的大氣校正在針對不同地物的反演中亦表現出不同的效果[11－12]，了解地物反演的參數貢獻有利于提高大氣校正的精度和效率。 目前，針對高光譜影像在二類水體葉綠素a 濃度反演方面的大氣校正方法相關研究較少，由于波段連續細微的高光譜數據對大氣干擾極為敏感，需要更精準的大氣校正處理，才能滿足水體葉綠素a 濃度反演的數據需求。

文章結合內陸水環境研究對影像光譜高精度的需求及高光譜衛星遙感器的特點，以南四湖為研究區，以(ZY－1 02D)高光譜遙感影像為數據源，按照6S、FLAASH、QUAC 大氣校正的特點，從影像外部大氣產品校正、影像內部大氣補償參數校正和影像的光譜均值校正3 個角度分別對影像進行大氣校正，進而提取多維光譜指數，并利用半分析和機器學習方法驗證反演模型的反演精度。

1 數據源與實驗區選取

1.1 研究區和水體樣本采集

南四湖位于微山縣境內，由微山、昭陽、獨山和南陽等4 個湖泊南北相連組成，是南水北調東線工程重要的水源地[13－14]，湖內水生生物與水禽種類眾多[15]，其水質狀況對東線工程以及生物多樣性影響巨大。 南四湖面積約為1 266 km2，平均水深為1.5 m[16]，其流域多年平均年降水量在700 mm 以上，入湖河流有50 多條[17]，且來水河流近九成集中于上級湖，外流入湖和湖內運作均對南四湖水質產生不同程度影響[18]。

根據衛星過境時間，按照均勻布點原則，在獨山湖周邊采集38 個表層水樣，分布如圖1 所示。 文章所用地圖審圖號為魯SG(2023)031 號。

圖1 南四湖位置及實測點分布示意圖

1.2 數據源及其預處理

遙感影像數據為2021 年9 月12 日的ZY－1 02D 高光譜數據，其搭載的可見短波紅外高光譜相機(Advanced Hyperspectral Imager， AHSI)傳感器比高分5 號的信噪比更高，但波段數量減少到166 個，光譜分辨率也隨之降低[19]。 原始影像數據經輻射定標處理，得到大氣校正需要的表觀反射率和輻射亮度值。

外業數據采集利用美國ASD 公司生產的FieldSpec 4 Hi－Res 便攜式地物光譜儀，按照白板、水體、天空、白板的順序，實測采樣點水體高光譜數據，每點重復5 次。 瓶裝湖面表層水樣[20]，放入黑色不透光袋中保存，并用GPS 定位儀記錄采集點位置信息。 返回室內避光處， 采用美國安諾牌ChloroTech 121A 手持式葉綠素測定儀測定水體葉綠素a 濃度。

2 理論和方法

2.1 大氣校正模型

大氣校正可以保留由影像地物重要組成成分差異導致的反射率的微小差別信息[8]。 針對不同地物的光譜異質性及周邊環境的差異，衍生出的大氣校正方法也具有各自的側重性。 6S 模型側重影像過境時的觀測及大氣條件參數，并基于此參數建立的查找表進行逐像元插值得到監測點的大氣校正系數；FLAASH 模型則側重于影像像素光譜特征，考慮周邊像元對目標像元造成的“鄰近效應”、光譜噪聲等因素對目標地物校正；QUAC 模型則是整幅影像的光譜信息收集，監測點的光譜特征即為目標地物的光譜均值。

6S 輻射傳輸模型由5S 模型發展而來[21]，波段處理范圍為0.25 ～4.0 μm，通過模擬信號傳輸過程的太陽輻射，并計算信號在進入傳感器前的輻射能量，得到校正參數進行大氣校正[7]。 ZY－1 02D 高光譜影像需進行可見光近紅外波段(Visible and Near-Infrared，VNIR)和短波紅外波段(Short Wave Infrared，SWIR)的數據整合，將表觀輻亮度轉換為表觀反射率，由式(1)表示為

式中ρTOA為表觀反射率； Lλ為波段λ 處的中心波長，nm；d 為日地距離，AU；Eλ為波段λ 處的太陽光譜輻照度，W/(m2·μm)；θ 為太陽天頂角，(°)。

6S 模型通過輸入表觀反射率和影像元文件，得到各波段大氣校正系數xa、xb、xc。 將表觀反射率轉化為真實的地表反射率ρr，可由式(2)和(3)表示為

FLAASH 模型基于MODTRAN4 ＋輻射傳輸模型，波段區間0.4 ～2.5 μm。 通過設定模型參數，逐像元反演校正參數。 針對不同的波段區間，進行特定的水汽反演，采用暗目標法進行氣溶膠厚度的反演。 影像像元光譜輻射亮度由式(4)表示為

式中L 為像元輻射亮度；ρ 和ρe分別為像元與相鄰像元地表反射率均值；S 為大氣球面反照率；Lα為大氣程輻射； A、B 是依賴于大氣(透過率) 和幾何狀況的系數。S、Lα、A、B 可由輻射傳輸模型MODTRAN4＋計算得到。

QUAC 是基于影像本身的統計方法，采集像元內的地物光譜值，取同一地物的光譜平均值作為判別地物的經驗值，利用端元平均反射率與參考物進行對比確定大氣的影響。 QUAC 大氣校正算法并不依賴影像獲取過程中的各類參數信息，大氣參數和儀器標定的誤差對校正結果的影響較小[7]。 端元平均反射率ρ′由式(5)表示為

式中ρ1，ρ2，…，ρn為視場內各物質的端元光譜反射率；n 為端元個數。

2.2 水體葉綠素a 濃度反演

(1) 三波段光譜指數

光譜指數是一種基于光譜反射率構建的指標函數，依據地物的光譜特性，通過波段組合度量地表參量[22]。 在嘗試了多種三波段組合后，選擇遍歷后相關系數較高的公式反演三波段指數(Three-Band Index，TBI)，由式(6)表示為

式中Rλ1、Rλ2、Rλ3分別為395 ～900 nm 范圍內波長為λ1、λ2和λ3處的遙感反射率。

(2) 四波段光譜指數

按照三波段光譜指數法選取波段組合的方法，運用水體反演葉綠素a 濃度的四波段公式反演四波段指數(Four-Band Index，FBI)，由式(7)[23]表示為

式中Rλ4為395～900 nm 范圍內波長為λ4處的遙感反射率。

(3) CatBoost 算法

CatBoost 算法是梯度提升決策樹(Gradient Boosting Decision Tree，GBDT)算法框架的改進算法，運用排序提升方法構建模型，在不同的迭代中會采用不同的排列順序，訓練集中的噪聲點被消除，梯度估計與濃度預測誤差得到了解決。 CatBoost 算法通過減少超參數調優來抵抗模型過度擬合，增強了算法的準確性和泛化能力，其計算過程添加了先驗值和權重參數，一些低頻率類型的數據和噪聲點對數據整體趨勢的影響得到有效控制。 CatBoost 算法可由式(8)表示為

式中xσj，k、xσp，k分別為第k 個訓練樣本的第j、p 個類別特征；xi，k為類別特征平均值；[]為指示函數運算，即括號內兩個量相等時取1，否則取0；Yj為第j 個樣本的標簽； P 為添加的先驗項；a 為先驗值的權重。

(4) 反演模型構建

對大氣校正的影像提取像元光譜，并基于葉綠素a 實測濃度，提取最優TBI 和FBI 用于CatBoost模型的特征組建。 CatBoost 模型采取相同的隨機種子進行訓練集(70%)與驗證集(30%)的劃分，網格搜索法在指定范圍內進行枚舉[24]，將效果最好的參數用于模型構建。

一元線性模型的建模參數分別為TBI 和FBI 的最優遍歷結果[25]，訓練集與驗證集的劃分同CatBoost 模型。

2.3 大氣校正及葉綠素a 反演精度評價

為定量比較大氣校正模型的校正效果，使用光譜角制圖(Spectral Angle Mapper，SAM)[26]、均方根誤差(Root Mean Square Error，RMSE)和相關系數來衡量不同模型校正后的影像光譜反射率與實測反射率之間的接近程度[21]。 SAM 光譜角α 的余弦值由式(9)表示為

式中N 為樣本總數；Ai、Bi分別為第i 個像元向量的光譜值。

相關系數r 由式(10)表示為

式中Xi、Yi分別為像元光譜值和實測光譜值；和分別為像元、實測的光譜均值。

葉綠素a 濃度反演精度采用決定系數R2和均方根誤差RMSE 評價。 RMSE 和R2分別由式(11)和(12)表示為

式中Xai、Yai分別為葉綠素a 濃度的實測值和反演值；yi、分別為實測值和模型反演值；為N 個yi的均值。

3 結果與分析

3.1 大氣校正

3.1.1 大氣校正下的光譜曲線對比

由于實測點位分布在獨山湖區域，因此僅針對南四湖的獨山湖區域展開分析。 不同大氣校正參數見表1，觀測日期為2021 年10 月22 日，過境時間為世界標準時間(Universal Time Coordinated， UTC)03:19:42，中心經緯度為35.28°N、116.80°E。 其中，FLAASH、QUAC 大氣校正參數均在ENVI 中完成，由于沒有相應的傳感器類型，所以統一選擇“unknown sensor”。

表1 大氣校正參數表

為保證光譜信息的可信度，任取一點(14 號點)的實測值和全部監測點光譜均值，分別與大氣校正光譜對比，結果如圖2 所示。 由于葉綠素和類胡蘿卜素的吸收，在440 和490 nm 附近形成兩個小反射谷；藻類色素的低吸收以及水中懸浮物的散射，使得光譜曲線在570 nm 附近形成了大反射峰；4 種光譜曲線均在675 和700 nm 附近分別形成了大反射谷和大反射峰，這與前人的研究結果一致[27]。 對比圖2(a)和(b)發現，除反射率略有升高外，光譜均值與14 號點光譜曲線的整體趨勢基本一致。 由于受到水面天空光等多種信號的影響，大氣校正曲線普遍高于實測光譜曲線。 校正方式的不同使得光譜反射率的值高于實測值的程度也不一樣，6S 借助大氣校正產品逐一計算得到每波段每像素的地表反射率值，光譜反射率最接近實測值；而FLAASH 是基于整張圖像[4]，且重采樣減少了光譜的波段數量，對于藍光波段的不完全校正使得FLAASH 曲線出現部分負值。 QUAC 則是基于影像本身，且水體占比較大，有利于進行水體的光譜校正。

圖2 不同大氣校正光譜對比

3.1.2 精度評價

38 個采樣點的實測光譜數據與大氣校正數據的指標值如圖3 所示。 6S 和QUAC 模型的光譜角余弦值都在0.9 以上，而FLAASH 模型的余弦值略低一些，且各點之間的差異較大。 同時6S 模型的r在3 種大氣校正模型中最高，這也說明6S 模型對地物光譜信息的保持度高。 6S 與FLAASH 模型的RMSE 值均約為0.15，QUAC 的值與實測值的離散程度較大。 綜上所述，6S 模型的評估結果均表現良好，FLAASH 模型的光譜信息保持較弱，QUAC 模型校正后的影像反射率缺乏代表性，校正效果不穩定。

圖3 大氣校正精度指標圖

3.2 葉綠素a 濃度反演

3.2.1 多維光譜指數提取

在文章研究的光譜波段范圍內(395～900 nm)，通過MATLAB 軟件實現基于TBI 和FBI 公式的波段反射率循環迭代。 基于實測葉綠素a 濃度，利用相關系數r 提取最優的波段組合，計算模型參數。選取相關性最高的前5 個光譜指數組合作為CatBoost 反演模型的參數，選取相關性最優的光譜指數進行一元線性回歸模型的反演。 TBI 和FBI 的相關結果分別見表2 和3。

表2 三波段光譜指數及相關系數

表3 四波段光譜指數及相關系數

3.2.2 葉綠素a 濃度反演

對影像光譜數據和實測葉綠素a 濃度，數據測試集和訓練集分別占30%、70%。 將基于TBI 和FBI 構建的一元線性反演模型應用到遙感影像上，分別得到12 個點的反演值。 將不同算法的反演值與實測值進行擬合，結果如圖4 所示。 6S 大氣校正后得到的葉綠素a 濃度反演值，在以三波段和四波段作為特征參數的線性模型中，均具有良好表現。QUAC 大氣校正與FLAASH 大氣校正后的線性模型反演精度略有差異，其中QUAC 四波段線性模型精度最高。 整體來說，反演值與實測值的R2最高能達到0.74，模型效果較好。

圖4 不同大氣校正算法反演結果與實測結果擬合圖

選擇與一元線性模型相同的訓練集，由于CatBoost 不需要過多調參，因此主要調節模型決策樹的數量(n＿estimators)、學習率(learning＿rate)和決策樹的深度(depth)，其余參數保持默認，按照參數調優標準，得到CatBoost 模型參數見表4。

表4 CatBoost 模型反演參數設置

AHSI 影像反演模型和精度評價結果見表5。CatBoost 模型的建模精度總體良好，測試集最高相關系數R2＝0.80、RMSE ＝4.97 μg/L。 由于CatBoost 模型在不同參數組合中取得的精度最高，且在6S 和QUAC 大氣校正后的參數反演結果均較好，因此選用CatBoost 模型進行南四湖水體葉綠素a 濃度反演。

表5 AHSI 影像反演模型和精度評價結果

對比不同大氣校正的模型建模結果，6S 大氣校正的模型擬合效果最好，其中四波段參數CatBoost模型的擬合精度最高。 基于QUAC 大氣校正的三波段參數和四波段參數組合的CatBoost 模型擬合精度均優于FLAASH 大氣校正。 QUAC 大氣校正效果缺乏穩定性，模型最低擬合精度R2＝0.35。

將模型應用到影像，得到南四湖葉綠素a 濃度反演結果如圖5 所示。 可以看出，質量濃度空間差異明顯，高值主要分布在南陽湖的東部沿岸和獨山湖的湖心，流動性較弱區域，這與前人的研究結果相符[28]。 然而，該組合方法下模型的泛化能力較差，選取的波段對該地區的水體葉綠素a 缺乏敏感性，影像模糊是由于反射信號包含高光譜噪聲。

圖5 不同大氣校正下的CatBoost 模型AHSI 影像南四湖葉綠素a 濃度

對比不同光譜指數的反演，四波段參數選擇的波段反射率更接近實測反射率，參數信息更能體現葉綠素a 濃度變化。 這一點在6S 與QUAC 大氣校正的影像中表現明顯。

6S 大氣校正的三波段與四波段參數模型反演，R2均穩定在較高水平，表明6S 大氣校正對于高光譜反演水體葉綠素a 濃度具有良好的適用性，可作為影像預處理時的優先選擇大氣校正方法。

考慮到QUAC 模型以影像本身作為對象，對同類地物光譜采集取均值的校正特性，對比葉綠素a濃度反演結果分析，特征影像受地物異物同譜特性的影響較大。 由于受大氣等光學條件的影響，不同水體葉綠素a 濃度監測點的光譜可能存在較大差異[5]，在求取均值后，不同葉綠素a 濃度對應的光譜值存在相近或相等的情況，這將對反演造成困難，三波段和四波段參數的光譜值出現誤差，造成最終反演精度出現不確定性。

FLAASH 大氣校正后的反演效果均不理想，可能是選取的TBI 與FBI 所用波段對葉綠素a 并不敏感[29]，FLAASH 模型以大氣校正產品作為系數，不同波段參數的組合對影像的泛化能力也不盡相同，導致構建的模型不能很好地預測水體中的葉綠素a濃度。

4 結論

基于南四湖ZY－1 02D 影像，分別進行了6S、FLAASH 和QUAC 的大氣校正處理，并對提取的光譜反射率進行多種組合，分別利用線性回歸模型和CatBoost 模型中反演水體中葉綠素a 濃度，得到主要結論如下:

(1) 6S 模型校正效果最佳，QUAC 次之，FLAASH 模型效果最差；

(2) CatBoost 模型能在一定程度上提高反演精度，排序提升的模型構建以及先驗值的添加消除了誤差噪聲，在不需要過多調參的情況下即可達到較好的回歸預測精度；

(3) 6S 模型的四波段組合CatBoost 模型反演結果的R2最高為0.80，更有利于提高光譜數據在南四湖中部的獨山湖葉綠素a 濃度反演中的真實性和反演精度。