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

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

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

0 引言

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

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

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

1 數(shù)據(jù)源與實驗區(qū)選取

1.1 研究區(qū)和水體樣本采集

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

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

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

1.2 數(shù)據(jù)源及其預(yù)處理

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

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

2 理論和方法

2.1 大氣校正模型

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

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

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

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

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

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

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

式中ρ1，ρ2，…，ρn為視場內(nèi)各物質(zhì)的端元光譜反射率；n 為端元個數(shù)。

2.2 水體葉綠素a 濃度反演

(1) 三波段光譜指數(shù)

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

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

(2) 四波段光譜指數(shù)

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

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

(3) CatBoost 算法

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

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

(4) 反演模型構(gòu)建

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

一元線性模型的建模參數(shù)分別為TBI 和FBI 的最優(yōu)遍歷結(jié)果[25]，訓(xùn)練集與驗證集的劃分同CatBoost 模型。

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

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

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

相關(guān)系數(shù)r 由式(10)表示為

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

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

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

3 結(jié)果與分析

3.1 大氣校正

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

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

表1 大氣校正參數(shù)表

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

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

3.1.2 精度評價

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

圖3 大氣校正精度指標(biāo)圖

3.2 葉綠素a 濃度反演

3.2.1 多維光譜指數(shù)提取

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

表2 三波段光譜指數(shù)及相關(guān)系數(shù)

表3 四波段光譜指數(shù)及相關(guān)系數(shù)

3.2.2 葉綠素a 濃度反演

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

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

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

表4 CatBoost 模型反演參數(shù)設(shè)置

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

(1) 6S 模型校正效果最佳，QUAC 次之，F(xiàn)LAASH 模型效果最差；

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

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