一種反演積雪面積的范數最小二乘算法

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

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

0 引言

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

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

1 范數最小二乘算法

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

(1)

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

(2)

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

(3)

(4)

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

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

2 實驗與分析

2.1 數據來源及預處理

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

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

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

圖1 研究區的MODIS和TM影像

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

圖2 參考的積雪覆蓋度

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

圖3 積雪端元光譜庫

圖4 求解范數后的積雪光譜庫

2.2 精度驗證

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

表1 3種模型的評價指標

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

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

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

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

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

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

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

3 結束語

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

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

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

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