基于MODIS時間序列的松嫩平原鹽漬地提取

花 錦 溪，臧 淑 英，那 曉 東

(黑龍江省普通高等學校地理環境遙感監測重點實驗室，哈爾濱師范大學，黑龍江 哈爾濱 150025)

基于MODIS時間序列的松嫩平原鹽漬地提取

花 錦 溪，臧 淑 英*，那 曉 東

(黑龍江省普通高等學校地理環境遙感監測重點實驗室，哈爾濱師范大學，黑龍江 哈爾濱 150025)

松嫩平原是我國內陸鹽漬土三大分布區之一，土壤鹽漬化是該區最主要的環境問題。以多時相中分辨率成像光譜儀(MODIS)的歸一化植被指數(NDVI)時間序列影像為主要數據源，通過Savizky-Golay濾波重構NDVI時序數據，依據研究區7種主要土地覆被類型的時間序列曲線差異性，應用分類回歸樹(Classification And Regression Tree,CART)方法確定像素歸屬類別，得到松嫩平原2013年鹽漬土的分布數據;并基于不同鹽漬化程度土壤的植被物候特征差異性建立CART決策樹區分不同程度鹽漬土。分類結果為：鹽漬地掩膜提取精度達98.13%，Kappa系數為0.83；不同程度鹽漬土識別的精度達到86.08%，Kappa系數為0.78。該研究表明多時相MODIS數據在大尺度鹽漬土信息識別中具有可行性。

松嫩平原；MODIS；NDVI；物候特征；鹽漬土

0 引言

松嫩平原是世界三大蘇打鹽漬土集中分布區之一，土壤鹽漬化導致土壤肥力下降，抑制了農作物正常生長，農作物產量下降，也使得許多土地撂荒，畜牧力下降，嚴重制約經濟和社會的發展，導致生態環境惡化[1]，因此對該區域鹽漬土的監測尤為重要。

目前，國內學者對鹽漬土遙感監測的研究主要運用高分辨率、高光譜以及雷達數據。彭望錄等[2]對TM影像進行K-T變換，利用3個分量(亮度、綠度、濕度)提高對鹽堿土的判讀分析效果。駱玉霞等[3]提取TM影像的光譜特征和紋理特征，比較角度分類器和距離分類器，對遙感信息單要素分類與遙感信息和地理信息相結合的綜合分類進行比較，實現土壤鹽漬化分級。陳實等[4]基于Landsat8 OLI遙感影像，引入土壤亮度指數、土壤濕度指數、土壤鹽蓋度指數和土壤輻射水平指數，結合野外實測樣點構建神經網絡,對新疆石河子農墾區土壤含鹽量的反演精度為66%。馬馳[5]以HSI高光譜影像為數據源，結合實地采樣點，建立主要鹽分參數與高光譜數據的偏最小二乘回歸模型，實現對土壤主要鹽分參數的反演。李彪等[5]基于雷達數據，以后向散射系數值與土壤含鹽量之間的反演關系為依據進行決策樹分類，區分河套灌區不同程度鹽漬土。

然而，土壤鹽漬化是一個土地退化的動態變化過程，它不僅表現在土壤、植被等要素某一時間的靜態特征，更呈現出植被等隨著時間周期性動態變化的特征。以上研究均基于單一時相的高分辨率、高光譜和雷達數據，缺少對植被覆蓋時間變化特征的分析。長時間序列的中分辨率成像光譜儀(MODIS)歸一化植被指數(NDVI)數據，能減少云覆蓋因素的影響，精確反映地表植被覆蓋的時間變化特征，在遙感信息提取中顯現出獨特的優勢。那曉東[6]利用多時相MODIS NDVI數據，依據植被物候特征，運用傅里葉組分相似度的監督分類方法提取三江平原濕地植被的分布信息，分類精度達到79.67%；謝相建等[7]基于多時相MODIS NDVI影像，對比BFAST模型分解出的關鍵物候特征和光譜特征對土地利用覆蓋的分類效果，發現物候特征對地物識別的精度較高。以上學者為多時相遙感分類提供了有效支撐，但將此方法應用到鹽漬土信息提取的研究剛剛開始。本研究利用經過Savizky-Golay濾波平滑后的23個時序MODIS NDVI特征進行CART分類，實現松嫩平原鹽漬地掩膜的提取,并利用土壤中植被物候特征的差異性識別不同程度鹽漬土。

1 研究區與數據處理

研究區松嫩平原位于黑龍江省西南部和吉林省西北部(119°45′～129°36′E,42°50′～49°18′N),行政區劃主要包括黑龍江省的35個縣和吉林省的20個縣，總面積約23.75萬km2。整個松嫩平原鹽漬土區涉及有黑龍江省的14個市縣、吉林省的14個市縣和內蒙古的7個市縣[8]。松嫩平原以大興安嶺東麓丘陵和臺地為界，北部和東部以小興安嶺及長白山地外緣山麓臺地為鄰，南抵達松遼分水嶺，由松花江和嫩江沖積而成。該研究區環境背景下形成的土壤類型比較復雜，主要有黑鈣土、淡黑鈣土、栗鈣土、草甸土、鹽土、堿土、風沙土、沼澤土等。本區地處半干旱-半濕潤的交錯地帶，屬于典型的溫帶大陸性氣候，年均降水量為400～500 mm，集中在夏季，并且由東南向西北遞減，年均蒸發量為1 250～1 650 mm，遠大于降水量,使得氣候干旱。海拔150～200 m，地勢低洼，排水不暢，為土壤積鹽提供了地形條件，是典型的蘇打鹽漬土代表區[9]。

本研究采用的數據源來自美國國家宇航局(NASA)提供的2013年全年MODIS NDVI 16 d合成數據(MOD13Q1)，空間分辨率為250 m。根據松嫩平原地理位置，選取軌道號為H26V04和H27V04的MODIS遙感影像，共46景。MOD13Q1產品屬于MODIS的4級遙感影像數據，已經過輻射校正和幾何校正，在此基礎上還需要進行拼接、投影、掩膜等預處理。利用MODIS數據處理工具對影像進行批量拼接，并投影為阿爾伯斯圓錐等面積投影。對23幅MODIS影像進行疊加，對應圖層序列號為1～23,并利用松嫩平原矢量邊界在ArcGIS中進行掩膜處理，提取出松嫩平原的范圍。

MODIS NDVI 16 d合成數據在合成期間難免會受到傳感器角度、氣溶膠變化等因素的影響而產生噪聲，使得NDVI曲線的季節變化趨勢及物候信息受到干擾而無法精確提取地表覆蓋信息[10,11]。本研究對MODIS數據在TIMESAT3.0中進行Savizky-Golay濾波，此方法采用最小二乘卷積平滑算法進行曲線擬合[12]，通過反復迭代使平滑后的曲線逼近原始曲線的上包絡線[13],從而去除損壞數據、降低噪聲，擬合NDVI曲線的季相信息，更好地再現各地類植被的季相生長規律。實際濾波中需設置窗口寬度和平滑次數，本研究根據MODIS NDVI數據集設置最佳濾波窗口寬度為4，平滑次數為2。

根據2015年4月初實地調查對2015年松嫩平原TM影像進行目視解譯，結合土地利用分類標準以及鹽漬土分類標準，選取研究區7種主要土地覆蓋類型(耕地、林地、草地、沼澤、水體、城鎮、鹽漬土)以及輕度、中度、重度鹽漬土的樣本。訓練樣本的選取越具有代表性分類結果越準確，不同等級鹽漬土在近紅外、紅光、綠光波段假彩色組合的TM遙感影像上特征明顯。重度鹽漬土在TM影像上為亮白色，常分布在湖泡邊緣，隨著季相變化水位上升和下降，可溶性鹽分在表層積累，形成鹽堿斑，形狀不規則，分布區內植被很少；中度鹽漬土在TM影像上為灰白色、淺黃色，形態不規則，呈斑塊狀，一般位于小湖泡周圍和地勢低洼的地方，該類型鹽漬地中有少量植被生長；輕度鹽漬土在TM影像上為粉色、淡藍色，分布在農田大沖溝邊緣，幾何形態規則，分布區內有較多的星點狀植被和塊狀農田[14]。將所選樣本(表1)隨機選取約2/3用于訓練，1/3用于驗證。

表1 樣本數量Table 1 Number of samples

2 研究方法

2.1 鹽漬地遙感信息提取方法

2.1.1 MODIS NDVI特征選取及可分性驗證 統計7種主要地類樣本的NDVI時間序列特征(圖1)，可以看出：林地植被NDVI值相對較高且波動平穩，生長季開始較早且持續時間長；耕地的時序曲線振幅大，農田管理、施肥、灌溉為作物創造了優越的生長條件，呈現出較高的峰值；沼澤地帶水分條件充足，植被生長季NDVI值高于草地，達到峰值的時期較早；水域幾乎沒有植被覆蓋，其NDVI值全年處于較低水平；城鎮和鹽漬土的時序曲線振幅較小且數值較低，但鹽漬土時序曲線具有明顯的季相變化特征，而城鎮的時序曲線并不體現季相變化信息。鹽漬土特有的NDVI時間序列特點表明基于MODIS NDVI時序數據集提取鹽漬土可行，因此,以23景MODIS NDVI影像為特征提取松嫩平原鹽漬地。

圖1 7種主要土地類型NDVI時間序列特征Fig.1 NDVI time series features of seven major land types

2.1.2 建模方法 以研究區各地類6 000個訓練樣本點的23個MODIS NDVI特征值作為預測變量，目標變量為7種土地覆被類型，采用salford-system公司的SPM軟件進行機器學習構建分類回歸樹。這種計算機自動分類方法能夠克服人為主觀判斷，同時可以快速處理高維數據。分類回歸樹(CART)算法采用經濟學中的基尼系數(Gini Index)作為選擇最佳測試變量和分割閾值的準則構建決策樹；并且用交叉驗證的方法對決策樹結構進行修剪以防止“過度擬合”的現象，即對某個樹枝的所有葉節點增加一個懲罰因子，如果增加該樹枝錯誤率降低則保留，否則予以剪除，最終得到一棵兼顧復雜度和錯誤率的最優二叉樹[14,15]。最初生成的決策樹節點為101個，可調誤差率為0.032,隨著決策樹節點數目的增加可調誤差率逐漸降低。當節點數目為20時可調誤差率為0.094，并且基本保持該值不變。因此，選取前20個根節點的CART決策樹結構，根據該分類規則得到松嫩平原鹽漬土掩膜。

2.2 鹽漬化程度的分類方法

2.2.1 物候特征選取及可分性驗證 松嫩平原區鹽漬土面積所占比例最大的是草地[16]，土壤鹽漬化使得植被生長狀況發生改變，植被指數時序曲線的差異有助于推測鹽漬化程度。如圖2，隨著土壤含鹽量的增加，植被生長季NDVI值逐漸減小。土壤含鹽量低時植被長勢好，NDVI值在160～200 d之間迅速增長，并且在較長的時間內保持較高的值。在含鹽量較高的土壤中幾乎沒有植被發育，NDVI值處于較低水平。從整個季相看，不同含鹽量土壤中的植被在萌芽、成熟、衰老時期NDVI時序曲線存在一定的相似性，但在植被生長旺盛時期曲線差異較大、可分性較高。為了進一步驗證不同程度鹽漬土中植被物候特征的差異性，對各級鹽漬土參考像素進行統計，將NDVI最大值這一物候特征用盒須圖表示。如圖3所示，重度鹽漬化土壤中的植被類內差異較大，但它與輕度、中度鹽漬化土壤中的植被類間差異更大，表明其可分性較好；輕度、中度鹽漬化土壤中植被類內差異較小、類間差異明顯，呈現出較好的可分性。表明植被物候特征的差異性指示不同含鹽量土壤的可行性，因此，以11個植被物候參數為特征識別不同程度的鹽漬土。植被物候特征分別是植被生長季開始、結束、長度、基值、中值、峰值、時序曲線振幅、左導數、右導數、小積分、大積分。各物候參數代表的含義如表2所示。

圖2 各級鹽漬化土壤NDVI曲線Fig.2 NDVI curves of different salinization soils

圖3 各級鹽漬土NDVI最大值的盒須圖Fig.3 Box plot for NDVI maximum value of different salinization soils

表2 物候參數代表含義Table 2 The meaning of each phenological parameter

2.2.2 建模方法 以各級鹽漬土310個訓練樣本點中植被的11物候特征作為預測變量，不同等級鹽漬土作為目標變量，在SPM中構建分類回歸樹。當節點數目為3時可調誤差率為0.194，并且隨著節點數的增加可調誤差率基本不變。因此，選取節點數為3時的規則在ENVI5.2中對遙感影像進行分類,實現不同程度鹽漬土的分類。

3 分類結果及精度驗證

3.1 鹽漬地信息提取結果與精度評價

基于MODIS NDVI特征的CART分類得到松嫩平原鹽漬地掩膜。第一步分類過程中最重要的變量是第14個時期的NDVI值，這個時期植被生長最旺盛，達到了NDVI的最大值。整體看，植被生長旺盛期各地類時序曲線差異較大,可分性好，NDVI值對鹽漬土提取的重要性較高。

對于低分辨率影像進行驗證的方法是實地驗證和基于高空間分辨率影像進行驗證，考慮到本研究區范圍較大，實地驗證難以實現，因此由野外實測結合TM影像目視解譯結果樣本的1/3進行驗證。利用隨機選取的2 800個非鹽漬土像素點和200個鹽漬土像素點驗證鹽漬地提取結果，混淆矩陣如表3所示，提取精度達到98.13%，Kappa系數為0.83，表明MODIS NDVI時間序列的差異性對鹽漬土同其他地類的區分效果較好，鹽漬土提取精度達到進一步分類的要求，能夠進行不同程度鹽漬土的提取。

表3 鹽漬土掩膜提取精度Table 3 Extraction accuracy of saline soil mask

3.2 鹽漬化程度遙感反演結果與精度評價

基于不同程度鹽漬土中植被物候特征的CART分類最終得到松嫩平原鹽漬化程度結果(圖4)，為了更清晰地表達不同鹽漬化程度空間分布格局，對鹽漬土集中分布區域大安、乾安等縣市放大顯示。

第二步分類過程中各變量的重要性如圖5，在不同程度鹽漬土識別中發揮最重要作用的變量是植被生長季峰值，這與提取鹽漬土掩膜的結果相一致。植被物候特征中時序曲線大積分、振幅也對鹽漬化程度遙感信息提取起到了重要作用。它們均表示植被在生長季期間的長勢強度，反映了生長季期間植被光合作用值。

圖4 不同程度鹽漬土分類結果Fig.4 Classification results of different degrees of saline soil

圖5 植被物候特征CART分類變量重要性得分Fig.5 Variable importance scores for CART classification of phenological parameters

基于植被物候特征識別不同鹽漬化程度土壤的混淆矩陣如表4所示，分類總體精度達到86.08%，Kappa系數為0.78，可以看出分布完整的輕度鹽漬土和重度鹽漬土分類精度較高，中度鹽漬土分布零散、分類精度較低。由于輕度鹽漬土與中度鹽漬土常交錯分布，容易造成錯分混分。Moncef Bouaziz等[17]應用線性光譜分解技術提高了電導率與18個MODIS指數的相關性，為降低鹽漬地的錯分、混分率提供了參考。一些學者嘗試運用能詳細表征地表信息的高分辨率影像提取鹽漬地信息，但提取效果并不理想，如：李彪等[6]基于雷達影像進行決策樹分類提取輕、中、重度鹽漬地的精度分別為79.48%、78.21%、86.27%，李晉等[18]采用定量定性相結合，基于TM影像進行決策樹分類提取不同程度鹽漬地的精度達到82.49%、77.05%、91.56%。與以往研究相比，本研究所采用的MODIS空間分辨率較低但時間分辨率較高，從中得到的植被物候特征對不同程度鹽漬地識別效果較好，充分顯示了植被物候特征在鹽漬地信息提取中的優勢。

表4 不同程度鹽漬土提取精度Table 4 Extraction accuracy of different degrees of saline soil

4 結論與討論

本研究首次基于多時相MODIS數據,應用CART方法實現鹽漬地遙感信息提取。首先將鹽漬土同其他主要土地覆被類型進行區分，鹽漬土掩膜的提取精度達到98%，Kappa系數為0.83。其次實現不同程度鹽漬土的識別，分類精度達到86.08%，Kappa系數為0.78。多時相數據反映了豐富的地表植被變化信息，很好地指示了不同程度含鹽量土壤。同時，結合CART分類方法克服了人為主觀因素的影響，使得分類結果更加客觀真實。研究發現植被物候參數對成片分布的輕度、重度鹽漬土的分類效果較好，對散布的中度鹽漬土分類效果不好，這主要是受到MODIS數據空間分辨率的限制。在兩步分類過程中，表征植被光合作用程度的變量發揮了重要作用，今后應考慮篩選出重要的變量,運用低維數據達到識別鹽漬土信息的目的。

本研究還應該注意的是，MODIS數據空間分辨率較低且存在混合像元問題，應嘗試運用混合像元分解等方法提高不同程度鹽漬地提取的精度。本研究基于MODIS NDVI時序數據和物候特征實現鹽漬土信息的提取，通過植被信息指示土壤含鹽量，沒有考慮地形、氣候等因素對土壤鹽分信息的影響，今后還應結合多種數據源挖掘鹽漬土信息。

[1] 鄧偉，裘善文，梁正偉.中國大安堿地生態試驗站區域生態環境背景[M].北京：科學出版社，2006.46.

[2] 彭望琭，李天杰.TM數據的Kauth-Thomas變換在鹽堿土分析中的作用——以陽高盆地為例[J].環境遙感，1989，4(3)：183-190.

[3] 駱玉霞，陳煥偉.GIS支持下的TM圖像土壤鹽漬化分級[J].遙感信息，2001(4)：12-15.

[4] 陳實，高超，徐斌,等.新疆石河子農區土壤含鹽量定量反演及其空間格局分析[J].地理研究，2014，33(11)：2135-2144.

[5] 馬馳.基于HJ1A-HIS反演松嫩平原土壤鹽分含量[J].干旱區研究，2014，21(2)：226-230.

[6] 李彪，王耀強.土壤鹽漬化雷達反演模擬研究[J].干旱區資源與環境，2015，29(8)：180-184.

[7] 那曉東，張樹清，李曉峰,等.MODIS NDVI時間序列在三江平原濕地植被信息提取中的應用[J].濕地科學，2007，5(3)：227-236.

[8] 謝相建，薛朝輝，王冬辰，等.顧及物候特征的喀斯特斷陷盆地土地覆蓋遙感分類[J].遙感學報，2015，19(4)：627-638.

[9] 王晶，肖延華，朱平，等.松嫩平原鹽漬土的發展演化與影響因素[J].吉林農業科學，1994，2：66-71.

[10] 吉林省土壤肥料總站.吉林省土壤[M].北京：中國農業出版社，1998.195-211.

[11] 侯東，潘耀忠，張錦水,等.農區MODIS植被指數時間序列數據重建[J].農業工程學報，2010，26(S1)：206-212.

[12] 于信芳，莊大方.基于MODIS NDVI數據的東北森林物候期監測[J].資源科學，2006，28(4)：112-117.

[13] 楊恒，沈潤平，吳立葉，等.基于S-G濾波的江西省植被覆蓋度時空變化遙感分析[J].科學技術與工程，2014，22(14)：101-106.

[14] SAVITZKY A,GOLAY M J E.Smoothing and differentiation of data by simplified least squares procedures[J].Analytical Chemistry，1964，36：1627-1639.

[15] 李俊杰，何隆華，戴錦芳,等..基于遙感影像紋理信息的湖泊圍網養殖區提取[J].湖泊科學，2006，18(4)：337-342.

[16] 那曉東，張樹清,孔博，等.基于決策樹方法的淡水沼澤濕地信息提取——以三江平原東北部為例[J].遙感技術與應用，2008，23(4)：365-372.

[17] 王晶，肖延華，朱平,等.松嫩平原鹽漬土的發展演化與影響因素[J].吉林農業學，1995(2)：66-71.

[18] 李晉,趙庚星,常春艷,等.基于HSI高光譜和TM圖像的土地鹽漬化信息提取方法[J].光譜學與光譜分，2014(2)：520-525.

Saline Land Extraction Using MODIS NDVI Time Series Data in Songnen Plain

HUA Jin-xi,ZANG Shu-ying,NA Xiao-dong

(KeyLaboratoryofRemoteSensingMonitoringofGeographicEnvironment，CollegeofHeilongjiangProvince,HarbinNormalUniversity,Harbin150025,China)

Songnen Plain is one of China′s inland saline soil distribution area and its salinization degree has aggravated,soil salinization is the main environmental problem in this area.It is necessary to monitor the saline soil information of Songnen Plain.Firstly,we applied Savitzky-Golay filter to multi-temporal MODIS NDVI data,the 16 day composite product with 250 meter resolution,to reconstruct NDVI time series curves.Based on the difference of NDVI curves of different land cover types,CART method was used to extract saline soil of Songnen Plain.Secondly,we also established a CART decision tree to identify different degrees of saline soil by using 11 phonological parameters.In the process of classification,the maximum value of NDVI plays the most important role to identify saline soil information which represents the strength of vegetation photosynthesis.It also explains that the significant difference of NDVI curves is valuable to distinguish saline soil information when NDVI at higher levels. Results showed that the accuracy for extract saline soil of Songnen Plain is 98.13%.At the same time,Kappa coefficient is 0.83.The accuracy for identify different degrees of saline soil is 86.08%,and Kappa coefficient is 0.78.This study indicates that the applying of CART method to MODIS time series data is feasible to acquire salt information in a large scale.

Songnen Plain;MODIS;NDVI;phonological character;saline soil

2015-11-09；

2016-01-31

國家自然科學基金面上項目(41571199);黑龍江省自然科學基金項目(D201409);國家自然科學基金青年項目(41001243);黑龍江省普通高校青年骨干學術項目(1253G034);哈爾濱師范大學碩士研究生創新基金重點項目(HSDSSCX2015-11)

花錦溪(1991-)，女，碩士研究生,主要從事生態環境遙感應用研究。*通訊作者E-mail：zsy6311@163.com

10.3969/j.issn.1672-0504.2016.02.013

TP75

A

1672-0504(2016)02-0067-05