段金亮，張瑞，2，李奎，龐家泰

(1.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，成都 611756；2.西南交通大學(xué) 高速鐵路運營安全空間信息技術(shù)國家地方聯(lián)合工程實驗室，成都 611756)

0 引言

積雪是地球表面覆蓋的重要組成部分，是地球表面最為活躍的自然元素之一，也是冰凍圈中地理分布最廣泛、季節(jié)與年際變化最顯著的部分[1]。積雪的特征包括積雪面積、積雪分布、雪深、雪粒徑等，其中積雪面積是全球能量平衡、氣候、水文、生態(tài)模型以及積雪定量遙感中的重要輸入?yún)?shù)之一[2]。當(dāng)前，積雪面積的反演主要集中在光學(xué)遙感影像上，方法主要有閾值法[3]、監(jiān)督分類法[4]、光譜混合分析法等。閾值法具有簡單快速的特點，但是不同影像的閾值不固定，每景影像都要重新選擇閾值，通用性差；監(jiān)督分類法雖然有一定的魯棒性，但是需要大量的樣本作為支撐條件，這在生產(chǎn)實踐中很難滿足。同時，閾值法和監(jiān)督分類法都是從像元尺度來獲取積雪面積，由于積雪覆蓋范圍具有較強(qiáng)的時空變異性特征，以及受影像空間分辨率的限制，積雪通常與土壤、巖石、植被等混合成混合像元[5-6]。如果單獨從傳統(tǒng)的像元尺度分類算法來考慮，很難解決積雪與其他地物混合的問題，往往難以獲得更有效、精度更高的積雪面積產(chǎn)品，所以從亞像元尺度考慮的光譜混合分析法在反演積雪面積上很有優(yōu)勢[7-8]。全約束最小二乘法(fully constrained least squares,FCLS)是遙感影像處理中常用的光譜混合分析模型[9]，其端元光譜庫被假定在整景遙感影像中是不變的，然后用其來解混整景遙感影像。但是在實際情況中，端元光譜存在異物同譜和同譜異物的現(xiàn)象，同時端元是沿著影像進(jìn)行變化的，這種變化導(dǎo)致了端元變化。這種端元變化減低了遙感影像的解混精度，端元變化已經(jīng)被國內(nèi)外學(xué)者證實是豐度估計的主要誤差源[10-11]。

針對遙感影像中端元變化的問題，國內(nèi)外學(xué)者提出了多種光譜混合模型[9-14]，其中比較經(jīng)典的是多端元光譜混合分析法(multiple endmember spectral mixture analysis，MESMA)。雖然MESMA算法能夠很好地處理端元變化問題，但是MESMA算法的運算效率存在一定的局限。同時，由于積雪的時空變化較大，造成了端元變化也隨之較大，同時還與其他地物進(jìn)行混合。基于以上考慮，本文對常規(guī)的線性混合模型全約束最小二乘法進(jìn)行改進(jìn)，提出范數(shù)最小二乘算法(norm least squares,NLS)，通過引入范數(shù)來降低積雪端元在可見光波段的相對差異，以減弱端元變化的影響，從而提高反演積雪面積的精度。

1 范數(shù)最小二乘算法

在遙感影像的特征提取和分類中，常用的光譜分析模型是線性混合模型(linear mixed model)。其通過提取研究區(qū)典型地物的光譜作為端元光譜庫，然后應(yīng)用于線性模型上獲取各個地類的豐度值。常用的線性混合模型是全約束最小二乘法[9]，表達(dá)如式(1)所示。

(1)

式中：y∈1×L為混合光譜向量；e∈1×L為端元向量；a為豐度值；ε∈1×L為噪聲向量；M是端元數(shù)目；L為波段的數(shù)目。其中，豐度值滿足非負(fù)、和為1 的限制，如式(2)所示。

(2)

因為光學(xué)影像在成像時受到地物物理特性的變化、光照條件和大氣環(huán)境的影響，存在同物異譜和同譜異物的情況，從而產(chǎn)生端元變化，導(dǎo)致光譜混合分析算法的精度降低。同時，積雪覆蓋的MODIS影像的端元變化較大，選擇合理的端元存在困難。此外，若應(yīng)用于解混的端元光譜庫的相關(guān)性過高，則遙感影像的解混精度也會降低，從而降低反演積雪面積的精度。此時，常規(guī)的FCLS算法就有一定局限性，特別是針對積雪這類光譜變化性較大的地物，常規(guī)的FCLS算法很難解決積雪的端元變化，導(dǎo)致其算法的反演積雪面積的精度不高。故此處利用范數(shù)最小二乘法來減弱端元變化，改寫式(1)為式(3)至式(4)所示。

(3)

(4)

式中：x∈1×L為被用于‖·‖p處理的值。

為減弱端元變化對解混結(jié)果的影響，這里每個端元光譜除以它們的2范數(shù)，這樣壓縮了數(shù)據(jù)空間，縮小了端元間相對差異，減弱了端元變化的絕對量，從而減弱了端元的變異性，提高了解混精度，同時保留了它們在高維特征空間的相對位置。

2 實驗與分析

2.1 數(shù)據(jù)來源及預(yù)處理

MODIS影像的光譜覆蓋度廣(36個光譜波段)，空間分辨率適中，且重訪周期短，獲取數(shù)據(jù)便捷，適合用于大范圍的積雪面積的提取。為此本文對MODIS數(shù)據(jù)進(jìn)行反演實驗，驗證NLS算法反演積雪面積的可行性，同時為驗證其反演的積雪面積精度，選擇同一時間段的TM影像進(jìn)行精度評價。TM影像編號為：LT51450362009307KHC00，其成像時間為：2009年11月3日。該影像云量小、地物類型單一，有大范圍的積雪覆蓋，適合積雪面積反演算法的研究。

首先，將HDF格式的MODIS影像轉(zhuǎn)換成ENVI的標(biāo)準(zhǔn)格式，重新投影MOD09GA產(chǎn)品，使其與TM影像的坐標(biāo)系統(tǒng)保持一致；其次，按照波長大小對波段進(jìn)行升序排列重新合并。由于MODIS影像存在大量的壞像元，同時第5波段存在明顯的條帶，為了減弱它們的影響，采用空間與光譜結(jié)合的去噪算法[15]，其參數(shù)設(shè)置如下：[0.1，0.2，0.2，40]，分別對應(yīng)λ，u，v，MaxIter[16]。然后，對TM影像進(jìn)行預(yù)處理，去除大氣誤差獲取地物反射率，其處理過程參考文獻(xiàn)[17]。

因為無法獲取研究區(qū)真正的積雪面積，此處采用空間分辨率高的TM影像制作研究區(qū)的積雪覆蓋度的參考值，通過幾何配準(zhǔn)讓MODIS影像與TM影像處于同一參考坐標(biāo)系下，然后用TM影像去裁剪MODIS影像得到研究區(qū)域如圖1所示(波長對應(yīng)于TM影像的7、4、3 波長)。

圖1 研究區(qū)的MODIS和TM影像

積雪面積驗證數(shù)據(jù)需要從分類后的TM影像中獲得。首先，對TM影像進(jìn)行輻射定標(biāo)和大氣校正后，結(jié)合天地圖和TM影像目視判讀研究區(qū)地物主要包括積雪、裸巖、水體、土壤和陰影5個大類；其次，通過對TM影像進(jìn)行目視判讀，同時參考文獻(xiàn)[18]選擇分類樣本和驗證樣本；再次，利用ENVI軟件中的支持向量機(jī)算法(support vector machine algorithm，SVM)進(jìn)行分類，其參數(shù)設(shè)置參考文獻(xiàn)[19]；然后，使用驗證樣本對分類結(jié)果進(jìn)行精度驗證，其分類結(jié)果的總體精度達(dá)到88.5%，Kappa系數(shù)達(dá)到0.86。獲得TM影像的分類結(jié)果后，通過目視判讀改正分類錯誤的積雪像元，判讀時參考天地圖的高分辨率影像逐一改正。最后,把經(jīng)過改正后的積雪類別保存出來。由于TM影像的空間分辨率(30 m)與MODIS影像的空間分辨率(500 m)不一樣，為滿足后期反演積雪面積算法的精度驗證實驗需要，對分類后的積雪類別影像進(jìn)行像元聚合，使積雪類別影像與MODIS影像的空間分辨率保持一致，從而獲得研究區(qū)的積雪覆蓋度的參考值。研究區(qū)的積雪覆蓋度參考影像如圖2所示。

圖2 參考的積雪覆蓋度

光譜混合分析算法在反演積雪面積前，需要提前構(gòu)建端元光譜庫，這里使用空間分塊自動端元提取算法[20]從研究區(qū)的MODIS影像上提取積雪、裸巖、水體、土壤和陰影5類端元，然后剔除異常的端元，最后每個類別得到15個端元。由于本文主要反演積雪面積，所以展示積雪光譜曲線如圖3所示，圖中積雪的光譜曲線變化較為明顯，主要是可見光譜波段，最大變化達(dá)到0.28，近紅外波段變化相對較小。通過對光譜進(jìn)行范數(shù)計算后，如圖4所示，引入范數(shù)計算，縮小端元變化的絕對量，減少端元間的相對差異，從而在一定程度上減弱積雪的光譜變化，提高反演的精度。

圖3 積雪端元光譜庫

圖4 求解范數(shù)后的積雪光譜庫

2.2 精度驗證

對FCLS、NLS和MESMA算法反演積雪面積的時間效率進(jìn)行對比分析，在電腦處理器為Intel(R)Core(TM)i5-4210M CPU@2.60 GHz、運行內(nèi)存為8.00 GB和64位windows7的操作系統(tǒng)下，對7 000行8 000列大小的數(shù)據(jù)進(jìn)行反演計算，最后得到3種模型的運行時間，如表1所示。

表1 3種模型的評價指標(biāo)

為驗證FCLS、NLS和MESMA算法反演的積雪面積的精度，使用TM分類影像作為參考值，同時使用絕對差異值、均方根誤差(RMSE)、相關(guān)系數(shù)(R)、總的積雪面積(TSCA)指標(biāo)來進(jìn)行精度評價，其FCLS、NLS和MESMA算法的精度評價指標(biāo)對比結(jié)果如表1所示。針對FCLS、NLS和MESMA算法反演出的結(jié)果，對積雪覆蓋度進(jìn)行求和與歸一化，其中NLS、MESMA和FCLS反演的積雪覆蓋度如圖5所示，其與參考值的絕對差異如圖6所示，其絕對差異值分段統(tǒng)計結(jié)果如表2所示。

圖5 NLS、MESMA和FCLS算法反演的積雪覆蓋度

圖6 NLS、MESMA和FCLS算法的絕對差異圖

通過將圖5(a)和圖5(b)與參考圖對比，發(fā)現(xiàn)圖5(a)十分接近參考影像，同時圖6(a)和圖6(b)的對比，進(jìn)一步說明圖5(a)的結(jié)果十分接近參考值，再結(jié)合表1的3種評價指標(biāo)來看，NLS的評價指標(biāo)參數(shù)都要優(yōu)于FCLS。而在表2中，F(xiàn)CLS的絕對差異值分布在0.0～0.1區(qū)間的只占有44.41%，且在0.2～0.4區(qū)間都超過了10%，而NLS的絕對差異值主要分布在0.0～0.1區(qū)間，其次分布在0.1～0.2區(qū)間，且0.0～0.1區(qū)間的值高達(dá)63.52%。綜上，可以得出NLS反演積雪面積明顯優(yōu)于FCLS，說明引入范數(shù)減弱了端元的變化，從而提高了算法的精度。

表2 3種算法的絕對差異值

通過分析表1和表2中MESMA和NLS的評價參數(shù)，可以發(fā)現(xiàn)其RMSE、R、TSCA的值相對一致，同時其與真值的絕對差異值在各個區(qū)間基本相似，可以證明NLS的反演積雪面積效果跟MESMA反演積雪面積的效果非常吻合。通過分析評價指標(biāo)的值，可以得出，85.71%的絕對差異值都集中在0.0～0.3范圍，且在0.0～0.1的范圍高達(dá)62.82%，最重要的是TSCA更是高達(dá)0.95以上，從而可以說明MESMA和NLS反演的積雪面積與真值十分接近，且可信度高，可以用MESMA和NLS來反演MODIS影像的積雪面積。雖然MESMA和NLS都可以用反演積雪面積，但是通過表2的時間指標(biāo)可以看出來，NLS的時間效率是MESMA的幾十倍，所以在這點上NLS要優(yōu)于MESMA。

綜上，NLS反演積雪面積的效果要遠(yuǎn)優(yōu)于FCLS，且與MESMA的結(jié)果保持一致，都能夠很好地來反演積雪面積，且與真值十分接近，說明減弱端元變化能夠顯著提高反演積雪面積的精度。同時驗證了MESMA算法通過直接將同類物質(zhì)的所有光譜參與光譜分解來減弱端元變化，但是增加算法的時間復(fù)雜度，消耗大量的時間。而NLS通過引入范數(shù)的方法同樣成功消除了端元變化的影響，且反演效果很理想，時間效率比MESMA高幾十倍。

3 結(jié)束語

本文基于MODIS影像反演積雪面積提出了一種范數(shù)最小二乘法，來解決端元變化和積雪光譜的異質(zhì)性問題，從而提高反演積雪面積的精度。實驗結(jié)果表明，NLS在反演MODIS影像上的積雪面積方面的可行性強(qiáng)且應(yīng)用前景廣闊。本文得出以下結(jié)論。

1)NLS反演積雪面積的精度遠(yuǎn)遠(yuǎn)高于FCLS，引入范數(shù)減低端元變化和光譜異質(zhì)性的影響，提高了反演精度。

2)NLS反演的積雪面積和MESMA保持高度一致，說明其解決端元變化的效果跟MESMA一樣，且時間效率比MESMA高，其更適用于反演MODIS影像的積雪面積。

NLS算法對于高效反演MODIS影像的積雪覆蓋度有較好的效果，但仍有改進(jìn)空間。此算法沒有考慮空間領(lǐng)域信息、波段噪聲權(quán)重對反演積雪面積的影響，同時其運行效率還有待進(jìn)一步提高[21-22]。