分層變端元混合像元分解的新疆北部積雪分量制圖研究

劉艷，楊耘，李楊

（1．中國氣象局烏魯木齊沙漠氣象研究所，新疆 烏魯木齊830002；2．長安大學地質工程與測繪學院，陜西 西安710054）

具有不同波譜屬性的物質出現在同一個像元內時會出現波譜混合現象，這類像元便稱為混合像元［1］。因傳感器空間分辨率的限制和被探測目標的多樣性，遙感影像中往往會產生混合像元［2］。理論上，形成混合像元的原因主要有以下幾方面：1）單一成分物質的光譜、幾何結構及其在像元中的分布；2）大氣傳輸過程中的混合效應；3）傳感器本身的混合效應［3］。其中，1）是線性效應，指混合像元中各端元間相互獨立互不影響時，混合像元光譜是該像元內各端元光譜的線性疊加［2］；2）和3）為非線性效應，它是由于端元（純像元）間的散射傳輸路徑和遙感儀器的混合效應所引入的光譜非線性疊加［2］。大氣校正可以對2）進行修正，而傳感器本身的混合效應可以通過儀器校準、定標加以部分克服。根據混合像元光譜的產生機理，研究人員建立了許多分解模型，主要包括線性、非線性分解模型、模糊監督和神經網絡等。線性分解模型因物理含義明確、建模簡單，得到了廣泛應用［4－9］。

新疆北疆地區冬季多云天氣較多，再加上地形和植被的影響，導致該地區積雪分布不均，這使得高精度的積雪遙感監測難以實現。從多光譜、高分辨率影像（ETM＋、QuickBird等）中雖能獲取中、高空間分辨率的積雪覆蓋數據，但這類影像具有如下特點：商業遙感影像價格昂貴且影像覆蓋范圍不大，時間分辨率有限。而新疆地域廣闊，冬天地面積雪隨時間變化頻繁。因此，高空間分辨率影像難以在氣象業務服務中得到廣泛應用。相對于高空間分辨率影像，中、低空間分辨率影像（MODIS、FY－3）可全球免費獲取，時間分辨率高，易獲取大范圍晴空條件下的影像數據，便于在業務中使用。但是，因傳感器空間分辨率的限制，MODIS等中低分辨率影像存在大量的混合像元，需要對混合像元進行分解，才能提高制圖精度。有學者利用中等分辨率成像光譜儀MODIS積雪數據對新疆或天山南、北地區進行積雪、冰川和融雪徑流研究［10－13］，但提取積雪面積以像元為單位，考慮混合像元問題的研究甚少。陳曉娜等［14］以MOD02HKM為基礎，通過線性分解模型對新疆天山中段MODIS影像進行分解，從中提取積雪面積；延昊和張國平［15］從甚高分辨率輻射計（advanced very high resolution radiometer，AVHRR）數據反演積雪蓋度；Painter等［16］利用機載可見光／紅外成像光譜儀（airborne visible infrared imaging spectrometer，AVIRIS）影像結合DISORT模型并考慮了亞像元對雪粒徑反演的影響，發展了MEMSCAG（MODIS snow－cover area and gain size）模型同步反演雪粒徑，取得了較好效果。通過遙感分析積雪覆蓋、草場植被和成片森林等非立體空間目標時，由于端元間散射效應弱，線性效應在混合光譜的成因中起絕對主導作用［2］。但是，線性分解模型采用固定端元對整個影像進行解混，在實際應用中一個像元可能由2、3或更多個端元組成，采用固定端元的解混方法會產生端元分解過剩的問題，特別是新疆北疆這種土地覆蓋和積雪分布狀況。

多端元混合像元分解法 （multiple endmember spectral mixture analysis，MESMA）是傳統混合像元分解法（spectral mixture analysis，SMA）的一個擴展，它的基本思想是基于像元為每類地物選取多條光譜，并以此生成多個端元組合（每個端元組合由不同地物的光譜組成），然后對每個像元搜索最小二乘誤差最小的端元組合，進而反演出每個像元的端元分量［17］。目前，MESMA被廣泛應用于高光譜遙感影像城市土地覆蓋分類中［17－22］。但是，MESMA仍存在如下問題：用于解混的端元模型中若端元數過少，導致模型難以達到滿意的解譯結果；如果模型中包含的端元個數過多，測量光譜和模擬光譜間的輕微偏離通常分配給模型所用端元，實際上這樣的端元并不存在［21］。在此背景下，針對新疆北疆地區特有的地理及氣候條件，本文利用MODIS數據，提出了新的分層變端元混合像元分解（hierarchical dynamic endmember spectral mixture analysis，DESMA）方法，進行積雪分量的遙感制圖研究。

1 材料與方法

1．1 研究區和數據

本文所用遙感影像是2012年1月6日獲取的、覆蓋新疆北疆地區的MODIS 1B數據，空間分辨率為500m，包含7個波段（MODIS 1～7）。影像中有一定量的積雪覆蓋，適合積雪分量遙感制圖研究（圖1）。數據進行了FLAASH（fast line－of－sight atmospheric analysis of spectral hypercube）大氣校正、輻射校正和幾何糾正等預處理。同時，收集同期研究區HJ－1BCCD數據和MODIS 8d合成最大積雪覆蓋產品（MOD10A2）。

圖1 研究區域位置Fig．1 Location of the study area in Northern Xinjiang

1．2 分層變端元混合像元分解方法基本思想

該方法的基本思想是根據地類特征對影像進行逐級分類，并逐級逐類地解混，將前一級解混所得A類地物的空間分布作為下一級其子類分解時的空間限制，采用既定模型僅在A區域內對影像進行解混（圖2）。在每一級分解時均假定像元在某一波段的反射率是以構成該像元的端元（純像元）反射率及其所占像元面積比例（分量）為權重的線性組合。其中，線性混合像元分解模型的表達式［23］為：

冬季研究區內耕地、草場和水體等均勻地表會被積雪完全覆蓋，云杉、荒漠植被和積雪混合分布，部分區域存在裸露地表。針對地物這種分布特征，第1級時將影像分為積雪－植被和土壤兩大類。然后僅對積雪－植被區采用積雪和植被兩端元模型進行解混，再消除山體陰影的影響。最后在植被區實施植被和積雪兩端元模型分解。綜合上述3級分解結果，最終獲取積雪分量空間分布。因濾去了第1級中與積雪無關地類，端元數減少，降低了地物錯分性。

圖2 分層變端元混合像元分解法流程圖Fig．2 Workflow of hierarchical dynamic endmember spectral mixture analysis（DESMA）

1．3 分層變端元混合像元分解制圖

1．3．1 建立端元庫 研究區主要地物有云杉、荒漠植被、耕地、裸土、草場、水體等，冬季大部分地物被積雪覆蓋。因此，針對研究區地物特征可建立植被－積雪－土壤的端元模型。可選端元庫包含影像端元和參考端元兩種。影像端元由研究區MODIS影像獲取。結合野外雪深GPS記錄點，在ENVI圖3環境中實現MODIS R（6）G（2）B（1）假彩色合成積雪（藍色）、植被（墨綠色）、裸土（紅色）的感興趣區域選擇，每類地物對應1～3個像元的矩形區域，由ENVI／VIPER TOOL［21］實現影像端元庫建立，包含91條光譜，其中積雪光譜19條、植被光譜49條、土壤光譜23條。參考端元源自全波段地物光譜儀測量積雪、荒漠植被、裸土光譜數據。測量光譜經MODIS傳感器響應函數（spectral response function，SRF）［22］轉換后利用 ENVI／Spectral Library Builder建立參考端元庫，包含13條積雪、2條植被和4條土壤光譜。

圖3 MODIS R（6）G（2）B（1）假彩色合成圖Fig．3 False color composite image of MODIS of the study［Red（6）Green（2）Blue（1）］

1．3．2 初選端元 利用EAR（Endmember Average RMSE）計算A端元光譜集內某一條光譜用A端元內其他光譜來分解產生的殘差，EAR越低，表明這條光譜的代表性越好，很高則證明這條光譜可能是離群點、沒有代表性［22］。因此，計算可選端元庫中各端元的EAR，挑選EAR較小的光譜。結果選出14條積雪光譜，6條植被光譜和8條土壤光譜。EAR計算公式：

式中，i為端元光譜；j為參考光譜；n為參考光譜總數；RMSE為均方根誤差。

1．3．3 優選端元 EAR初選后的端元光譜雖能很好地模擬其所在類的其他端元光譜，能否很好地對影像中所屬地類進行分解，需進一步實驗分析［22］。以積雪端元為例，首先，對積雪端元庫中各條光譜的EAR進行升序排列，然后，采用EAR最低的端元1對整個影像進行2端元分解，得到該端元分解影像像元數，其后引入端元2，利用端元1和端元2對影像解混，依次類推，得到各端元光譜組合分解像元數（表1），判斷各組合中每條光譜分解像元數是否大于44022（即分解像元數占整個影像像元總數的20%）。

1．3．4 逐級逐類解混和制圖 為了解釋地物因衛星傳感器角度、地形起伏和其他一些陰影，所有模型都包含陰影端元，用來表征端元亮度變化，而陰影端元并非地物組成成分。因此，需要對陰影進行歸一化以消除陰影端元，即用每個像元的非陰影端元分量除以全部非陰影端元分量的總和，最后得到每個像元實際組成端元的分量。

2 結果與分析

2．1 典型端元庫

大氣校正后MODIS通道1存在部分像元反射率大于1的反射率超飽和現象和通道5存在壞數據的情況。因此，本文選用MODIS 2、3、4、6和7通道反射光譜數據構建端元光譜庫。依據上述端元庫建立原則，最終選取6條積雪（表1中1、7、9、10、13和14）、5條土壤和4條植被光譜建立了適合研究區的積雪－植被－土壤端元庫（圖4a，b，c）。

表1 積雪端元優選原則Table 1 Criteria used to select snow endmembers for DESMA library

圖4 優選積雪（a），土壤（b），植被（c）端元的光譜Fig．4 Spectra contained in optimal endmember library：snow（a），soil（b），vegetation（c）

2．2 分級分解結果

由于僅收集到研究區同期HJ－1BCCD數據，缺少可用的近紅外數據IRS，僅利用可見光波段采用非監督分類和手動校正方法，提取積雪覆蓋，同時，利用夏季HJ－1BCCD獲取植被分布圖（圖5），將其與分級結果進行比對。

2．2．1 第1級分解結果 第1級分解時，2端元或3端元分解模型對分量最大和最小值不做限制。首先，將包含6條積雪、5條土壤和4條植被光譜的端元輸入第1光譜庫中，分別進行積雪－陰影、土壤－陰影和植被－陰影的2端元分解，端元分解次數為15；將包含積雪和植被的端元光譜輸入第1光譜庫，土壤端元光譜輸入第2光譜庫，進行積雪－土壤－陰影和植被－土壤－陰影的3端元分解，端元分解次數為50。然后，比較2端元和3端元分解模型的RMSE，如果RMSE變化大于0．1，選擇3端元分解結果，否則選擇2端元。比較發現，2端元模型RMSE整體較小，表明選擇2端元即可以完成分解，分解像元數達100%。合并2端元模型生成的積雪－陰影和植被－陰影類，生成積雪－植被區類型圖（圖6），吉木乃、黑山頭和福海附近存在裸露地表，阿勒泰和富蘊山區部分山坡無積雪覆蓋，與HJ－1BCCD分類圖（圖7）整體分類結果一致。但是，由于地形起伏和陰影影響，陰坡積雪區域呈現出明顯的非積雪特性，因此，僅利用HJ－1BCCD可見光通道信息進行分類，這部分區域被分成了非積雪區，造成了漏分。

圖5 研究區HJ－1BCCD假彩色合成圖Fig．5 False color composite image of HJ－1BCCD

圖6 第1級：2－端元和3－端元模型解混分類圖Fig．6 Merged result of 1st layer derived from 2－endmember and 3－endmember model of DESMA

2．2．2 第2級分解結果 第1級分解后得到積雪－植被的空間分布，將其對整個MODIS影像進行掩膜運算，得到積雪－植被區的MODIS數據，將包含積雪和植被的端元光譜輸入第1光譜庫中，進行積雪－陰影和植被－陰影的2端元分解，模型分解像元最小和最大分量控制在（－0．05，1．05），端元分解次數為10（圖8），吉木乃和黑山頭附近存在荒漠植被，哈巴河、阿勒泰和富蘊中山帶有云杉分布和山體陰影，前山帶靠近古爾班通古特沙漠邊緣區域有荒漠植被，與 HJ－1BCCD R（4）G（3）B（2）假彩色合成圖（紅色區域為植被）（圖5）分類結果一致。發現在積雪和植被混合區出現漏分現象。

圖7 研究區HJ－1BCCD積雪覆蓋圖Fig．7 Snow－cover map of HJ－1BCCD of the study area

圖8 第2級：掩膜后積雪－植被區MODIS影像2端元分解分類圖Fig．8 Classification of 2nd layer derived from 2－endmember model of DESMA for MODIS image in snow and vegetation area

2．2．3 第3級分解結果 第2級分解后得到植被區主要植被類型有荒漠植被和天山云杉兩類，考慮山體陰影與云杉易錯分，引入山體陰影端元，對掩膜后的植被區MODIS數據進行植被－陰影和山體陰影－陰影2端元模型，模型分解像元最小和最大分量控制在（－0．05，1．05）且滿足RMSE≥0．0025，將山體陰影從植被區中去除（圖9）。

圖9 修正了部分山體陰影的第2級分類圖Fig．9 Classification result of level 2excluded part of the mountain shadows

修正了部分山體陰影的植被區，還會存在部分積雪，將積雪光譜輸入第1光譜庫，植被光譜輸入第2光譜庫，對植被區MODIS數據進行積雪－植被－陰影的3端元分解，端元分解次數為24。

2．3 積雪分量制圖與精度驗證

疊加第2級積雪分量圖和第3級山體陰影修正后植被區積雪－植被－陰影3端元分解歸一化后所得積雪分量，得到研究區積雪分量圖。2012年1月在北疆地區進行了雪深、積雪覆蓋度的野外觀測（圖10和表3），輔助記錄了樣點周圍1km場景。將積雪分量與29個實測點進行比對，其中25個點正確，分量準確率為87%（表3）。

同時，收集同期MODIS 8d合成最大積雪覆蓋產品（MOD10A2）與研究區積雪分量圖進行比對分析（圖11），結果發現，29個實測點上MOD10A2均監測為積雪（200表征像元被積雪覆蓋），但是MOD10A2無雪區面積（白色區域）比本文提出的分層變端元混合像元分解生成的積雪分量圖中無雪區域面積大很多（白色區域），MOD10A2中無雪區比例為8%，本文計算積雪分量圖中無雪區域占3%。冬季研究區內湖泊被積雪覆蓋，MOD10A2未對湖泊和湖冰進行積雪識別，占研究區面積比例為1%。同時，分層變端元混合像元分解法考慮區域地貌特征和積雪分布狀況，以積雪為核心地類，各分級分類時，與積雪相關的地類利用2端元或3端元模型進行解混，為下一級積雪分量的定量分析提供相對純凈的數據環境，在山區云杉、山體陰影和積雪以及荒漠植被、積雪交錯分布的地貌環境下，與MOD10A2僅從像元尺度上提供積雪和非積雪空間信息相比，較大程度地提高了積雪識別率，真實反映了積雪的空間分布特征。

3 結論與討論

針對研究區土地覆蓋類型，提出了基于最小EAR和少端元分解模型的MESMA方法來實現以積雪為核心的地類分級分類方法。該方法通過逐層刪除非積雪地類，最終挑選出積雪，其適用條件寬泛，積雪覆蓋區均能適用。實驗結果表明，變端元解混時選擇較少的端元個數更有利于提高積雪制圖精度。相對于MESMA［17－22］，本文采用的最小EAR和選擇2－端元分解模型構建典型端元庫的方法，可使積雪制圖精度高達87%，即選擇較少的端元個數更有利于提高。

表3 模型計算值與實測值比較Table 3 Comparison between the result from our model and field data

圖10 研究區積雪分量與積雪蓋度測量點位置分布圖Fig．10 Snow fraction map generated from DESMA of the study area and distribution of the measurement point of snow cover

圖11 研究區MODIS積雪識別產品（MOD10A2）Fig．11 MODIS snow products of MOD10A2of the study

同時，本文方法具有很好的擴展性，可推廣應用到我國風云三號中分辨率光譜成像儀（FY－3／MERSI）數據，有利于提升中低分辨率影像積雪覆蓋精度，為氣象業務工作中的融雪徑流模擬、雪災監測等提供數據和技術支持。在上述應用中，MODIS和FY－3／MERSI數據中積雪占混合像元比例多少會直接影響到積雪、植被和土壤等端元光譜的提取，進而影響到積雪分量求解結果。因此，當使用上述中分辨率遙感影像數據時，該方法適用于新疆北疆地區冬季積雪累積期的混合像元分解，這是由于該區域覆蓋了積雪、植被、裸土和戈壁等地貌，地貌類型相對簡單。隨著遙感影像空間分辨率的提高，如采用ETM＋等影像時，積雪占混合像元的比例對于積雪等端元光譜的提取精度影響會隨之降低，該方法適用性也會有所提高。

［1］鄧書斌，陳秋錦．基于MTMF的混合像元分解方法研究［A］．第十七屆中國遙感大會摘要集［C］．南京：中國遙感應用協會，2010：69－73．

［2］潘明忠，元洪興，肖功海，等．基于高精度端元的混合像元線性分解模型研究［J］．紅外與毫米波學報，2010，（5）：357－361．

［3］閻國倩，趙云升，寧艷玲，等．混合像元光譜特征影響因子分析［J］．光譜學與光譜分析，2009，29（12）：3358－3361．

［4］Chabrillat S，Pinet P C，Ceuleneer G，etal．Ronda peridotite massif：Methodology for its geological mapping and lithological discrimination from airborne hyper spectral data［J］．International Journal of Remote Sensing，2000，21：2363－2388．

［5］李慧，陳健飛，余明．線性光譜混合模型的ASTER影像植被應用分析［J］．地球信息科學，2005，7（1）：103－106．

［6］Small C．Estimation of urban vegetation abundance by spectral mixture analysis［J］．International Journal of Remote Sensing，2001，22（7）：1305－1334．

［7］張熙川，趙英時．應用線性光譜混合模型快速評價土地退化的方法研究［J］．中國科學院研究生院學報，1999，16（2）：169－176．

［8］萬軍，蔡運龍．應用線性光譜分離技術研究喀斯特地區土地覆被變化——以貴州省關嶺縣為例［J］．地理研究，2003，22（4）：440－443．

［9］鄒蒲，王云鵬，王志石，等．基于ETM＋圖像的混合像元線性分解方法在澳門植被信息提取中的應用及效果評價［J］．華南師范大學學報（自然科學版），2007，（7）：131－136．

［10］王瑋，黃曉東，呂志邦，等．基于MODIS和AMER－E資料的青藏高原牧區雪被制圖研究［J］．草業學報，2013，22（4）：227－238．

［11］馮琦勝，張學通，梁天剛．基于MOD10A1和AMSR－E的北疆牧區積雪動態監測研究［J］．草業學報，2009，18（1）：125－133．

［12］裴歡，房世峰，覃志豪，等．基于遙感的新疆北疆積雪蓋度及雪深監測［J］．自然災害學報，2008，17（5）：52－57．

［13］李寶林，張一馳，周成虎．天山開都河流域雪蓋消融曲線研究［J］．資源科學，2004，26（6）：23－29．

［14］陳曉娜，包安明，張紅利，等．基于混合像元分解的MODIS積雪面積信息提取及其精度評價——以天山中段為例［J］．應用氣象學報，2004，15（6）：665－671．

［15］延昊，張國平．混合像元分解法提取積雪蓋度［J］．應用氣象學報，2004，15（6）：655－661．

［16］Painter T H，Dozier J，Robiters D，etal．Retrieval of subpixel snow－covered area and grain size from imaging spectrometer data［J］．Remote Sensing of Environment，2003，25（1）：64－77．

［17］Roberts D A，Gardner M，Church R，etal．Mapping chaparral in the Santa Monica mountains using multiple endmember spectral mixture models［J］．Remote Sensing of Environment，1998，65（3）：267－279．

［18］Roberts D，Halligan K，Dennison P．VIPER Tools User Manual Version 1．5［EB／OL］．（2007－12－11）．http：／／www．vipertools．org．

［19］Dennison P E，Halligan K Q，Roberts D A．A comparison of error metrics and constraints for multiple endmember spectral mixture analysis and spectral angle mapper［J］．Remote Sensing of Environment，2004，93（3）：359－367．

［20］Song C．Spectral mixture analysis for subpixel vegetation fractions in the urban environment：How to incorporate endmember variability［J］．Remote Sensing of Environment，2005，95（2）：248－263．

［21］馬孟莉，朱艷，李文龍，等．基于分層多端元混合像元分解的水稻面積信息提取［J］．農業工程學報，2012，8（2）：154－159．

［22］Frankea J，Robertsb D，Halliganc K，etal．Hierarchical multiple endmember spectral mixture analysis（MESMA）of hyperspectral imagery for urban environments［J］．Remote Sensing of Environment，2009，93（3）：1712－1719．

［23］饒萍．EOS－MODIS像元組分分解中端元的選擇與改進［A］．《測繪通報》測繪科學前沿技術論壇摘要集［M］．北京：測繪出版社，2008：1－10．