NDVI相似性分區(qū)下天山地區(qū)草地總產(chǎn)草量遙感估算

劉 艷，聶 磊，楊 耘

(1.中國氣象局烏魯木齊沙漠氣象研究所,新疆 烏魯木齊 830002; 2.中亞大氣科學(xué)研究中心,新疆 烏魯木齊 830002; 3.武漢大學(xué)測繪遙感信息工程國家重點實驗室,湖北 武漢 430079; 4.長安大學(xué)地質(zhì)工程與測繪學(xué)院,陜西 西安 710054)

新疆現(xiàn)有牧草地總面積5 116.07萬hm2，占全國牧草地總面積的19.52%[1]。北疆地區(qū)天然草地有效面積約2 377.52萬hm2，占全疆46.47%，是新疆草原畜牧業(yè)的重要基地[2]。建立草原產(chǎn)草量監(jiān)測模型以準確、及時地估算草原產(chǎn)草量對確定載畜量、合理安排畜牧業(yè)生產(chǎn)和維護草原生態(tài)平衡等具有重要意義。新疆牧草產(chǎn)量估算研究主要為90年代初開展的NOAA/AVHRR衛(wèi)星數(shù)據(jù)和天然草地牧草產(chǎn)量關(guān)系模型研究[2-7]；地面測產(chǎn)數(shù)據(jù)和同步MODIS遙感數(shù)據(jù)集合的天山北坡烏魯木齊草甸、草原和荒漠草原草地牧草估產(chǎn)模型[8]；遙感分類結(jié)果和產(chǎn)量估測模型集合的阜康市1990-2008年近20年草地總面積和總產(chǎn)量趨勢變化[9]；1991-1996年新疆天山北坡阜康市內(nèi)不同草地類型實測草地可食產(chǎn)量、環(huán)境與遙感資料等綜合的多重相關(guān)分析和遙感估產(chǎn)技術(shù)研究[10]。這些研究主要集中在1990-2000年；研究區(qū)域有限且集中分布在天山北坡的烏魯木齊(南郊)、阜康和阿勒泰；模型采樣點少且時空分布不均，致使模型代表性不足；模型主要針對一個觀測年進行，時間序列性不強。僅僅使用一個觀測年內(nèi)高寒草甸、山地草甸、山地草原等7種不同牧草類型逐月產(chǎn)量觀測數(shù)據(jù)，進行產(chǎn)量變化特征分析和實現(xiàn)單一牧草類型產(chǎn)量計算。新疆草地分布錯綜復(fù)雜，上述模型代表性和時效性都非常有限。NOAA/NDVI植被指數(shù)產(chǎn)品在空間和光譜分辨率上劣于MODIS/NDVI植被指數(shù)產(chǎn)品，相對獲取不便，特別是近年來的植被指數(shù)產(chǎn)品[11]。MODIS/NDVI植被指數(shù)產(chǎn)品空間分辨率有250 m、500 m和1 km，共3種，有單日、16日和月合成3類數(shù)據(jù)，產(chǎn)品時空分辨率能很好地反映草地植被時空變化和估算草地生物量[12-18]。因此，根據(jù)前人已有研究成果不足，考慮到如何有效合理地利用有限樣點數(shù)據(jù)，本研究選取MOD13Q1歸一化植被指數(shù)(NDVI)16 d合成產(chǎn)品為試驗數(shù)據(jù)，從遙感監(jiān)測植被指數(shù)分布特征入手，基于計算機圖像直方圖識別手段定量評價遙感監(jiān)測植被指數(shù)的區(qū)域相似性，以縣(市)為研究單元，計算多年植被指數(shù)均值直方圖巴氏距離，定量評價各單元及其臨近單元NDVI頻數(shù)直方圖相似性以得到有效遙感分區(qū)，以此分區(qū)為基礎(chǔ)，建立NDVI-草地產(chǎn)草量統(tǒng)計模型，將其推廣至臨近具有相似性NDVI的區(qū)域，從而獲取整個天山山區(qū)牧草(鮮重)總產(chǎn)草量遙感估算模型。

1 數(shù)據(jù)來源與研究方法

1.1 研究區(qū)概況

按草地植被型組歸納新疆草地有草原、荒漠、草甸、沼澤4個類組。草原草地除阿爾泰山南麓沖積扇上部和天山北麓山前傾斜平原存在平原荒漠草原片段外都發(fā)育在山地?；哪莸卦谛陆植紡V、面積大，包括溫性草原化荒漠、溫性荒漠和高寒荒漠類。草甸草地分布面積僅次于荒漠和草原草地，包括低地草甸、山地草甸和高寒草甸3種類型。沼澤草地分布面積不大，分別占新疆草地面積和可利用土地面積的0.47%和0.50%，主要分布在南北疆各大河流下游或河口、河漫灘、湖泊周圍、山前潛水溢出帶山間低位盆地的湖濱和山地高位盆地等，植被是由濕中生多年生根莖禾草和莎草類組成。由表1和圖1可見，本研究區(qū)域草地在類型上具有一定的典型性和完整性，在空間分布上具有一定的連續(xù)性，可作為草地總產(chǎn)草量遙感估算的典型試驗區(qū)域。

表1 研究區(qū)4類典型草地統(tǒng)計面積及比例Table 1 Area and proportion of four kinds of typical grasslands in the study area

圖1 研究區(qū)4類典型草地空間分布圖Fig. 1 Spatial distribution of four kinds of typical grasslands in the study area

1.2 數(shù)據(jù)來源

數(shù)據(jù)來源包含地面監(jiān)測和衛(wèi)星遙感監(jiān)測數(shù)據(jù)。地面監(jiān)測主要包括草地總產(chǎn)草量觀測，觀測集中在2009-2015年7月底至8月初(圖2)。樣地選擇在相應(yīng)群落典型地段且生境條件、植物群落種類組成、群落結(jié)構(gòu)、利用方式和利用強度等具有相對一致性；樣地間具有異質(zhì)性，每個樣地能控制的最大范圍內(nèi)地貌、植被等具有同質(zhì)性，即地貌及植被生長狀況相似。草原植被樣地面積不小于100 hm2。矮小草本和小半灌木草地樣方面積為草原1 m2，荒漠草原2 m2，草甸0.5 m2；株叢小的灌木和高大草本測產(chǎn)面積為4 m×4 m，灌叢或半灌木、高大草本分布均勻的地段設(shè)20 m×5 m或10 m×10 m樣方。設(shè)置一個灌木樣方，灌木樣方內(nèi)設(shè)置一個草本樣方。每個樣方要求測定植物名稱、植物蓋度和草地總產(chǎn)草量。遙感監(jiān)測植被指數(shù)數(shù)據(jù)選用美國國家宇航局250 m的MOD13Q1歸一化植被指數(shù)(NDVI)16 d合成產(chǎn)品，時間是2009年1月至2015年12月，軌道號h23v04、h23v05、h24v04和h24v05，共計644景影像。

1.3 研究方法

1.3.1植被指數(shù)處理 利用MODIS數(shù)據(jù)重投影工具(MODIS reprojection tools，MRT)對MOD13Q1數(shù)據(jù)集進行批量HDF-TIF格式和定義WGS84投影，應(yīng)用GIS平臺計算2009-2015年7月底至8月初多年NDVI均值和最大值合成法(maximum value composition，MVC)合成2009-2015年7月底至8月初最大NDVI，以表征對應(yīng)年草地總產(chǎn)草量實測期內(nèi)牧草平均和最好長勢情況，基于GPS測點信息提取地面采樣點對應(yīng)NDVI。

圖2 2009-2015年7月底至8月初牧草總產(chǎn)量實測位置分布圖Fig. 2 Spatial distribution of the measured locations of the total production of the herbage from the end of July to the beginning of August during the years 2009-2015

1.3.2巴氏距離計算 以縣(市)為研究單元，統(tǒng)計2009-2015年7月底至8月初時段多年NDVI均值頻數(shù)直方圖，利用巴氏距離(Bhattacharyya distance)測量各縣(市)及其臨近縣(市)NDVI直方圖分布的相似性。巴氏距離定義為[19-21]：

(1)

式中：H1與H2為縣(市)尺度多年NDVI均值頻數(shù)直方圖，N為頻數(shù)直方圖組數(shù)，I為頻數(shù)直方圖組編號，H1(I)為縣(市)頻數(shù)直方圖中組數(shù)I對應(yīng)的頻數(shù)。d(H1,H2)越接近1說明兩個直方圖越相似，即所在縣(市)被劃分為一個遙感分區(qū)。本研究中，定義d>0.5，即為相似區(qū)域。

1.3.3模型精度驗證 應(yīng)用留一交叉驗證(leave-one-out cross validation，LOOCV)對草地總產(chǎn)草量遙感模型模擬結(jié)果進行驗證[22]。假設(shè)有n個樣本，從中選擇一個觀測值作為驗證數(shù)據(jù)，其他n-1個樣本作為訓(xùn)練樣本來建立回歸模型，用回歸模型模擬選出的單個驗證數(shù)據(jù)的預(yù)測值來檢驗?zāi)Ｐ途?，如此重?fù)n次，用n個結(jié)果平均值來衡量模型模擬精度。每個模型預(yù)測能力由均方根誤差(RMSE)及擬合系數(shù)R2決定。RMSE用來量化模型精度，R2用來評估模型的準確性。

(2)

式中：Fw(yi)表示第i個樣方草地總產(chǎn)草量(鮮重)實測值，yi為模型模擬第i個預(yù)測值，n是觀測樣本總數(shù)。RMSE越小表明回歸模型越精確，R2越接近于1表示模型精度越高。

2 基于巴氏距離的NDVI均值相似性分析和建模分區(qū)

各縣(市)間NDVI均值巴氏距離計算結(jié)果如表2所列，將巴氏距離大于0.5的縣市定義為相似縣(市)，研究區(qū)域被分為7個建模區(qū)域(圖3)。

各建模區(qū)域所含縣(市)多年NDVI均值-頻數(shù)直方圖如圖4所示，可見在地形、氣候條件、土壤環(huán)境等方面具有很大相似性的區(qū)域應(yīng)該具有相似的牧草生物量，從而使得該區(qū)域NDVI值也相近。而對于地形、氣候條件、土壤環(huán)境等方面空間差異較大的地區(qū)，其牧草生物量也呈現(xiàn)較大的差異性，對應(yīng)區(qū)域NDVI值也具有大的離散度。因此，本研究依據(jù)NDVI值的局部相似性準則進行牧草生物量建模分區(qū)是合理的。

對各個分區(qū)分別采用指數(shù)、線性、多項式以及冪函數(shù)4類函數(shù)建立植被指數(shù)與草地總產(chǎn)草量回歸模型(表4)。結(jié)果顯示：Ⅰ、Ⅲ、Ⅵ分區(qū)選擇冪指數(shù)函數(shù)、Ⅱ、Ⅳ、Ⅶ分區(qū)選擇指數(shù)函數(shù)、Ⅴ分區(qū)采用多項式進行建模具有最小的擬合誤差。這可能是因為不同分區(qū)，氣候、地形、牧草生長環(huán)境等因子存在差異，使得不同分區(qū)NDVI空間分布差異較大，從而使得產(chǎn)草量隨NDVI的變化速率及趨勢也不同。因此，利用NDVI進行建模分區(qū)，不同分區(qū)采用不同的擬合模型才能更精確地估算整個研究區(qū)的牧草產(chǎn)草量。

如圖5a所示，分區(qū)Ⅰ采用冪指數(shù)模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系，效果最佳，能體現(xiàn)牧草鮮重隨NDVI呈明顯的快速遞增現(xiàn)象。該分區(qū)內(nèi)NDVI在0.25～0.85變化，表明該分區(qū)內(nèi)大部分區(qū)域草地長勢較好、產(chǎn)草量較高。圖5b所示，分區(qū)Ⅱ采用指數(shù)模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系效果較好，能夠很好地反映牧草鮮重隨NDVI值增長的變化趨勢。該分區(qū)內(nèi)NDVI主要在0.04～0.65間變化。當少數(shù)子區(qū)域NDVI增加至0.7左右時，對應(yīng)產(chǎn)草量驟然變大。當NDVI在0.04～0.50內(nèi)時牧草鮮重隨NDVI呈近似線性變化。圖5c所示，分區(qū)Ⅲ采用冪指數(shù)模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系效果較好，能夠反映牧草鮮重隨NDVI的變化趨勢，該分區(qū)內(nèi)NDVI主要在0.10～0.47變化，但是，少數(shù)子區(qū)域產(chǎn)草量超出分區(qū)產(chǎn)量均值，異常大。如圖5d所示，分區(qū)IV采用指數(shù)模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系，誤差最小，該分區(qū)內(nèi)NDVI值大都大于0.25。但是，存在一些子區(qū)域的產(chǎn)草量異常大。圖5e所示，分區(qū)Ⅴ采用多項式擬合模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系，誤差最小。該分區(qū)內(nèi)NDVI值都大于0.05，牧草鮮重隨NDVI值變化速度較慢。圖5f所示，分區(qū)Ⅵ采用冪指數(shù)模型擬察布查爾Chabhar County，伊寧縣Yining County，霍城縣Huocheng County，伊寧市Yining City，尼勒克縣Nilka County，鞏留縣Gongliu County，特克斯縣Tekesi County，新源縣Xinyuan County，昭蘇縣Zhaosu County，烏蘇市Wusu City，博樂市Bole City，溫泉縣Wenquan County，精河縣Jinghe County，石河子市Shihezi City，沙灣縣Shawan County，瑪納斯縣Manas County，呼圖壁縣Hutubi County，奎屯市Kuitun City，昌吉市Changji City，米泉市Miquan City，阜康市Fukang City，奇臺縣Qitai County，木壘自治縣Mulei Autonomous County，吉木薩爾縣Jimsar County，巴里坤縣Barlikun County，伊吾縣Yiwu County，和靜Hejing County.

圖3 基于巴氏距離的NDVI時間序列相似性分區(qū)示意圖Fig. 3 Diagram of similarity zoning of NDVI time serial data based on the Bhattacharyya distance

圖4 歸一化NDVI均值-頻數(shù)直方圖Fig. 4 The normalized mean-frequency histogram of NDVI data in the study area

合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系，誤差最小。該分區(qū)內(nèi)NDVI值都大于0.1，當NDVI在0.1～0.3牧草鮮重隨NDVI呈近似的線性增加。圖5g所示，分區(qū)Ⅶ采用指數(shù)模型擬合草地總產(chǎn)草量與NDVI之間的函數(shù)關(guān)系是合適的，該分區(qū)內(nèi)NDVI都大于0.3，主要集中在0.3～0.8，牧草鮮重隨NDVI上升的速度較快。

3 討論

3.1 與不分區(qū)建模對比

忽略NDVI的空間異質(zhì)性，采用研究區(qū)內(nèi)所有采樣點建立NDVI-草地總產(chǎn)草量回歸模型(表5)，與考慮NDVI空間異質(zhì)性采用相似性分區(qū)的回歸模型對比，以探究相似性分區(qū)回歸模型在各個區(qū)域建模的提升效果?？梢钥闯觯€性回歸模型決定系數(shù)R2為0.389，其余回歸模型的決定系數(shù)R2均在0.570～0.600，最高的冪指數(shù)模型為0.600。同時，除線性回歸模型的估計誤差RMSE大于2 000 kg·hm-2，為2 205.219 kg·hm-2外，其余回歸模型的估計誤差RMSE均在2 000 kg·hm-2以下，最低為冪指數(shù)模型，為1 478.640 kg·hm-2。結(jié)合R2與RMSE兩個指標，可知適合全區(qū)域產(chǎn)草量估算的回歸模型為冪指數(shù)模型。

表4 NDVI-草地總產(chǎn)草量回歸模型及交叉檢驗結(jié)果列表Table 4 List of regression models and the results of cross validation of total production of the herbage-NDVI

圖5 各個建模分區(qū)草地總產(chǎn)草量-NDVI擬合關(guān)系圖Fig. 5 The fitting relation curves of the total production of the herbage-NDVI for each modeling zone

回歸模型Modeling type回歸系數(shù)Regression coefficientsR2調(diào)整R2 Adjusted correlationcoefficientRMSE/(kg·hm-2)樣點數(shù)Sample numbery=aebx+ca=57.185; b=0.703; c=425.1560.5970.5961 492.482y=ax+ba=1 330.417; b=-2 619.2600.3890.3892 205.219y=a+bx+cx2a=3 820.722; b=-2 682.155; c=482.3200.5730.5711 680.967y=axb+ca=0.430; b=5.035; c=750.2530.6000.5991 478.640600

對比表4中各分區(qū)的最優(yōu)模型以及表5中不分區(qū)情況下的最優(yōu)模型，不分區(qū)回歸模型的決定系數(shù)R2低于各個分區(qū)回歸模型。說明采用NDVI相似性分區(qū)后，回歸模型自變量能解釋更多因變量的變異，更適合用于牧草產(chǎn)量的估算。同時，除分區(qū)Ⅰ外，其余分區(qū)回歸模型的估計誤差RMSE均在1 000 kg·hm-2以內(nèi)，最低可到266 kg·hm-2，均優(yōu)于不分區(qū)回歸模型。但分區(qū)Ⅰ回歸模型的估計誤差達到了2 951.910 kg·hm-2，這主要是因為分區(qū)Ⅰ內(nèi)實測草地總產(chǎn)量在10 000～30 000 kg·hm-2的樣點多；而不分區(qū)回歸模型估計誤差低于分區(qū)Ⅰ的原因是研究區(qū)內(nèi)其余區(qū)域的實測草地總產(chǎn)草量較低，低于10 000 kg·hm-2，由交叉驗證原理可知，較低的實測值可將較高的模型估計誤差平均。因此，分區(qū)回歸模型能在各個區(qū)域內(nèi)提升牧草產(chǎn)量的估算。

3.2 與隨機分區(qū)建模對比

分區(qū)建模能夠在各個區(qū)域內(nèi)很好提升牧草產(chǎn)量的估算精度。根據(jù)采樣點的聚集程度進行隨機分區(qū)(圖6)，并在各個隨機分區(qū)內(nèi)建立NDVI-草地總產(chǎn)草量回歸模型(表5)，以探究相似性分區(qū)回歸建模對牧草產(chǎn)量估算的提升?？梢钥闯?，隨機分區(qū)Ⅰ的最佳回歸模型為冪指數(shù)模型，其決定系數(shù)R2為0.549，估計誤差RMSE為2 508.562 kg·hm-2(表6)；隨機分區(qū)Ⅱ的最佳回歸模型為拋物線模型，其R2為0.524，RMSE為479.584 kg·hm-2；隨機分區(qū)Ⅲ最佳回歸模型為指數(shù)模型，其模型決定系數(shù)R2達到了0.692，模型的估計誤差為604.101 kg·hm-2。

結(jié)合表5和表6，隨機分區(qū)后模型對因變量變異的解釋程度提升不大，只有隨機分區(qū)Ⅲ回歸模型決定系數(shù)大于不分區(qū)回歸模型，這表明分區(qū)建模并不總是帶來好的估算結(jié)果。同時，模型估計精度RMSE也只在隨機分區(qū)Ⅱ與隨機分區(qū)Ⅲ低于不分區(qū)回歸模型，表明分區(qū)建模對產(chǎn)草量估算的提升具有區(qū)域性。結(jié)合表4和表6，采用相似性分區(qū)后，回歸模型在模型決定系數(shù)R2、模型估計精度RMSE兩個指標上都有顯著提升，除了分區(qū)相似性分區(qū)Ⅰ回歸模型的RMSE。

圖6 NDVI時間序列隨機分區(qū)示意圖Fig. 6 Diagram of random zoning of NDVI time serial data

隨機分區(qū) Random zone回歸模型Modeling type回歸系數(shù)Regression coefficientsR2調(diào)整R2 Adjusted correlationcoefficientRMSE/(kg·hm-2)樣點數(shù)Sample numberⅠy=aebx+ca=333.549; b=0.485; c=-690.4080.5470.5442 516.493y=ax+ba=2 383.587; b=-7.565E+030.4700.4682 951.910y=a+bx+cx2a=4 929.915; b=-3 038.495; c=529.2840.5450.5412 552.954y=axb+ca=5.675; b=3.782; c=459.5520.5490.5452 508.562251Ⅱy=aebx+ca=21.182; b=0.720 6; c=686.1250.5290.522501.783y=ax+ba=366.143; b=-82.7850.4040.400540.456y=a+bx+cx2a=1 530.595; b=-662.416; c=133.3700.5240.517497.584y=axb+ca=0.533; b=4.449; c=776.7870.5300.524500.581149Ⅲy=aebx+ca=0.036 0; b=1.802; c=860.0360.6920.689604.101y=ax+ba=615.016; b=-389.0690.3220.318915.760y=a+bx+cx2a=2 441.925; b=-1 643.751; c=316.8540.5200.515821.791y=axb+ca=8.301E-07; b=11.965; c=882.2260.6890.686610.138200

4 結(jié)論

不同類型草地因其種類構(gòu)成、長勢、蓋度、產(chǎn)量差異引起衛(wèi)星植被指數(shù)的變化，使得充分利用衛(wèi)星資料進行大面積周而復(fù)始地監(jiān)測草地生產(chǎn)力動態(tài)變化成為可能[5]。遙感估算天山地區(qū)草地總產(chǎn)量是一項具有重要實用價值的工作。天山北坡草地垂直分異明顯，類型復(fù)雜，給遙感識別分類和監(jiān)測帶來了很大的困難，本研究提出的NDVI相似性分區(qū)下天山山區(qū)草地總產(chǎn)草量遙感估算技術(shù)實用可行，為后期高級估算方法開發(fā)打下了良好基礎(chǔ)，將研究區(qū)域分為若干個分區(qū)可以減少估產(chǎn)模型的異質(zhì)性，使得回歸模型中的誤差項完全符合理想假設(shè)(誤差獨立同分布期望為零，從而提高模型對數(shù)據(jù)的整體擬合度)。根據(jù)NDVI直方圖顯著的空間差異性特征和鄰近區(qū)域特征相似性，分區(qū)域?qū)?009-2015年內(nèi)7-8月牧草總產(chǎn)量(鮮重)和NDVI進行單一響應(yīng)因子建模，發(fā)現(xiàn)不同建模區(qū)域擬合方程形式也不同，總產(chǎn)草量與植被指數(shù)具有線性、冪指數(shù)、指數(shù)關(guān)系、一元二次回歸方程多種擬合關(guān)系?？傮w來看，7個建模分區(qū)NDVI-草地總產(chǎn)草量擬合相關(guān)系數(shù)在0.836～0.784。交叉檢驗除天山北坡西段-伊犁河谷草原畜牧業(yè)區(qū)RMSE值為2 951 kg·hm-2外，其他分區(qū)RMSE值在266～928 kg·hm-2，原因在于伊犁河谷草原畜牧業(yè)區(qū)實測草地總產(chǎn)量在10 000～30 000 kg·hm-2的樣點多。

后續(xù)研究中，會在NDVI分區(qū)建模的基礎(chǔ)上，綜合考慮NDVI分布的巴式距離分區(qū)和NDVI空間分布(NDVI空間結(jié)構(gòu)信息)，增加確定性自變量的一般回歸模型和帶有觀測誤差自變量的擴展回歸模型(遙感反演誤差到回歸統(tǒng)計誤差的傳播)的比較分析，加入草地長勢其他有關(guān)解釋變量進行建模，設(shè)計區(qū)域差異、區(qū)域效應(yīng)在回歸系數(shù)中，如空間濾波和空間面板回歸模型里的空間或個體效應(yīng)變量。