基于TM影像的博斯騰湖葉綠素a濃度反演研究*

牛 婷，文 方*，常夢迪，高 潔，李蘊輝，程 帆，杜姿影，陳曉斐

（1.新疆維吾爾自治區(qū)環(huán)境保護(hù)科學(xué)研究院，新疆 烏魯木齊 830011；2．新疆環(huán)境污染監(jiān)控與風(fēng)險預(yù)警重點實驗室，新疆 烏魯木齊 830011；3．新疆清潔生產(chǎn)工程技術(shù)研究中心，新疆 烏魯木齊 830011；4．國家環(huán)境保護(hù)準(zhǔn)噶爾荒漠綠洲交錯區(qū)科學(xué)觀測研究站，新疆 烏魯木齊 830011；5．成都信息工程大學(xué) 資源環(huán)境學(xué)院，四川 成都 610225）

富含大量有機污染物的工業(yè)、農(nóng)業(yè)以及居民生活廢水排入湖水中，造成湖泊水質(zhì)嚴(yán)重富營養(yǎng)化，從而導(dǎo)致藻類以及其他浮游生物的大量爆發(fā)，形成水華，而葉綠素作為藻類的主要組成部分，且葉綠素濃度便于人工測定。因此葉綠素濃度常常作為反映水體富營養(yǎng)化程度的重要評價指標(biāo)[1]。常規(guī)的水質(zhì)監(jiān)測是通過采集水樣進(jìn)行實驗室分析[2]，比較費時費力，對大面積水體的水質(zhì)狀況難以全面反應(yīng)，遙感技術(shù)可以實現(xiàn)湖泊水質(zhì)的高頻、大范圍、準(zhǔn)實時監(jiān)測[3]。葉綠素主要分為葉綠素a和葉綠素b，其中葉綠素a又是藻類植物中最主要的色素，所以葉綠素a濃度是水色遙感監(jiān)測中最重要的參數(shù)之一。利用遙感數(shù)據(jù)反演葉綠素的濃度主要有三種模型，即經(jīng)驗?zāi)Ｐ?、半分析模型和分析模型，不同的反演模型因為水體組成成分的不同而使得反演精度也不同。經(jīng)驗?zāi)Ｐ突诮y(tǒng)計分析，算法比較成熟、過程簡單，但缺乏物理依據(jù)，反演出的模型精度較低；半分析模型通過水體組成成分的光譜特征與統(tǒng)計分析相結(jié)合，有一定的物理依據(jù)，模型精度較高，外推性適宜性較好；分析模型基于大量的地面實時采集數(shù)據(jù)，物理機理復(fù)雜，模型精度高，但模型建立難度大。

Gordon在1975年提出了一種新型模型稱為分析模型，它主要依賴于水體組分、固有光學(xué)量、表觀光學(xué)量之間的關(guān)系，通過代數(shù)方程直接求解葉綠素濃度，具有極高的物理性與相關(guān)性。但是相較于經(jīng)驗/半經(jīng)驗?zāi)Ｐ头?，分析模型需要大量水體各種組分的固有光學(xué)特性數(shù)據(jù)，具體實施較為困難，建立算法的難度比較大[4]。劉忠華[5]、RUNDQUIST等[6]通過大量的實測數(shù)據(jù)，表明葉綠素a濃度在特定波長處于地表反射率具有較高的相關(guān)性。祝令亞[7]和溫新龍等[8]以太湖為例采用波段組合算法，分別建立了葉綠素a與MODIS數(shù)據(jù)、環(huán)境一號衛(wèi)星CCD數(shù)據(jù)的反演模型。李旭文等[9]基于Landsat5 TM數(shù)據(jù)和地表實測數(shù)據(jù)建立了經(jīng)驗?zāi)Ｐ?，對梅梁湖區(qū)藍(lán)藻生物量進(jìn)行了估算，證明葉綠素a濃度和差異植被指數(shù)的相關(guān)性較高，同樣也有諸多研究表明TM不同波段組合與葉綠素濃度具有較好的相關(guān)性[10-16]。

博斯騰湖位于內(nèi)陸河口，是焉耆盆地各種河流的聚集地，受人類活動影響，屬于典型的二類水體。本文通過對博斯騰湖水體葉綠素a濃度的反演，研究博斯騰湖中水體葉綠素的濃度分布，從而更好地應(yīng)對水體富營養(yǎng)化對博斯騰湖水資源的影響，解決博斯騰湖水體富營養(yǎng)化問題。基于以上分析，本課題選用經(jīng)驗?zāi)Ｐ秃桶敕治瞿Ｐ蛯Σ┧跪v湖水體葉綠素a濃度進(jìn)行反演建模研究。

1 研究區(qū)概況

博斯騰湖（Bosten Lake），維吾爾語意為“綠洲”，位于中國新疆維吾爾族自治區(qū)焉耆盆地東南面博湖縣境內(nèi)，是中國最大的內(nèi)陸淡水吞吐湖（見圖1）。湖區(qū)介于東經(jīng)86°19′-87°28′，北緯41°46′-42°08′之間，湖面的平均海拔高度約為1048米，總面積1646平方公里，水域面積約為800多平方公里。其既是位于其上游區(qū)域開都河、黃水溝等水系的尾閭，又是位于其下游區(qū)域孔雀河水系的源頭，相當(dāng)于一個巨大的“調(diào)節(jié)水庫”，是開孔河流域的“心臟”。湖區(qū)南邊和北邊海拔較高，受湖泊水位變化影響較小，西側(cè)和東側(cè)地勢平坦更容易受到湖泊水位變化的影響，周圍生長著廣茂的蘆葦，是中國四大集中產(chǎn)葦區(qū)之一。

圖1 博斯騰湖地理位置圖

2 數(shù)據(jù)與方法

2.1 水體葉綠素濃度實測數(shù)據(jù)

本研究采樣點的布設(shè)沿用巴州環(huán)境監(jiān)測站采樣點分布（見圖2），采樣點位坐標(biāo)見表1。對照影像數(shù)據(jù)，選用2010年6月9日、2010年9月1日采樣點實測數(shù)據(jù)進(jìn)行反演研究。

表1 博斯騰湖水質(zhì)監(jiān)測點位坐標(biāo)

圖2 博斯騰湖水環(huán)境監(jiān)測采樣點分布圖

2.2 TM數(shù)據(jù)預(yù)處理

TM成像儀是Landsat4和Landsat5攜帶的傳感器，它是一種改進(jìn)的多光譜掃描儀。為保證反演的精度，選擇天氣晴朗、研究區(qū)上方天空少云或者無云及衛(wèi)星過境時間最靠近實測數(shù)據(jù)時間的TM影像。根據(jù)地面水質(zhì)監(jiān)測數(shù)據(jù)，選取2010年6月10日、2010年8月29日的影像進(jìn)行建模。

遙感器觀測目標(biāo)物的輻射或者反射的電磁能量時，不可避免地受到大氣分子、氣溶膠和云粒子等大氣成分的吸收與散射的影響，再加之遙感器本身的光電系統(tǒng)特征、地形以及太陽高度等因素，所得到的測量值與目標(biāo)物真實的光譜轄射亮度或者反射率等是有一定差異的。地表參數(shù)的遙感定量反演需要糾正目標(biāo)輻射的不確定性信息。

對遙感影像數(shù)據(jù)的預(yù)處理主要包括幾何校正、輻射定標(biāo)、大氣校正。本研究針對這三個方面的要求對獲取的遙感影像進(jìn)行相應(yīng)的預(yù)處理。

2.2.1 TM影像的幾何精校正

利用已知坐標(biāo)的遙感圖像數(shù)據(jù)，選擇適當(dāng)?shù)目刂泣c，然后基于多項式糾正模型的校正。本研究以谷歌影像數(shù)據(jù)為底圖，在ENVI軟件中對遙感影像進(jìn)行幾何精糾正，糾正誤差控制在0.5個像元以內(nèi)。

2.2.2 TM影像的輻射定標(biāo)

本研究基于ENVI軟件，利用ENVI的Landsat Calibration模塊，讀入TM數(shù)據(jù)各個波段的定標(biāo)系數(shù)，將原始數(shù)據(jù)DN值轉(zhuǎn)換為大氣外層表面反射率（輻射亮度值）。

2.2.3 TM影像的大氣校正

大氣校正就是將轄射亮度或者表觀反射率轉(zhuǎn)換為地表實際反射率，目的是消除大氣散射、吸收、反射引起的誤差。本文采用ENVI軟件的FLAASH大氣校正模塊對TM輻射定標(biāo)后的數(shù)據(jù)進(jìn)行大氣校正，計算影像的真實反射率。需要輸入的參數(shù)中幾何參數(shù)可以從影像數(shù)據(jù)的元數(shù)據(jù)獲得，包括影像中心經(jīng)緯度、太陽高度角、傳感器類型和遙感影像獲取具體時間；大氣組分參數(shù)包括能見度（氣溶膠光學(xué)厚度，由FLAASH根據(jù)暗目標(biāo)法，根2200nm暗目標(biāo)的反射率與660nm暗目標(biāo)的反射率之間的比值關(guān)系，反演出TM影像的一個平均的能見度數(shù)據(jù)）、氣溶膠模式和大氣模式，大氣模式和氣溶膠類型參數(shù)則根據(jù)博斯騰湖的地理位置（經(jīng)緯度）和遙感影像數(shù)據(jù)成像時間進(jìn)行設(shè)置。

2.2.4 采樣點光譜信息的提取

首先對需要建模的影像數(shù)據(jù)采取3*3窗口的滑動平均處理，基于IDL-ENVI軟件，通過編程將采樣點的經(jīng)緯度信息轉(zhuǎn)化為TM反射率數(shù)據(jù)的文件坐標(biāo)，從滑動平均后的影像中提取出對應(yīng)坐標(biāo)的柵格像元值，得到TM數(shù)據(jù)各波段對應(yīng)的采樣點處的光譜信息。

2.3 研究方法

2.3.1 APPEL模型法

APPEL模型是El－Alem等針對Modis波段通道范圍與三波段模型不完全匹配的問題而提出[17]。葉綠素a在綠光波段的光譜特征表現(xiàn)為高反射，根據(jù)此特征可以用來獲取最大的葉綠素信息量；有色可溶性有機物（CDOM）在藍(lán)光波段表現(xiàn)為強反射，因此可以使用藍(lán)光波段來消除CDOM的影響；同時根據(jù)水體懸浮物在紅光波段表現(xiàn)為反射的光譜特征，選用紅光波段來最小化懸浮物的影響；近紅外波段是葉綠素的敏感波段，可以通過近紅外波段來去除藍(lán)光波段的葉綠素信息、藍(lán)光波段的后向散射影響和紅光波段的葉綠素信息；因為在紅、近紅外波段，水體的光譜特征表現(xiàn)為強吸收，因此后向散射的影響可以忽略不計。

APPEL模型光譜指數(shù)：

以此光譜指數(shù)來構(gòu)建APPEL模型的一般回歸模型：

公式（1）中bN、bB和bR分別為近紅外、藍(lán)光以及紅光波段；公式（2）中A、B為回歸模型的相關(guān)系數(shù)，Ca為反演的葉綠素濃度。

2.3.2 經(jīng)驗?zāi)Ｐ头?/p>

經(jīng)驗?zāi)Ｐ椭饕且匀~綠素a濃度和遙感參數(shù)之間的統(tǒng)計關(guān)系為基礎(chǔ)來實現(xiàn)對水體葉綠素a濃度的遙感反演，是一種較為廣泛的葉綠素a濃度反演模型[18]。

利用葉綠素a濃度在TM數(shù)據(jù)波段范圍內(nèi)的光譜特征差異，采用不同的波段組合方式，以擴大葉綠素a吸收峰與反射峰之間的差異，從而最大化地獲取葉綠素a濃度信息量，用于提取葉綠素a濃度。

本研究通過對TM數(shù)據(jù)多個波段組合的反射率光譜值與博斯騰湖葉綠a濃度實測數(shù)據(jù)進(jìn)行Pearson相關(guān)性分析，選擇相關(guān)系數(shù)較高的波段組合方式，將最佳波段組合方式計算出的反射率光譜值與博斯騰湖葉綠素a濃度實測數(shù)據(jù)進(jìn)行一元線性回歸分析，最終得出反演葉綠素a濃度的經(jīng)驗?zāi)Ｐ汀?/p>

2.3.3 對比分析

根據(jù)經(jīng)驗?zāi)Ｐ秃虯PPEL模型反演出的葉綠素濃度模擬值分別與實測數(shù)據(jù)計算，得出模型反演精度的評價指標(biāo)均方根誤差（RMSE），綜合兩種模型的各方面信息，選擇最佳模型作為博斯騰湖水體葉綠素a濃度含量反演模型。

3 結(jié)果分析

3.1 葉綠素a濃度反演模型

3.1.1 經(jīng)驗?zāi)Ｐ头囱萑~綠素a濃度結(jié)果

為提升經(jīng)驗?zāi)Ｐ头囱菥?，全面研究TM數(shù)據(jù)波段光譜值與葉綠素a濃度的關(guān)系，本研究參考國內(nèi)外學(xué)者對基于TM影像的葉綠素a濃度經(jīng)驗?zāi)Ｐ头囱莸难芯砍晒x取TM影像中的TM1-TM4單波段及其各波段的線性組合，選擇采樣點數(shù)據(jù)中的23個采樣點數(shù)據(jù)（2010年6月10日的P1-P11點、2010年8月29日的P1-P12點）作為反演數(shù)據(jù)建立回歸模型，利用SPSS軟件分析其光譜值與博斯騰湖采樣點實測數(shù)據(jù)之間的PEARSON相關(guān)系數(shù)，統(tǒng)計結(jié)果見表2。

表2 波段光譜值與葉綠素a濃度之間的相關(guān)系數(shù)

由表2可以看出，單波段中相關(guān)系數(shù)最高的兩個波段分別為b3和b4，達(dá)到了0.717、0.715，波段組合中相關(guān)系數(shù)最高的三個組合分別為b3*b4、b3+b4和b2+b3，其相關(guān)系數(shù)分別為0.81、0.76和0.718?？梢钥闯龆嗖ǘ谓M合的相關(guān)系數(shù)高于單波段，選擇b3*b4、b3+b4和b2+b3這三個波段組合方式作為因子分別與實測數(shù)據(jù)通過SPSS軟件做線性回歸分析，建立3個葉綠素a濃度的反演模型如下：

模型一：

模型二：

模型三：三個公式中Y代表葉綠素a濃度，波段組合代表相應(yīng)波段組合的真實地表反射率值。

3.1.2 APPEL模型反演葉綠素a濃度結(jié)果

水體中的葉綠素a濃度、懸浮物、可溶性有機物（CDOM）以及后向散射對水體的光譜特征有著較大影響。APPEL模型是通過分析上述因子對水體光譜特征反射、吸收的相互關(guān)系得出。APPEL模型最早是由ELAlem等提出，原模型是根據(jù)葉綠素在近紅外波段表現(xiàn)為強反射來提取最大的葉綠素信息量，但是作者在實際實驗過程中發(fā)現(xiàn)，使用近紅外波段提取葉綠素最大信息量的反演模型精度并沒有使用綠光波段的反演模型精度高，其相關(guān)系數(shù)（表3）相差較大。

表3 APPEL模型光譜值與葉綠素a濃度之間的相關(guān)系數(shù)

研究表明，當(dāng)水體中的葉綠素a濃度低于100ug/L時水體反射率紅光波段高于近紅外波段，這種現(xiàn)象隨著葉綠素a濃度的降低而更加明顯。博斯騰湖中的葉綠素a濃度較低，綠光波段的反射率大于紅光波段和近紅外波段，所以用綠光波段能提取到更大的葉綠素信息，因此本研究根據(jù)葉綠素在綠光波段強反射的特征來提取葉綠素的最大信息量。

綠光波段APPEL模型光譜值：

公式（6）中，模型中bG、bN、bB和bR分別為綠光、近紅外、藍(lán)光以及紅光波段的反射率。根據(jù)綠光波段APPEL模型的光譜值與實測數(shù)據(jù)建立線性回歸模型如下：

APPEL模型：

公式（7）中，X為APPEL模型的光譜指數(shù)Fa，Yappel為博斯騰湖的葉綠素a濃度。

3.1.3 最佳反演模型選擇

本研究因為2010年6月10日P14點位處于陸地上，故剔除了2010年6月10日P14這一個異常點，所以用剩下的10個采樣點數(shù)據(jù)進(jìn)行模型的精度評價。將10個點的光譜值代入上述四個模型中得到葉綠素a濃度的模擬值，將實測值與模擬值進(jìn)行比較（如圖3）。通過圖3可以看出，APPEL模型驗證的擬合效果較好。

圖3 葉綠素a濃度反演模型與檢驗

四個模型中，三個經(jīng)驗?zāi)Ｐ偷姆匠叹€性擬合度較高，從表4中可以看出b3*b4波段組合模型的R2最大為0.6638，可是其均方根誤差較其他模型最大為3.1；估測值與實測值之間相關(guān)系數(shù)最大的兩個分別為APPEL模型和（b2+b3）模型，分別為0.8887和0.8889，幾乎相同；但是APPEL模型的均方根誤差最低為0.8883，（b2+b3）模型的均方根誤差為1.4514；且APPEL模型原理是基于水體中葉綠素a、有色可溶性有機物（CDOM）、懸浮物等光譜特征的半經(jīng)驗半分析模型，模型外推適宜性較強，因此綜合考慮選擇APPEL模型作為反演葉綠素a濃度的最佳模型。

表4 模型擬合效果對比

3.2 葉綠素a濃度反演結(jié)果

在可見光到近紅外波段范圍內(nèi)，水體的吸收強度逐漸增大，在近紅外波段的吸收強度非常大，反射率極低，而植被在近紅外和綠光波段的反射強度非常大，所以本研究采用綠光波段與近紅外波段的比值來最大化提取水體信息。

通過ENIV軟件，導(dǎo)入大氣校正后的反射率數(shù)據(jù)，反演歸一化水體指數(shù)。歸一化水體指數(shù)（NDWI）的反演公式為（8）：

公式（8）中B2、B4分別為TM影像的綠光波段和近紅外波段，根據(jù)水體條件指數(shù)設(shè)立閾值，其中閾值大于0則為水體。利用Arcgis的按屬性提取模塊提取出NDWI值大于0的像元作為博斯騰湖水體的范圍，提取結(jié)果如圖4所示。

圖4 博斯騰湖水域范圍示意圖

利用IDL編程建立APPEL數(shù)學(xué)模型，并根據(jù)APPEL模型計算葉綠素a濃度。在2010年6月-2010年10月選取博斯騰湖天空無云的TM遙感影像，反演其葉綠素a濃度，分析葉綠素a濃度的時空變化情況。反演2010年6月10日、2010年8月29日、2010年9月5日、2010年9月30日和2010年10月16日的博斯騰湖水體葉綠素濃度分布圖（圖5）。

圖5 2010年6月、8月、9月、10月葉綠素a濃度分布圖

葉綠素是藍(lán)藻、水華的主要組成成分，從葉綠素a濃度反演結(jié)果分布圖（如圖6）上可以看出，葉綠素a濃度的空間分布受到風(fēng)向、水流方向的影響，而呈現(xiàn)出條紋狀，所有時間段內(nèi)葉綠素a在湖區(qū)靠近岸邊的地方分布較大，湖中心的葉綠素a分布較少。同時表現(xiàn)出6月份至8月份隨著時間增長，北部靠近岸邊的區(qū)域葉綠素a濃度分布逐漸減少，而南部和東部靠近湖岸區(qū)域葉綠素a濃度分布逐漸增加，這種現(xiàn)象的可能原因是博斯騰湖位于焉耆盆地，地形總體趨勢是北高南低，由于水體的流動，葉綠素a濃度分布自北向南隨著時間逐漸增加。9月初博斯騰湖水體中葉綠素a濃度達(dá)到2010年6月至10月的頂點。到了10月初湖岸區(qū)域葉綠素分布逐漸消失，湖區(qū)葉綠素含量逐漸減少。該反演結(jié)果與博斯騰湖對應(yīng)時間的葉綠素a濃度實測數(shù)據(jù)濃度分布趨勢基本相同，能夠反映出博斯騰湖葉綠素a在同一時間下的空間分布狀況以及不同時間下的時空分布狀況。

圖6 博斯騰湖不同時期葉綠素a濃度平均值

4 結(jié)論

本研究基于葉綠素a濃度實測數(shù)據(jù)與同時期多期TM影像數(shù)據(jù)，分別通過經(jīng)驗?zāi)Ｐ秃虯PPEL模型建立反演值與實測值關(guān)系模型，經(jīng)過精度評價以及對比分析，選擇最優(yōu)模型用于博斯騰湖TM影像葉綠素a濃度的反演，得出如下結(jié)論：

（1）TM影像單波段與實測葉綠素a濃度相關(guān)性較高的為b3、b4波段，為0.72左右，但b2、b3、b4的多波段組合與實測葉綠素a濃度相關(guān)性高于單波段，尤其是b3、b4的乘法組合，達(dá)到0.82左右。

（2）在本研究中，采用APPEL模型在綠光波段反射率大于紅光波段和近紅外波段，因此采用綠光波段APPEL模型的光譜值建立與實測數(shù)據(jù)線性回歸模型。

（3）通過對比三個多波段組合經(jīng)驗?zāi)Ｐ团c一個APPEL模型發(fā)現(xiàn)，四個模型估測值與實測值之間相關(guān)系數(shù)均較高，在0.85-0.89之間，多波段組合經(jīng)驗?zāi)Ｐ停╞3*b4）R2（0.6638）最高，但均方根誤差（3.100）也最大，APPEL模型均方根誤差最小為0.8883，同時估測值與實測值之間相關(guān)系數(shù)也較高（0.8887），同時APPEL模型外推適宜性較強，因此綜合考慮選擇APPEL模型作為反演葉綠素a濃度的最佳模型。

（4）湖岸區(qū)域的葉綠素濃度較高，湖心區(qū)域較低；在6月至10月葉綠素由湖區(qū)北部向湖區(qū)東南部轉(zhuǎn)移，并伴隨著時間葉綠素濃度逐漸減小。