應(yīng)用混合像元分解提取胡楊覆蓋度信息1)

郭春蕾 解濰嘉 黃華國(guó)

(省部共建森林培育與保護(hù)教育部重點(diǎn)實(shí)驗(yàn)室(北京林業(yè)大學(xué))，北京，100083)

責(zé)任編輯:王廣建。

胡楊(Populus euphratica)屬楊柳科(Salicaceae)胡楊亞屬木本植物，是第三紀(jì)遺留下來(lái)的孑遺植物，生長(zhǎng)緩慢，異性葉，樹(shù)冠不規(guī)則。胡楊林在我國(guó)荒漠區(qū)有大面積天然分布［1］，主要集中于西北地區(qū)的新疆、內(nèi)蒙古自治區(qū)西部、甘肅的內(nèi)陸河流沿岸的沖積平原［2］。由于胡楊具有涵養(yǎng)水源、保持水土、防風(fēng)固沙、調(diào)節(jié)地方氣候的生態(tài)防護(hù)效益，已成為抵御風(fēng)沙、遏制沙化、維護(hù)區(qū)域生態(tài)平衡、保護(hù)生物多樣性和保障綠洲農(nóng)牧業(yè)生產(chǎn)的重要屏障。因其數(shù)量稀少，長(zhǎng)期以來(lái)受到有關(guān)國(guó)家和國(guó)際組織的廣泛關(guān)注，也成為很多學(xué)者的研究熱點(diǎn)［3］。

目前，遙感技術(shù)已應(yīng)用于自然災(zāi)害監(jiān)測(cè)、農(nóng)作物估產(chǎn)、森林調(diào)查、土地利用調(diào)查等方面，特別是在資源調(diào)查上應(yīng)用較為成熟，直接為政府決策提供科學(xué)依據(jù)［4］。植被覆蓋率作為反映地表信息的重要參數(shù)，一直是植被遙感領(lǐng)域的重要研究課題，但是對(duì)胡楊林這種極度干旱地區(qū)植被的覆蓋度研究很少。胡楊林的特殊性在于樹(shù)冠不規(guī)則、孤立木較多、沙土背景反射率高，樹(shù)冠反射信息較弱，因此，提取覆蓋度有一定難度，有關(guān)干旱地區(qū)植被覆蓋度的應(yīng)用研究相對(duì)薄弱。

中分辨率影像(如Landsat TM/ETM+)，由于數(shù)據(jù)免費(fèi)，分辨率適中，時(shí)空連續(xù)性強(qiáng)，質(zhì)量穩(wěn)定，研究基礎(chǔ)好，仍然是植被監(jiān)測(cè)最常用的數(shù)據(jù)。就覆蓋度監(jiān)測(cè)而言，混合像元分解(SU)，尤其是線性混合像元分解(LSU)方法應(yīng)用最廣泛，在干旱半干旱區(qū)［5］、城市及周邊等地區(qū)［6］得到了成功應(yīng)用。因此，本文將采用LSU 方法對(duì)額濟(jì)納胡楊林自然保護(hù)區(qū)的Landsat TM 影像，進(jìn)行以胡楊為主體的植被覆蓋度的提取研究，分析LSU 方法的適用性，為大尺度動(dòng)態(tài)監(jiān)測(cè)研究提供技術(shù)支撐。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古額濟(jì)納胡楊林國(guó)家級(jí)自然保護(hù)區(qū)。額濟(jì)納旗土地面積11.46 萬(wàn)km2，其中90%以上為戈壁、沙漠和低山丘陵等裸露地面。綠洲地帶僅占全旗總面積的5.6%，主要植被有胡楊、沙棗、檉柳、梭梭等。額濟(jì)納旗地表水資源匱乏，發(fā)源于祁連山的額濟(jì)納河，是全旗林農(nóng)牧業(yè)的生命線。

保護(hù)區(qū)位于額濟(jì)納旗的中心位置——額濟(jì)納綠洲，西鄰額濟(jì)納旗政府駐地達(dá)來(lái)呼布鎮(zhèn)，北臨居延海，地理坐標(biāo)為97°10'23″～103°1'11″E、39°52'～42°47'N，屬野生植物類(lèi)型自然保護(hù)區(qū)。2003年1月經(jīng)國(guó)務(wù)院批準(zhǔn)建立國(guó)家級(jí)自然保護(hù)區(qū)。保護(hù)區(qū)總面積26 253 hm2，其中核心區(qū)面積為8 774 hm2，緩沖區(qū)面積為10 018 hm2，實(shí)驗(yàn)區(qū)面積為7 461 hm2。額濟(jì)納胡楊林覆蓋面積約占45 萬(wàn)畝，是當(dāng)今世界僅存的三處天然河道胡楊林之一，是我國(guó)典型荒漠地區(qū)天然胡楊林的主要分布區(qū)之一，也是內(nèi)蒙古西部荒漠區(qū)唯一的喬木林區(qū)，是阻止巴丹吉林沙漠向北擴(kuò)散的重要屏障，是中國(guó)西部生態(tài)的天然寶庫(kù)?，F(xiàn)今由于地下水資源的減少、不適宜的生產(chǎn)方式、現(xiàn)有保護(hù)措施的無(wú)力、旅游資源的無(wú)序利用等原因，以胡楊林為主的綠洲生態(tài)體系正遭到嚴(yán)重破壞。額濟(jì)納胡楊林自然保護(hù)區(qū)示意圖(見(jiàn)圖1)。

圖1 胡楊林自然保護(hù)區(qū)示意圖(右側(cè)為部分胡楊照片)

1.2 遙感數(shù)據(jù)

主要的遙感數(shù)據(jù)(見(jiàn)圖2)。其中，2009年9月的Landsat TM5 數(shù)據(jù)(圖2a)用于開(kāi)展LSU 覆蓋度提取。圖2b 為同期的經(jīng)過(guò)全色和多光譜融合以及幾何糾正的Quick Bird 圖像，用于手工樹(shù)冠分割，提取覆蓋度的相對(duì)真值，以實(shí)現(xiàn)對(duì)LSU 方法的驗(yàn)證。

應(yīng)用ENVI 軟件對(duì)TM 數(shù)據(jù)進(jìn)行了輻射定標(biāo)和FLAASH 大氣校正。校正前后胡楊林的像元波譜曲線(見(jiàn)圖3)，可以看到校正后的植被波譜曲線更符合真實(shí)植被波譜規(guī)律。然后，以QuickBird 圖像為基準(zhǔn)，采用二次多項(xiàng)式和最鄰近點(diǎn)內(nèi)插法對(duì)TM 影像進(jìn)行了幾何精校正(偏差＜0.5 個(gè)像元)，重采樣分辨率為30 m。最后，應(yīng)用ENVI 軟件對(duì)圖像進(jìn)行MNF 變換處理，分離出了數(shù)據(jù)的噪聲。

2 基于混合像元分解的胡楊覆蓋度提取

2.1 線性光譜混合模型

線性光譜混合模型是指像元在某一光譜波段的反射率(或亮度值)是由構(gòu)成像元的終端單元的反射率(或亮度值)及其所占像元面積比例為權(quán)重系數(shù)的線性組合［7－11］，公式如下:

式中:Ri為第R 波段第i 像元的光譜反射率，fki為對(duì)應(yīng)于i 像元的第k 個(gè)基本組分所占的比例，Ck為第k 個(gè)基本組分在第R 波段的光譜反射率，εi為殘余誤差值(即光譜的非模型化部分，值越小分解結(jié)果越好)。根據(jù)公式，利用最小二乘法求得每個(gè)終端單元的豐度，計(jì)算的結(jié)果表現(xiàn)為各終端單元的份量值圖像和殘余誤差圖像。

2.2 終端端元的確定

選取合適的端元是混合像元分解成功的關(guān)鍵［8－9］。端元選取包括確定端元數(shù)量以及各端元光譜。理論上，只要端元數(shù)量m 小于等于L +1(L 表示波段數(shù))，線性方程組就可以求解。然而，由于波段間存在相關(guān)性，選取過(guò)多的端元雖然會(huì)使殘差減少［10］，但會(huì)導(dǎo)致分解結(jié)果產(chǎn)生更大的誤差。因此，在能夠描述一個(gè)場(chǎng)景內(nèi)光譜的大部分方差的前提下，越少的端元數(shù)量是越好的選擇［11］。在線性混合像元分解模型的一般應(yīng)用中，選擇3～4 個(gè)端元較為適宜。本研究選擇3 個(gè)端元，分別應(yīng)用兩種端元提 取方法。

圖2 研究區(qū)假彩色合成遙感圖像

圖3 大氣校正前后植被的波譜曲線

2.2.1 基于幾何頂點(diǎn)的端元提取

將MNF 變換結(jié)果的前面兩個(gè)相關(guān)性很小的波段，作為X、Y 軸構(gòu)成二維散點(diǎn)圖［12］。在理想情況下，散點(diǎn)圖是三角形狀，根據(jù)線性混合模型數(shù)學(xué)描述，純凈端元幾何位置分布在三角形的三個(gè)頂點(diǎn)，而三角形內(nèi)部的點(diǎn)則是這三個(gè)頂點(diǎn)的線性組合，也就是混合像元。根據(jù)這個(gè)原理，在二維散點(diǎn)圖上選擇端元波譜，根據(jù)三角形三個(gè)頂點(diǎn)對(duì)應(yīng)于影像上的地物類(lèi)型，目視識(shí)別出散點(diǎn)圖上各端元所對(duì)應(yīng)的影像上的地物類(lèi)型。三角形頂點(diǎn)紅色部分代表裸地，綠色部分代表植被區(qū)域，藍(lán)色部分代表水體(見(jiàn)圖4)。

2.2.2 基于PPI 的端元提取

該方法是借助純凈像元指數(shù)(PPI)和n 維可視化工具用于端元波譜收集。

首先，在MNF 變化的結(jié)果上計(jì)算純凈像元指數(shù)，得到像元PPI 空間分布圖。PPI 值越大表示入選純凈指數(shù)的次數(shù)越多，也表示該像元在多波段影像空間里表現(xiàn)越獨(dú)特，光譜也越純凈，像元地物也越單一。在PPI 計(jì)算結(jié)果基礎(chǔ)上選擇感興趣區(qū)域，將結(jié)果值大于50 的點(diǎn)作為感興趣區(qū)域輸出。然后，利用ENVI 的多維可視化工具(n －D Visualizer)把上述選擇的感興趣像元點(diǎn)投影到MNF 變換主成分空間，選擇MNF 變化的前3 個(gè)波段，構(gòu)造3 維散點(diǎn)圖。如圖5a 所示，在每個(gè)頂點(diǎn)選擇幾個(gè)像元作為該類(lèi)地物終端單元。同樣選擇三類(lèi)地物，輸出三類(lèi)地物的端元波譜曲線，如圖5b 所示。參考原始影像中各類(lèi)地物的波譜曲線，識(shí)別每條端元波譜曲線對(duì)應(yīng)的地物類(lèi)型，可知藍(lán)色代表植被波譜，紅色代表水體波譜，綠色代表裸地波譜。

圖4 基于散點(diǎn)圖的端元選取

圖5 基于PPI 的端元選取圖

2.3 線性混合像元分解

通過(guò)以上步驟獲得各端元地物的波譜曲線，將生成的波譜曲線文件在端元波譜收集器中打開(kāi)，選擇LSU 方法，即可生成三個(gè)端元的覆蓋度圖像。

2.4 精度檢驗(yàn)

以Quickbird 為基準(zhǔn)影像，根據(jù)胡楊林的覆蓋程度從低到高，設(shè)立20 個(gè)檢驗(yàn)樣本(見(jiàn)圖6)。每個(gè)樣本面積為3 600 m2，對(duì)應(yīng)2 ×2 個(gè)TM 像元和100 ×100 個(gè)Quickbird 像元。在每個(gè)Quickbird 樣本內(nèi)，進(jìn)行人工解譯，將樹(shù)冠分割，并統(tǒng)計(jì)分割后樹(shù)冠面積比例作為該樣本覆蓋度真值。將每個(gè)樣本對(duì)應(yīng)的2 ×2 個(gè)TM 像元覆蓋度取平均值，作為預(yù)測(cè)值。通過(guò)比較真值和預(yù)測(cè)值，來(lái)評(píng)價(jià)覆蓋度提取精度。

3 結(jié)果與分析

3.1 精度評(píng)價(jià)

兩種端元提取方法得到的覆蓋度驗(yàn)證圖(見(jiàn)圖7)。由圖7可知，基于幾何頂點(diǎn)的端元提取的植被覆蓋度精度較高(R2=0.893，RMSE=0.12)?；赑PI 方法得到的覆蓋度相關(guān)性顯著偏低(R2=0.402)。

3.2 植被指數(shù)與植被覆蓋度的關(guān)系

植被指數(shù)(VI)同混合像元分解一樣，是目前用于提取植被覆蓋度的主要方法之一。迄今為止，國(guó)內(nèi)外學(xué)者已經(jīng)研究發(fā)展了幾十種不同的植被指數(shù)，其中應(yīng)用最廣泛的有歸一化植被指數(shù)(NDVI)、考慮土壤背景干擾的土壤調(diào)節(jié)植被指數(shù)(SAVI)、修正的土壤調(diào)節(jié)植被指數(shù)(MSAVI)等。通過(guò)對(duì)干旱地區(qū)的各類(lèi)植被指數(shù)的運(yùn)算分析，得到反演干旱地區(qū)植被覆蓋度精度最高的是歸一化植被指數(shù)(NDVI)［13］。

20 個(gè)驗(yàn)證樣本的NDVI 和真值覆蓋度的比較(見(jiàn)圖8)。由圖8可知，NDVI 和覆蓋度雖然具有一定相關(guān)性，但是偏差較大(R2=0.574)。精度不如基于幾何頂點(diǎn)的混合像元分解方法(R2=0.893)。

圖6 精度檢驗(yàn)的樣本分布圖(底圖為T(mén)M 的NDVI 圖像)

圖7 TM 提取的覆蓋度驗(yàn)證

圖8 NDVI 與覆蓋度真值的關(guān)系

3.3 保護(hù)區(qū)植被覆蓋度分布

將額濟(jì)納胡楊林自然保護(hù)區(qū)矢量圖疊加到研究區(qū)域上，分別對(duì)研究區(qū)內(nèi)的核心區(qū)、實(shí)驗(yàn)區(qū)、緩沖區(qū)以及河流兩岸300 m 范圍內(nèi)區(qū)域進(jìn)行植被覆蓋度提取。

核心區(qū)是額濟(jì)納自然保護(hù)區(qū)的核心，是額濟(jì)納自然景觀精華所在;實(shí)驗(yàn)區(qū)內(nèi)大部分為次生胡楊林，又有耕地、牧區(qū)和旅游區(qū);緩沖區(qū)左部分處于河流分界線，右部分則是核心區(qū)與實(shí)驗(yàn)區(qū)過(guò)渡地段。通過(guò)計(jì)算得到核心區(qū)、實(shí)驗(yàn)區(qū)、緩沖區(qū)以及河流兩岸的植被覆蓋度均值分別為:0.279、0.343、0.305、0.278?？梢钥闯?，核心區(qū)和河流兩岸覆蓋度較為接近;實(shí)驗(yàn)區(qū)和緩沖區(qū)覆蓋度反而略高，這主要是來(lái)自農(nóng)田和草地覆蓋度的貢獻(xiàn)。

4 結(jié)論與討論

研究表明，采用混合像元分解方法可以從中等分辨率的TM 數(shù)據(jù)中較好地提取胡楊林覆蓋度(R2=0.893，RMSE=0.120)?；趲缀雾旤c(diǎn)的端元提取結(jié)果優(yōu)于基于PPI 的端元提取結(jié)果和基于植被指數(shù)的統(tǒng)計(jì)回歸結(jié)果。由于美國(guó)陸地衛(wèi)星Landsat 擁有近30 a 的歷史數(shù)據(jù)，2013年又發(fā)射并正常運(yùn)行的8 號(hào)衛(wèi)星，是唯一全球覆蓋、免費(fèi)提供、業(yè)務(wù)化運(yùn)行的數(shù)據(jù)，結(jié)合本文方法可以為胡楊群落覆蓋度動(dòng)態(tài)變化研究提供數(shù)據(jù)支持。當(dāng)然，本文方法理論上也可以應(yīng)用于其他復(fù)雜類(lèi)型的森林提取植被覆蓋度。

本文需要深入研究和探討的問(wèn)題:首先，為什么幾何頂點(diǎn)的端元提取結(jié)果優(yōu)于基于PPI 的端元提取結(jié)果?PPI 方法對(duì)像元純度要求高，精度依賴于端元選擇精度。在較為稀疏的胡楊林冠層中，提取得到的純像元很難去除林冠間隙中的土壤背景的影響。而幾何頂點(diǎn)方法中每個(gè)頂點(diǎn)一般對(duì)應(yīng)水體、裸地和植被，具有一定的物理含義。根據(jù)幾何光學(xué)模型的定義，在稀疏林冠中，光譜為光照樹(shù)冠、陰影樹(shù)冠、光照土壤和陰影土壤四個(gè)組分的加權(quán)和，林冠間隙的土壤是可以分離的。因此，幾何頂點(diǎn)方法得到的端元應(yīng)該是樹(shù)冠本身，其覆蓋度可視為光照樹(shù)冠與陰影樹(shù)冠垂直可視面積比之和。所以，幾何頂點(diǎn)方法含義較為明確，操作也簡(jiǎn)單，適合胡楊覆蓋度提取。

其次，NDVI 和覆蓋度的關(guān)系為什么不高?在沙土高反射背景和稀疏林冠的共同作用下，NDVI取值較小，對(duì)樹(shù)冠敏感性降低。相比之下，混合像元分解是基于光譜反射模型中像元內(nèi)不同組分豐度的最優(yōu)化求解，物理意義明確，受林冠密度影響較低。

最后，覆蓋度提取仍然存在一定的誤差(RMSE=0.120)。誤差主要來(lái)自于樣地選取、QB 圖像分割問(wèn)題以及混合像元分解算法的誤差三個(gè)方面。樣地選擇雖兼顧高中低三種情況，分布也較為均勻，具有代表性，但是人工選擇中的主觀因素仍會(huì)有一定影響。另外，人工解譯QB 圖像時(shí)，分割的胡楊樹(shù)冠可能會(huì)比實(shí)際的不規(guī)則形狀略微平滑，在統(tǒng)計(jì)面積上也會(huì)有一定的誤差。最后，在基于幾何頂點(diǎn)的端元選取中，每個(gè)頂點(diǎn)選取的平均范圍大小會(huì)影響到端元平均光譜，進(jìn)而影響像元分解誤差。

［1］ 張萬(wàn)儒，主編.中國(guó)造林樹(shù)種土壤條件［M］.北京:中國(guó)科學(xué)技術(shù)出版社，1998.

［2］ 王世績(jī)，陳炳號(hào)，李護(hù)群.胡楊林［M］.北京:中國(guó)環(huán)境科學(xué)出版社，1995.

［3］ 趙文智，常學(xué)禮，李秋艷.人工調(diào)水對(duì)額濟(jì)納胡楊荒漠河岸林繁殖的影響［J］.生態(tài)學(xué)報(bào)，2005，25(8):1987 －1993.

［4］ Hugh W，Murray J，Peter S.Land use mapping using satellite data of Magawatu［J］.New Zealand ITC Journal 1996(2):149 －156.

［5］ 張良培，李德仁，童慶禧，等.鄱陽(yáng)湖地區(qū)土壤、植被光譜混合模型的研究［J］.測(cè)繪學(xué)報(bào)，1997，26(1):72 －76.

［6］ 萬(wàn)軍，蔡運(yùn)龍.應(yīng)用線性光譜分離技術(shù)研究喀斯特地區(qū)土地覆被變化:以貴州省關(guān)嶺縣為例［J］.地理研究，2003，22(4):439－446.

［7］ 張熙川，趙英時(shí).應(yīng)用線性光譜混合模型快速評(píng)價(jià)土地退化的方法研究［J］.中國(guó)科學(xué)院研究生院學(xué)報(bào)，1999，16(2):169 －175.

［8］ Tompkins S，Mustard J F，Pieters C M，et al.Optimization of endmembers for spectral mixture analysis［J］.Remote Sensing of Environment，1997，59(3):472 －489.

［9］ Elmore A J，Mustard J F，Manning S J，et al.Quantifying vegetation change in semiarid environments:Precision and accuracy of spectral mixture analysis and the Normalized Difference Vegetation Index［J］.Remote Sensing of Environment，2000，73(1):87 －102.

［10］ Radeloff V C，Mladenoff D J，Boyce M S.Detecting jack pine budworm defoliation using spectral mixture analysis:separating effects from determinants［J］.Remote Sensing of Environment，1999，69(2):156 －169.

［11］ Sabol D E，Adams J B，Smith M O.Quantitative subpixel spectral detection of targets in multispectral images［J］.Journal of Geophysical Research，1992，97(E2):2659 －2672.

［12］ 李君，范文義，李明澤.線性混合光譜模型在多光譜遙感影像分類(lèi)中的應(yīng)用［J］.東北林業(yè)大學(xué)學(xué)報(bào)，2008，36(8):45 －48.

［13］ 郭玉川，何英，李霞.基于MODIS 的干旱區(qū)植被覆蓋度反演及植被指數(shù)優(yōu)選［J］.國(guó)土資源遙感，2011(2):115 －118.