基于Landsat和MODIS數據融合的農牧區NPP模擬

尹小君 祝宏輝 GAO Jerry 高 軍 郭麗潔 茍貞珍

(1.石河子大學信息科學與技術學院， 石河子 832003； 2.兵團空間信息工程技術研究中心， 石河子 832003；3.石河子大學經濟管理學院， 石河子 832003； 4.圣何塞州立大學計算機工程系， 圣何塞 CA 95192；5.安徽四創電子股份有限公司， 合肥 230094)

0 引言

植被NPP是指綠色植物經過光合作用在單位時間、單位面積所獲取的有機物質總量與維持自養呼吸所消耗量的差值[1-3]，它是植被自身生理特性與外界環境共同作用的結果[4-5]，能夠直觀反映地表植被在外界環境作用下的生產能力[6]，是評判陸地生態系統健康狀況以及研究碳源/匯過程的重要指標[7-9]。天山北坡是我國西北干旱區的主體，同時也是我國重要的農牧業發展基地[10]，由于近年來氣候變化、人類活動以及土地利用變更等因素的影響，已造成天山北坡區域內草地退化、林地萎縮、地下水位下降等一系列生態環境問題[11-13]。因此，適時進行區域內植被NPP的動態模擬，掌握植被生長健康狀況，對研究植被動態變化、合理分配農牧場草地資源至關重要。

隨著遙感技術的不斷發展，遙感數據以其成像周期短、覆蓋范圍廣、成本低廉等特點，可快速、準確地獲取植被NPP的動態變化信息，為植被生長健康狀況監測提供了有效途徑。目前，基于遙感的植被NPP模擬模型主要分為3類[14]：氣候生產力模型、光能利用率模型和植被生理生態過程模型。其中，基于光能利用率的CASA(Carnegie-Ames-Stanford Approach)模型可操作性強，在不同地域、不同研究尺度的NPP模擬中得到了廣泛應用[15-17]。自2000年以來，MODIS數據因其時間分辨率高，能夠較好地反映區域內NPP時間序列變化信息，而成為研究植被NPP動態變化的主要數據源[18-20]。但MODIS數據受到空間分辨率的限制，導致植被變化的空間細節信息表達不夠精細，難以滿足在精細尺度下的NPP模擬。為此，學者們開始尋求中分辨率遙感數據對區域植被NPP進行模擬。池源等[21]基于Landsat 8 OLI影像探討了廟島群島北五島景觀格局與植被NPP的關系；管小彬等[22]基于TM/ETM+影像分析了武漢市冬季植被NPP時空變化特征；楊會巾等[23]基于Landsat 8 OLI對瑪納斯河流域植被NPP進行了模擬。然而，由于中分辨率衛星傳感器重訪周期較長，且易受到云、雨、雪等自然天氣的影響，使得獲取連續時間序列的NPP數據受到了限制。

為滿足大范圍、高精度及快速變化的地表信息遙感監測需求，多源遙感數據時空融合技術為解決上述問題提供了有效的途徑。目前，國內外學者已將多源遙感數據時空融合技術廣泛應用于不同領域[24-26]，因此，利用該方法將不同時間分辨率、空間分辨率和光譜分辨率的遙感數據進行時空融合，獲取同時具有高時間分辨率遙感數據的時相信息和中空間分辨率數據的空間紋理信息，可為陸地植被NPP的時空動態變化研究提供有效的數據源。本文以天山北坡中段區域為研究區，探討利用多源遙感數據時空融合技術獲取中空間分辨率遙感影像，模擬研究區2016年的植被NPP，并對融合數據的有效性進行精度檢驗，以期為多源遙感影像融合技術與光能利用率模型協同模擬NPP提供新的思路。

1 研究區概況及數據獲取

1.1 研究區概況

研究區位于新疆天山北坡中段綠洲核心區，地理位置為43°40′～45°33′N，84°72′～86°56′E，海拔260～4 900 m，總面積為2.38×104km2，是新疆典型的農牧業發展基地，如圖1所示。研究區地處亞歐大陸腹地，遠離海洋，屬于典型的溫帶大陸性氣候，區域內四季分明，夏季炎熱干燥，溫差較大，蒸發強烈，冬季寒冷漫長，降水量稀少，生態系統比較脆弱。

圖1 研究區行政區劃圖Fig.1 Administrative map of study area

研究區植被類型分布主要與海拔有關，因不同的海拔梯度，植被整體呈現垂直結構分布。在研究區低海拔山麓平原區，主要分布荒漠耐旱植被。海拔2 000 m左右主要分布低山草原，是研究區主要的植被類型，也是家庭牧場的主要聚集區，其中森林與草甸相間出現。高海拔區域，主要分布一些高山草甸和苔蘚地衣。

1.2 數據及處理

Landsat8是NASA于2013年成功發射的衛星，有2個主要載荷：OLI和TIRS。其中OLI為陸地成像儀，TIRS為熱紅外傳感器。選取2016年研究區內數據質量較好，云層覆蓋度小于5%的Landsat 8 OLI影像，空間分辨率為30 m，軌道重訪周期16 d。影像行列號為144/29和144/30，數據來源于中國科學院計算機網絡信息中心的“地理空間數據云平臺”(http:∥www.gsclud.cn)。利用Landsat 8 OLI 近紅外波段反射率和紅光波段反射率，獲取NDVI數據[27]。

MODIS是美國進行對地觀測系統戰略實施發射的重要光譜傳感器，該數據不僅可以免費獲取，而且數據時間序列范圍廣，更新頻率快，能夠較好地獲取快速變化的自然植被等信息。MODIS搭載在Terra和Aqua 2顆衛星上，每天獲取最少2次白天和2次黑夜MODIS數據。MODIS陸地標準數據產品具有每日產品、8 d合成產品、16 d合成產品、月合成產品、季度合成產品和年合成產品。本文使用的MODIS數據為MODIS13Q1產品中的NDVI數據，空間分辨率為250 m，時間分辨率為16 d。數據來源于美國陸地過程分布式數據檔案中心(https:∥lpdaac.usgs.gov/)。影像獲取時間為2016年1月1日—12月18日，共23期，行列號為h23v04、h24v04。數據預處理首先使用MRT(MODIS reprojection tool)軟件對影像進行鑲嵌、格式轉換和重投影等操作，再利用ArcGIS軟件對研究區范圍進行裁剪。氣象數據選取研究區內及周邊氣象站點的日平均溫度、日降水總量和日太陽總輻射數據資料，該數據來源于中國氣象數據網(http:∥data.cma.cn/)，將所得日數據信息合成月數據信息，利用澳大利亞科學家Hutchinson基于薄盤樣條理論編寫的專業氣象數據插值工具ANUSPLIN[28]，以經緯度為自變量，高程作為協變量，對研究區內氣象數據進行空間插值，獲取氣象柵格數據。

研究區地形數據采用SRTM 30 m DEM產品，數據來源于“地理空間數據云平臺”(https:∥www.gscloud.cn)。土地利用覆被數據來源于國家地球系統科學數據共享網服務平臺(http:∥www.geodata.cn/)。所有柵格數據的行列數、像元大小及投影方式均與Landsat 8 OLI數據一致。

于2016年7月18日、7月26日、8月9日、8月18日和8月26日，在天山北坡中段農牧區，即在牧場草地生長期對草地地上生物量進行采集，根據DEM數據掌握當地的地形，分別在陽坡和陰坡設置采樣點140個，如圖2所示。根據實際情況選取一定數量的樣方，樣方選取為1 m×1 m，在實驗區每隔300 m重復的點采樣3次，其平均值作為樣本值，共計46個樣本。齊地剪取樣方地上生物量，剪取樣本的同時用手持GPS定位儀，記下單位樣方的經緯度，樣本裝袋并編號統一帶回實驗室稱量。樣本數據的詳細信息主要包括生物量干質量、濕質量、樣點坐標、高程、覆蓋度、經緯度、海拔、時間、草地群落類型等。

圖2 采樣點空間分布圖Fig.2 Spatial distribution map of sampling points

2 研究方法

基于Landsat數據與MODIS數據，采用時空融合技術實現對研究區2016年NPP的模擬。首先將研究區遙感數據進行預處理，獲取研究區內不同數據源的NDVI數據，然后利用STARFM數據融合算法得到研究區中空間分辨率時間序列的NDVI(30 m)數據，其次結合已處理好的氣象數據(溫度、降水、太陽總輻射)、土地利用覆被數據采用CASA模型對研究區2016年植被NPP進行模擬，最后根據研究區實測數據對估測NPP結果進行精度驗證，其具體流程如圖3所示。

2.1 遙感數據時空融合

采用的遙感數據時空融合算法為時空適應性反射率融合算法模型(Spatial and temporal adaptive reflectance fusion model，STARFM)[25]，該方法是基于t0時刻的中空間分辨率影像與低空間分辨率MODIS影像，結合t1時刻低空間分辨率MODIS影像，通過計算像元間的空間分布差異及波譜差異，模擬預測t1時刻中空間分辨率數據的方法。基于t0時刻獲取的中空間分辨率Landsat 8 OLI NDVI數據(2016年6月24日)和低空間分辨率MODIS NDVI數據(2016年6月25日)以及t1時刻獲取的低空間分辨率MODIS NDVI數據(2016年4月6日和2016年8月26日)，預測t1時刻Landsat 8 OLI NDVI 數據(2016年4月6日和2016年8月26日)，預測過程通過使用滑動窗口法減少低分辨率遙感數據像元邊界的影響[29]。其計算公式為

(1)

其中

(2)

式中L(xw/2,yw/2,t1)——預測t1時刻的Landsat 8 OLI NDVI像元值

w——移動窗口尺寸

xw/2、yw/2——窗口中心的像元坐標

M(xi,yi,t1)——窗口位置(xi,yi)處在t1時刻的MODIS像元NDVI值

L(xi,yi,t0)——t0時刻給定位置(xi,yi)中Landsat像元的NDVI值

M(xi,yi,t0)——t0時刻給定位置(xi,yi)中MODIS像元的NDVI值

Wijk——移動窗口內各像元對預測中心像元值權重

cijk——由移動窗口內的中心像元與窗口內其他像元的光譜距離、時間距離和空間距離共同確定的系數

2.2 CASA模型

CASA(Carnegie-Ames-Stanford Approach)模型是由遙感、氣象、土地覆被數據共同驅動的光能利用率模型，是目前NPP模擬模型中應用最廣且較為成熟的模型[30]。采用朱文泉等[31]改進的CASA模型對天山北坡中段植被NPP進行模擬，計算公式為

NPP(x,t)=APAR(x,t)E(x,t)

(3)

式中NPP(x,t)——像元x在t月的NPP,g/(m2·month)

APAR(x,t)——像元x在t月植被所吸收的光合有效輻射，MJ/(m2·month)

E(x,t)——像元x在t月的光能利用率，g/MJ

植被吸收的光合有效輻射(APAR)取決于太陽總輻射量(SOL)和植被對APAR的吸收比例(FPAR)，計算公式為

APAR(x,t)=0.5SOL(x,t)FPAR(x,t)

(4)

式中SOL(x,t)——像元x在t月的太陽總輻射量，MJ/m2

FPAR(x,t)——植被冠層對光合有效輻射的吸收比例

在一定范圍內FPAR與歸一化植被指數(NDVI)、比值植被指數(SR)都存在良好的線性關系[32]，可以通過遙感影像提取植被的NDVI來對FPAR進行模擬，計算公式為

(5)

FPAR(x,t)SR=

(6)

(7)

式中NDVIi,max——對應第i種植被類型的NDVI最大值

NDVIi,min——對應第i種植被類型的NDVI最小值

SRmax——對應第i種植被類型的SR最大值

SRmin——對應第i種植被類型的SR最小值

FPARmax、FPARmin——常數，分別為0.950、0.001

FIELD等[33]研究發現通過這2種方法模擬所得的結果各有不同，由NDVI模擬的FPAR結果比實測值偏高，而由SR模擬所得的FPAR結果比實測值偏低，因此，可以利用調整系數0.5，將兩者取平均值的方式作為對FPAR的模擬結果，從而使其得到的FPAR模擬誤差最小，計算公式為

(8)

天山北坡中段為洪積沖積扇，呈荒漠-綠洲-山地景觀格局，植被覆蓋稀疏，把調節系統調整為0.52。

光能利用率是指綠色植被把所吸收的光合有效輻射轉換為有機碳的效率，它主要與溫度及水分脅迫系數有關，其計算公式為

E(x,t)=TE1(x,t)TE2(x,t)WE(x,t)εmax

(9)

式中TE1(x,t)——低溫作用下對植被光能轉換率的脅迫系數

TE2(x,t)——高溫作用下對植被光能轉換率的脅迫系數

WE(x,t)——水分條件對植被光能轉換率的影響系數

εmax——最大光能轉換率

3 結果與分析

3.1 數據融合精度

為驗證Landsat與MODIS數據采用STARFM模型融合生成NDVI數據的精度，本文以2016年6月24日的Landsat 8 OLI NDVI數據和2016年6月25日的MODIS NDVI數據為基礎，分別與2016年4月6日和2016年8月12日的MODIS NDVI數據進行融合，得到融合后的2016年4月6日與2016年8月12日的NDVI數據，空間分辨率為30 m。然后，將融合后的兩期NDVI數據分別與真實的Landsat 8 OLI NDVI數據進行對比分析。如圖4所示，通過分析可知基于STARFM融合算法融合的NDVI數據相比于原MODIS NDVI數據，空間分辨率得到顯著提高，且具有良好的空間細節信息，能夠較好地表達較小地物之間空間差異，與真實的Landsat 8 OLI NDVI數據相比，具有良好的空間分辨率，且空間分布趨勢基本一致。

圖4 MODIS、Landsat 8 OLI與融合的NDVI對比分析Fig.4 Comparative analysis of NDVI under different data sources

為驗證融合數據與真實數據之間的差異，除了通過目視判別融合效果外，還可通過相關系數r、偏差Bias、均方根誤差RMSE等指標對影像融合精度進行定量評價[34-35]，其偏差計算公式為

(10)

式中N——像元個數i——像元序號

xi——融合數據像元值

yi——真實數據像元值

利用2016年多期融合的空間分辨率為30 m的NDVI數據，與對應時刻的空間分辨率為30 m的Landsat 8 OLI NDVI真實影像進行對比(表1)，檢驗結果表明，融合生成的NDVI數據與真實NDVI數據之間的r不小于0.759，Bias在0.006 2～0.009 4之間，RMSE在0.074～0.135之間，因此，說明利用該融合算法生成的NDVI數據精度較高，與真實影像具有較好的一致性，能夠用于后續植被NPP的模擬。

3.2 NPP模擬結果及精度驗證

3.2.1NPP模擬結果

利用CASA模型模擬研究區2016年植被NPP空間分布，結果如圖5所示，基于MODIS NDVI數據模擬的植被NPP平均值為327.28 g/(m2·a)，基于融合的NDVI數據模擬的NPP平均值為305.47 g/(m2·a)，兩者空間上整體呈現南高北低的趨勢，沿天山山脈一側植被的NPP平均值最高，原因主要為天山上的融雪積水為植被的生長提供了必備的有機質，相對較高的海拔使得人類活動對植被的干擾強度減弱，因此NPP均值較高?；诓煌瑪祿措m在模擬研究區內植被NPP空間分布狀況整體相差不大，但是由于空間分辨率的較大差異，導致在反映植被內部變化信息時差異較為明顯，整體上，基于MODIS NDVI數據模擬的NPP結果由于空間分辨率較低，導致研究區內存在大量的混合像元，難以清晰地表達不同地物間的空間紋理細節信息，而基于融合的NDVI數據模擬的NPP結果空間分辨率相對較高，影像內存在混合像元的數量大大減少，空間紋理信息較清晰，能夠很好地表達微小地物間的空間差異?；贛ODIS NDVI數據模擬的NPP像元大部分集中在240～660 g/(m2·a)之間，變化幅度不明顯，而基于融合數據模擬的NPP像元大部分集中在320 g/(m2·a)左右，整體近似于正態分布，能夠清晰地表達農田、草地、建筑物、道路、河流等不同地物間的NPP差異。

表1 融合數據與真實數據間的相關系數、偏差及 均方根誤差Tab.1 Correlation, bias and RMSE between fused data and true data

圖5 研究區2016年NPP模擬結果Fig.5 Simulation of NPP in 2016

3.2.2融合數據模擬NPP精度驗證

利用2016年在研究區內采集的46個樣本生物量數據，驗證基于融合后的NDVI數據模擬研究區植被NPP的精度，由于實測NPP難度較大，因此通常采用實測生物量換算為NPP數據，根據生物量與NPP之間的轉換系數0.475，進而得到采樣點實測NPP數據[31]。由于采樣時間為7、8月，年度植被最為旺盛的時期，其NPP可以代表年度NPP[36]，將實測數據的樣本生物量位置信息與CASA模型模擬的NPP位置信息進行對比，如圖6所示。 結果表明，利用融合數據模擬的NPP模擬值與NPP實測值基本吻合(R2=0.860 1，P<0.01)。因此利用Landsat與MODIS數據的融合結果能夠較好地對研究區植被NPP進行模擬。

圖6 NPP實測值與NPP模擬值對比分析Fig.6 Comparison analysis between simulated NPP and observed NPP

4 結論

(1)利用STARFM時空融合算法，將MODIS NDVI數據與Landsat 8 OLI NDVI數據時空信息進行融合，獲得中空間分辨率NDVI時間序列數據，既能在時間上較好地保留MODIS數據時間變化信息，又能在空間上反映Landsat數據的空間細節信息，與8個時期真實的NDVI數據相比，相關系數r不小于0.759，偏差在0.006 2～0.009 4之間，均方根誤差RMSE在0.074～0.135之間，能夠較好地模擬植被NPP。

(2)運用CASA模型模擬NPP，MODIS數據模擬的NPP像元大部分集中在240～660 g/(m2·a)之間，變化幅度不明顯，而基于融合數據模擬的NPP像元大部分集中在320 g/(m2·a)左右，整體近似于正態分布，能夠清晰表達農田、草地、建筑物、道路、河流等不同地物間的NPP差異。

(3)利用野外實測NPP值進行遙感數據融合與CASA模型協同精度驗證，研究區融合后模型植被NPP與實測NPP數據具有較好的相關性，R2=0.860 1(P<0.01)。