基于BEPS模型的塞罕壩植被凈初級生產(chǎn)力時空變化分析研究

包志意，范文義

摘要：葉面積指數(shù)（leaf area index ，LAI）是北部生態(tài)系統(tǒng)生產(chǎn)力模擬模型（boreal ecosystem productivity simulator，BEPS）的關(guān)鍵驅(qū)動數(shù)據(jù)，獲取高精度LAI對區(qū)域森林生態(tài)系統(tǒng)碳循環(huán)十分重要，然而當(dāng)前大多研究采用的MODIS LAI產(chǎn)品缺乏可信度。為此，基于LAI動態(tài)模型、PROSAIL輻射傳輸模型和層狀貝葉斯網(wǎng)絡(luò)（Hierarchical Bayesian Network， HBN）構(gòu)建數(shù)據(jù)同化系統(tǒng)，獲得空間分辨率為20 m的LAI數(shù)據(jù)，驅(qū)動BEPS模型，模擬塞罕壩機(jī)械林場2011—2021年的植被凈初級生產(chǎn)力（Net Primary Productivity， NPP），并對NPP時空變化特征及影響因子進(jìn)行分析。結(jié)果表明，基于貝葉斯同化方法獲得的高分辨率LAI數(shù)據(jù)極大提高了MODIS LAI產(chǎn)品的精度；基于同化后的LAI數(shù)據(jù)驅(qū)動BEPS模型獲取模擬森林NPP，與樣地實(shí)測數(shù)據(jù)計(jì)算NPP間相關(guān)性較高（R2=0.77）；2011—2021年塞罕壩機(jī)械林場植被NPP平均值為307.4 g/（m2·a），森林NPP呈現(xiàn)平穩(wěn)增長趨勢；不同植被類型模擬NPP存在較大差異，針葉林、落葉林及混交林模擬NPP分別為484.9、402.4、287.9 g/（m2·a）；植被NPP與溫度因子相關(guān)性較高，偏相關(guān)系數(shù)為0.2～0.8，而植被NPP與降水量的相關(guān)性總體而言相對較低，其偏相關(guān)系數(shù)為-0.3～0.4，在該地區(qū)降水量對植被NPP的影響較低，溫度為該地區(qū)NPP變化的主導(dǎo)因子。研究結(jié)果可獲取高空間分辨率的LAI數(shù)據(jù)，為森林生態(tài)系統(tǒng)碳循環(huán)的精準(zhǔn)時空模擬提供依據(jù)。

關(guān)鍵詞：塞罕壩；葉面積指數(shù)；BEPS模型；植被凈初級生產(chǎn)力

中圖分類號：S718.55文獻(xiàn)標(biāo)識碼：A文章編號：1006-8023（2024）01-0074-11

Analysis of Spatial and Temporal Variation of Vegetation Net Primary Productivity in Saihanba Based on BEPS Model

BAO Zhiyi1，F(xiàn)AN Wenyi1，2*

（1.College of Forestry， Northeast Forestry University， Harbin 150040， China；2.Key Laboratory of Forest Ecosystem Sustainable Management of Ministry of Education， Northeast Forestry University， Harbin 150040， China）

Abstract：

Leaf area index （LAI） is the key driving data of BEPS model， and it is important to obtain high accuracy LAI for regional forest ecosystem carbon cycle， however， the MODIS LAI products used in most current studies lack credibility. To this end， this study constructed a data assimilation system based on LAI dynamic model， PROSAIL radiative transfer model and Hierarchical Bayesian Network （HBN） to obtain LAI data with a spatial resolution of 20 m to drive the BEPS model and simulate the vegetation net primary productivity （NPP） of Saihanba Mechanical Forest during 2011—2021， and the spatial and temporal variation of NPP and the influencing factors of NPP were analyzed. The results showed that the high-resolution LAI data obtained based on Bayesian assimilation method greatly improved the accuracy of MODIS LAI products; the correlation between the simulated forest NPP obtained from BEPS model driven by assimilated LAI data and the NPP calculated from the sample plots was high （R2=0.77）; the mean value of vegetation NPP in Saihanba Mechanical Forest during 2011—2021 was 307.4 g/（m2·a）， and the NPP of forest showed a steady growth trend; the simulated NPP of different vegetation types were different， and the simulated NPP of coniferous， deciduous and mixed forests were 484.9， 402.4， 287.9 g/（m2·a）; the correlation between vegetation NPP and temperature factor was high， and the bias correlation coefficient was 0.2-0.8， while the correlation between vegetation NPP and precipitation was relatively low in general， with bias correlation coefficients of -0.3-0.4. The influence of precipitation on vegetation NPP was low in this region， and temperature was the dominant factor of NPP variation in this region. In this study， high spatial resolution LAI data were obtained to provide a basis for accurate spatial and temporal simulation of the carbon cycle in forest ecosystems.

Keywords：Saihanba; leaf area index; BEPS model; vegetation net primary productivity

0引言

植被凈初級生產(chǎn)力（Net Primary Productivity，NPP）定義為在單位時間、單位面積內(nèi)植被利用光合作用所固定的有機(jī)物總量與自養(yǎng)呼吸消耗的差值，通過植被與外界環(huán)境因子相互作用下所產(chǎn)生[1]。NPP是反映植被的生產(chǎn)能力和生態(tài)過程的指標(biāo)，用以評價(jià)生態(tài)系統(tǒng)植被質(zhì)量狀況和區(qū)域生態(tài)系統(tǒng)功能[2-3]。

模擬NPP的方法包括實(shí)測法和模型模擬法。實(shí)測法的優(yōu)點(diǎn)是可以準(zhǔn)確估算出樣地NPP，但難以實(shí)現(xiàn)對大區(qū)域NPP的估算。而模型模擬法可以結(jié)合遙感數(shù)據(jù)，實(shí)現(xiàn)對大區(qū)域NPP的估算[4]。模型模擬可歸納為統(tǒng)計(jì)模型、生態(tài)系統(tǒng)過程模型和光能利用率模型3大類[5]。統(tǒng)計(jì)模型估算NPP比較簡單但誤差相對較大[6]，光能利用率模型無法解釋NPP變化的生理生態(tài)學(xué)機(jī)制，而生態(tài)過程模型能準(zhǔn)確詳細(xì)地描述植被的生理機(jī)制問題，模擬結(jié)果也較為準(zhǔn)確。生態(tài)過程模型中目前應(yīng)用廣泛且具有代表性的為BEPS模型（boreal ecosystem productivity simulator）。BEPS模型是在FOREST-BGC模型基礎(chǔ)上進(jìn)行不斷改進(jìn)，發(fā)展形成的一種應(yīng)用廣泛的模型，最初使用1 km分辨率的MODIS數(shù)據(jù)為主要驅(qū)動數(shù)據(jù)模擬了加拿大陸地凈初級生產(chǎn)力，解決了時間和空間上的尺度轉(zhuǎn)換問題[7]。近年來，BEPS模型在植被NPP模擬方面得到了廣泛應(yīng)用，如Feng等[8]利用BEPS模型對全國范圍內(nèi)的植被生態(tài)環(huán)境凈初級生產(chǎn)力進(jìn)行模擬。王培娟等[9]利用地形修正對BEPS模型改進(jìn)后，對長白山林區(qū)凈初級生產(chǎn)力進(jìn)行了模擬估測。Mo等[10]利用集合卡爾曼濾波改進(jìn)BEPS模型參數(shù)，成功降低了模型模擬的誤差。毛學(xué)剛等[11]通過平滑處理的MODIS LAI產(chǎn)品結(jié)合日步長的氣象數(shù)據(jù)以及土壤數(shù)據(jù)驅(qū)動BEPS模型，模擬了東北林區(qū)森林生態(tài)系統(tǒng)碳循環(huán)，并在研究中以不同方法對模型模擬精度進(jìn)行驗(yàn)證[12-13]。

在以往的森林生產(chǎn)力模擬研究中，由于用于驅(qū)動模型的LAI數(shù)據(jù)缺乏可信度，導(dǎo)致模型模擬結(jié)果存在較大誤差[14]。葉面積指數(shù)（Leaf Area Index，LAI）定義為單位地表面積綠葉總面積的一半，對植物生理過程以及生態(tài)系統(tǒng)生產(chǎn)力的形成有著重要影響[15-16]，是碳循環(huán)模擬的重要驅(qū)動數(shù)據(jù)[17-18]，獲得高分辨率時間序列LAI數(shù)據(jù)是準(zhǔn)確模擬區(qū)域碳循環(huán)的關(guān)鍵。但從單一時相遙感數(shù)據(jù)中反演LAI，難以獲取長時間序列LAI[17]，MODIS衛(wèi)星遙感提供LAI產(chǎn)品又存在噪聲大、在時間序列上波動較大等缺陷。因此，獲取高精度的LAI數(shù)據(jù)對于區(qū)域生態(tài)系統(tǒng)生產(chǎn)力的精確模擬變得尤為重要。越來越多的學(xué)者采用數(shù)據(jù)同化方法獲取高分辨率LAI數(shù)據(jù)，如李雪建等[19]通過雙集合卡爾曼濾波同化MODIS LAI時間序列數(shù)據(jù)，極大地提高了MODIS LAI的產(chǎn)品精度。Xing等[20]利用HBN同化1 km分辨率的MODIS LAI和250 m分辨率的MODIS反射率數(shù)據(jù)，獲得了高精度多分辨率的LAI數(shù)據(jù)。

綜上所述，LAI數(shù)據(jù)作為BEPS模型的一個關(guān)鍵輸入?yún)?shù)，在以往研究中大多采用MODIS LAI產(chǎn)品，其存在空間分辨率低、混合像元較多以及區(qū)域應(yīng)用中出現(xiàn)明顯高估或低估現(xiàn)象等問題，如何獲取可信度和時空分辨率更高的LAI數(shù)據(jù)變得尤為重要。因此，本研究以塞罕壩機(jī)械林場為研究對象，以2021年MODIS LAI時間序列產(chǎn)品和Sentinel-2反射率數(shù)據(jù)為數(shù)據(jù)源，基于LAI動態(tài)模型、PROSAIL輻射傳輸模型和層狀貝葉斯網(wǎng)絡(luò)（Hierarchical Bayesian Network， HBN）構(gòu)建數(shù)據(jù)同化系統(tǒng)，獲得LAI高時空分辨率數(shù)據(jù)，然后將同化的高時空分辨率LAI作為輸入數(shù)據(jù)驅(qū)動BEPS模型，實(shí)現(xiàn)對塞罕壩地區(qū)植被凈初級生產(chǎn)力時空變化及影響因素的精準(zhǔn)模擬與分析，為植被生態(tài)系統(tǒng)碳循環(huán)的時空模擬提供更精確的葉面積指數(shù)數(shù)據(jù)，為今后塞罕壩森林合理經(jīng)營管理、科學(xué)規(guī)劃提供技術(shù)支撐。

1研究區(qū)與數(shù)據(jù)源

1.1研究區(qū)概況

塞罕壩機(jī)械林場位于河北省承德市圍場滿族蒙古族自治縣最北部區(qū)域，東邊與圍場縣4個鄉(xiāng)相連；南邊與河北省御道口牧場相接壤；北邊與克什克騰旗，西邊與多倫縣相連。中心地理位置（42°22′～42°31′N，116°53′～117°31′E），海拔1 010 ～1 939.9 m。塞罕壩地處典型的森林草原交錯帶，主要樹種為樟子松（Pinus sylvestris var. mongolica）、碩樺（Betula costata）、油松（Pinus tabuliformis）、云杉（Picea asperata）、柞樹（Quercus mongolica）、椴樹（Tilia tuan）、華北落葉松（Larix gmelinii）、山楊（Populus davidiana）和色木槭（Acer pictum）等，如圖1所示。

1.2遙感數(shù)據(jù)收集與處理

1.2.1MODIS數(shù)據(jù)

MOD15A2是500 m空間分辨率的MODIS LAI產(chǎn)品，每8 d合成1次，每年有46個數(shù)據(jù)點(diǎn)。在NASA網(wǎng)站（https：//ladsweb.modaps.eosdis.nasa.gov/）下載了2021年全年的MODIS LAI產(chǎn)品。通過使用MODIS重投影工具（MRT）軟件，將其重投影至WGS-84坐標(biāo)系，空間分辨率為500 m×500 m。使用ENVI 5.3對研究區(qū)域進(jìn)行裁剪，并提取了LAI值。最后，使用Savitzky-Golay（SG）濾波算法，以減少數(shù)據(jù)噪聲和消除異常。SG平滑后的LAI（MODIS_SG_LAI）被用作貝葉斯層狀網(wǎng)絡(luò)（HBN）的初始輸入數(shù)據(jù)。

1.2.2Sentinel-2數(shù)據(jù)

Sentinel-2是由2顆衛(wèi)星組成的星座，帶有多光譜儀器（MSI），MSI對13個光譜波段進(jìn)行采樣。4個具有10 m空間分辨率的波段，6個具有20 m空間分辨率的波段和3個具有60 m空間分辨率的波段，3個植被紅邊波段在準(zhǔn)確監(jiān)測植被生長和獲得相關(guān)參數(shù)方面起著主導(dǎo)作用。哨兵-2可以在哥白尼開放存取網(wǎng)站上免費(fèi)下載。本研究中下載了2021年共計(jì)16組數(shù)據(jù)。并由Sen2cor工具（http：//step.esa.int/main/snap-supported-plugins/sen2cor/）進(jìn)行大氣校正，以減少采集時大氣條件的影響。采集時的大氣條件的影響降到最低，利用SNAP（http：//step.esa.int/main/download/snap-download/）工具對數(shù)據(jù)進(jìn)行重新取樣，用最近的鄰居方法將數(shù)據(jù)取為20 m的像素大小，并重新投影到WGS-84坐標(biāo)系。

1.2.3氣象數(shù)據(jù)

本研究所使用的氣象數(shù)據(jù)為逐日降水量、最高溫度、最低溫度、太陽輻射以及相對濕度等。以上數(shù)據(jù)由“國家青藏高原科學(xué)數(shù)據(jù)中心”下載。對氣象數(shù)據(jù)插值采用的是克里金插值的方法，得到塞罕壩機(jī)械林場區(qū)域2010、2015、2021的1 km空間分辨率的逐日氣象數(shù)據(jù)插值數(shù)據(jù)。

1.2.4土壤有效持水量數(shù)據(jù)

數(shù)據(jù)采用的是南京土壤研究所（http：//www.soil.csdb.cn）制作的土壤類型數(shù)據(jù)，投影方式為等積圓錐投影。首先將中圖土壤類型數(shù)據(jù)進(jìn)行裁切獲得塞罕壩土壤類型矢量圖，然后把土壤類型矢量圖根據(jù)中國土壤類型分類編碼與美國土壤分類系統(tǒng)編碼的對應(yīng)關(guān)系轉(zhuǎn)換為美國土壤分類系統(tǒng)分類標(biāo)準(zhǔn)的土壤質(zhì)地矢量圖，其次，將轉(zhuǎn)化后的土壤質(zhì)地類型圖根據(jù)土壤有效持水量與土壤質(zhì)地間的聯(lián)系把土壤質(zhì)地圖轉(zhuǎn)換為土壤有效持水量圖，最后將土壤有效持水量矢量數(shù)據(jù)轉(zhuǎn)換成空間分辨率為1 km 柵格數(shù)據(jù)。

1.2.5模型精度驗(yàn)證數(shù)據(jù)

驗(yàn)證數(shù)據(jù)來自塞罕壩機(jī)械林場2020年二類調(diào)查數(shù)據(jù)以及2021年外業(yè)測量數(shù)據(jù)。對樣地內(nèi)樹木進(jìn)行每木檢尺調(diào)查，測量并記錄胸徑和樹高等樣地信息。

1.2.6土地覆蓋數(shù)據(jù)

研究區(qū)土地覆蓋類型數(shù)據(jù)是模型重要輸入數(shù)據(jù)，其決定著不同植被類型的生理參數(shù)。本研究中的土地覆蓋類型數(shù)據(jù)采用中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心（https：//www.resdc.cn/）公開的土地利用現(xiàn)狀遙感監(jiān)測數(shù)據(jù)，分辨率為30 m。將研究區(qū)地類分為耕地、草地、水域、結(jié)社用地、裸露土地、落葉林、針葉林和混交林，如圖2所示。

2研究方法

本研究采用的北部生態(tài)系統(tǒng)生產(chǎn)力模擬模型（boreal ecosystem productivity simulator， BEPS）是在FOREST-BGC[21]模型基礎(chǔ)上發(fā)展起來的遙感機(jī)理模型[22-23]，BEPS模型主要由土壤光合作用模型、氣孔導(dǎo)度模型、水分平衡模型和自養(yǎng)呼吸作用模型等4個部分組成。Liu等[18，24-25]和Bruand等[22]后來對BEPS模型繼續(xù)完善。時間尺度擴(kuò)展是利用葉片尺度的瞬時光合作用將分為陽葉和陰葉的森林冠層葉片通過日積分計(jì)算日光合作用實(shí)現(xiàn)， 空間尺度擴(kuò)展是基于將葉片分為陽葉和陰葉并分別模擬葉片的輻射收支實(shí)現(xiàn)，獲取植被總初級生產(chǎn)力（GPP）與植被自養(yǎng)呼吸的差值等于植被NPP[26-28]。

研究中LAI數(shù)據(jù)的同化是以空間分辨率為500 m的MODIS LAI 和空間分辨率為20 m的sentinel-2反射率數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，由LAI動態(tài)模型，PROSAIL模型與層狀貝葉斯網(wǎng)絡(luò)算法組成的同化系統(tǒng)實(shí)現(xiàn)高空間分辨率LAI數(shù)據(jù)的同化。然后將同化后的LAI數(shù)據(jù)結(jié)合氣象數(shù)據(jù)、土地覆蓋類型數(shù)據(jù)、土壤有效持水量數(shù)據(jù)輸入并驅(qū)動BEPS模型，模擬塞罕壩植被凈初級生產(chǎn)力。

2.1LAI同化方法

2.1.1LAI動態(tài)模型

本研究中使用LAI動態(tài)模型獲得LAIt+1來描述LAI隨時間的動態(tài)變化，LAI動態(tài)模型是由Dickinson等[29]提出的半經(jīng)驗(yàn)?zāi)Ｐ停Ｐ偷贸鼋Y(jié)果用于模擬樹冠反射率并參與多尺度LAI同化[19，30-31]。

LAIt+1=LAIt+∫t+1tdLdtdt－Lt·LAIt。（1）

dLdt=λ0·R（L′AIt）·L0·（1－exp（－c·LAIt））。（2）

式中：LAIt和LAIt+1表示t和t+1時刻的葉面積指數(shù)；R（L′AIt）是平滑函數(shù)；L′AIt表示LAI歸一化值；L0為最大LAI；Lt為葉片凋落率；參數(shù)c取值為0.5；參數(shù)L0、λ0和Lt由多次擬合得出。

2.1.2PROSASIL模型

PROSAIL模型是先進(jìn)的冠層反射率模型，該模型將PROSPECT模型與SAIL模型相結(jié)合。模型能精確地模擬冠層反射率[19-20，30-31]。

（ρ，τ）=PROSAIL（N，Cab，Car，Cm，Cw）。（3）

ρc=4SAIL（LAI，ρ，τ，H，Psoil，θV，θt，φ）。（4）

式中：ρ為葉片反射率；τ為葉片透光率；ρc為模擬冠層反射率；4SAIL代表4尺度模型，其他參數(shù)見表1。

2.1.3構(gòu)建貝葉斯同化方法

貝葉斯同化方法包括數(shù)據(jù)模型、過程模型和參數(shù)模型組成。利用條件獨(dú)立、思想復(fù)雜的問題化解為若干個相對簡單的模型，并通過條件概率將其聯(lián)系起來。復(fù)雜的聯(lián)合概率推理問題也會被轉(zhuǎn)換成相對簡單的后驗(yàn)概率推理問題來解決。具體如下。

p（過程，參數(shù)|數(shù)據(jù)）∝ 。（5）

（數(shù)據(jù)模型）：p（數(shù)據(jù)|過程，參數(shù)）? 。（5a）

（過程模型）：p（過程|參數(shù)） 。（5b）

（參數(shù)模型）：p（參數(shù)）? 。（5c）

式中，p[a|b]表示已知b的條件下a的條件概率分布。

數(shù)據(jù)模型定義了數(shù)據(jù)和狀態(tài)變量之間的條件概率模型。在研究中有3個量表定義數(shù)據(jù)模型所需的比例，第一尺度為500 m分辨率的MODIS_SG_LAI；第二尺度為100 m分辨率的過渡尺度，第三尺度為20 m分辨率的sentinel-2反射率。2個相鄰的尺度之間為5倍的關(guān)系，一個父節(jié)點(diǎn)對應(yīng)25個子節(jié)點(diǎn)。過程模型定義了真實(shí)過程在3個尺度上的條件依賴關(guān)系，從而得到所有尺度狀態(tài)的后驗(yàn)概率。參數(shù)模型定義了所有參數(shù)的先驗(yàn)分布，采用分辨率相關(guān)似然推斷（ resolution-specific likelihood inference，RESL）和分辨率相關(guān)限制似然推斷（ resolution-specific restricted-likelihood inference，RESREL）實(shí)現(xiàn)參數(shù)極大似然推斷。

在貝葉斯的計(jì)算和推理過程中，涉及到由下而上的過濾和由上而下平滑2個基本過程。向上濾波是將高分辨率的數(shù)據(jù)信息傳遞給低分辨率數(shù)據(jù)，得到節(jié)點(diǎn)在不同尺度下的概率分布；向下平滑是將低分辨率的數(shù)據(jù)信息傳遞給高分辨率數(shù)據(jù)，在3個尺度上更新所有節(jié)點(diǎn)的后驗(yàn)概率分布。

2.1.4向上濾波

通過初始化每一層的初始概率分布，根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)參數(shù)和概率分布函數(shù)，計(jì)算從下到上任一層節(jié)點(diǎn)的后驗(yàn)概率分布，從而將較高分辨率的數(shù)據(jù)信息（底層）傳遞給較低分辨率數(shù)據(jù)（頂層），實(shí)現(xiàn)預(yù)測精細(xì)尺度的信息和更新其他尺度節(jié)點(diǎn)的概率分布。

頂層數(shù)據(jù)初始概率分布和所有層、所有節(jié)點(diǎn)的概率分布均服從正態(tài)分布，見式（6）—式（9）。

p（Ypa（i，r）|Ych（i，r））∝－12（Ypa（i，r）－b（i，r）Ych（i，r））′U－1（i，r）（Ypa（i，r）－b（i，r）Ych（i，r））。（6）

b（i，r）=U（i，r）（1′W－1ch（i，r））。（7）

U（i，r）=（1′W－1ch（i，r）1+W－1ch（i，r））－1。（8）

Wch（i，r）=Qi，rW（i，r）Q′（i，r）。（9）

式中：Ypa（i，r）表示定標(biāo)器中節(jié)點(diǎn)i 觀測數(shù)據(jù)的真值；Ych（i，r）表示對應(yīng)于節(jié)點(diǎn)Ych（i，r）的子節(jié)點(diǎn)（i ，j）的真實(shí)值；式（8）中的1表示所有元素都為1的列向量，-1為逆矩陣，1′表示轉(zhuǎn)置；W（i，r）從參數(shù)模型中獲取；Q′（i，r）是一個正交矩陣。

2.1.5向下平滑

將頂層的后驗(yàn)概率分布作為計(jì)算的初始概率分布，將較低分辨率的數(shù)據(jù)信息傳遞給較高分辨率數(shù)據(jù)，即將所有信息從上往下傳遞，直至到達(dá)最底層，以更新各個尺度所有節(jié)點(diǎn)的后驗(yàn)概率分布。

將自下而上濾波得到最頂層的概率分布作為計(jì)算的初始概率分布，然后計(jì)算任意一點(diǎn)的概率分布

p（Y（i，r）|z（I，0））=∫p（Y（I，R）|Ypa（i，r），Z（i，r））

p（Ypa（i，0））dYpa（ir）。（10）

式中：Z（i，r）是節(jié)點(diǎn)（i，j）的觀測數(shù)據(jù)[30-31]。具體細(xì)節(jié)請參考文獻(xiàn)[32-34]。

2.2.1模型輸入?yún)?shù)的確定

模型中所輸入生理參數(shù)的確定是基于分析文獻(xiàn)查找以及實(shí)測數(shù)據(jù)相結(jié)合，具體參數(shù)見表2。

2.3模型驗(yàn)證方法

固定樣地實(shí)測數(shù)據(jù)為2020年國家森林資源連續(xù)清查數(shù)據(jù)以及2021年野外樣地調(diào)查數(shù)據(jù)，共計(jì)67塊樣地。森林NPP主要包括生物量增長量（ΔBtree）、森林枯損量（ΔMtree）、葉片及細(xì)根更新量（ΔLtree）和林下植被（灌木和草本）NPP貢獻(xiàn)量（ΔNPPund）4部分。其中，生物量增長量為干、枝、葉和根生物量增長量之和，森林枯損量依據(jù)不同森林類型在不同齡組下的平均蓄積枯損率（Mortality Rate ，MR％）估算得出，葉片及細(xì)根更新量計(jì)算公式如下。

Lf=Bf×Tf×Cf 。 （11）

Lfr=Lf×e? 。 （12）

式中：Lf為每年掉落物量；Bf為葉生物量；Tf表示葉周轉(zhuǎn)率；Cf表示葉片碳含量與生物量的轉(zhuǎn)換系數(shù)；Lfr為細(xì)根處的更新量；e表示分配到新細(xì)根的碳含量與新葉碳含量的比例。

2.4變化趨勢分析

采用一元線性回歸分析對本研究得出結(jié)果進(jìn)行分析。變化趨勢分析能更好地反映出塞罕壩植被NPP（式中記為NPP）在近十年的變化趨勢。計(jì)算公式如下[35-36]。

k=n×∑nj=1×NPPj－∑nj=1j∑nj=1NPPjn∑nj=1j2－（∑nj=1j）2。（13）

式中：n為研究時間段的年數(shù)；k為趨勢線的斜率；NPPj為對應(yīng)像元第j年的NPP值；k>0表示NPP呈增加的變化趨勢，k<0則是減少。

2.5植被NPP與氣因子相關(guān)性分析

為了探究NPP對氣象因子的響應(yīng)采用基于像元的空間分析法[37]，本研究中著重計(jì)算了植被NPP對降水與溫度的響應(yīng)。相關(guān)系數(shù)計(jì)算公式如下[35]。

Rab=∑ni=1［（ai－a-）（bi－b-）∑ni=1（ai－a-）2∑ni=1（bi－b-）2。（14）

式中：n為年數(shù)；Rab表示相關(guān)系數(shù)；ai 為第i年NPP；a-為NPP的多年平均值；bi為第i年的溫度或降水量；b-為溫度或降水量的多年平均值[37]。

偏相關(guān)系數(shù)剔除第3個變量影響之后計(jì)算另外2個變量之間的相關(guān)系數(shù)，能更好地反映單一氣候因子對 NPP 的影響，因此本研究的相關(guān)性分析均采用偏相關(guān)系數(shù)。計(jì)算公式如下[37]。

r123=r12－r13r23（1－r213）+（1－r223）。（15）

式中：r123為將變量3固定后變量1與變量2之間的偏相關(guān)系數(shù)；r12、r23、r13分別表示變量1與變量2、變量2與變量3、變量1與變量3的相關(guān)系數(shù)。

3結(jié)果與分析

3.1LAI同化

由圖3可知，MODIS_LAI在生長季LAI值為0～4，且頻繁出現(xiàn)異常值，結(jié)果存在較大誤差。在春季，HBN_LAI_20m LAI值由0.25緩慢增長至1.85，夏季達(dá)到全年最大值3.66，在秋季LAI值由2.63下降至0.77，冬季出現(xiàn)全年最低值0.18，HBN_LAI_20 m符合LAI實(shí)際趨勢。綜上所述，將基于動態(tài)模型模擬的LAI輸入PROSAIL模型對冠層反射率進(jìn)行模擬，并結(jié)合Sentinel-2反射率數(shù)據(jù)，利用層狀貝葉斯網(wǎng)絡(luò)算法得到同化LAI，可以獲得高分辨率、低誤差的MODIS_LAI。

3.2BEPS模型精度驗(yàn)證

對BEPS模型模擬NPP進(jìn)行精度驗(yàn)證，結(jié)果見表3。模型模擬得出的結(jié)果與樣地估算得出結(jié)果較相近，各樣地BEPS模型模擬NPP精度較高（R2=0.77），樣地實(shí)測數(shù)據(jù)估算的NPP與模型模擬得出的NPP進(jìn)行比較，如圖4所示，結(jié)果一致。

3.3植被NPP時空變化分析

研究區(qū)十年間植被NPP均值空間分布如圖5所示。在空間上，由土地覆蓋類型圖結(jié)合圖5可知，塞罕壩機(jī)械林場區(qū)域由西到中，植被NPP值呈現(xiàn)依次降低的情況，這是由于該區(qū)域受到氣候、土壤質(zhì)地及土地覆蓋類型的影響所呈現(xiàn)的獨(dú)特分布格局?；贐EPS模型對2011—2021年進(jìn)行模擬，結(jié)果表明，塞罕壩地區(qū)植被NPP值在0.77～779.0 g/（m2·a）。森林NPP平均值為307.4 g/（m2·a），模擬得出各植被類型NPP值如下。針葉林均值最高，為484.983 8 g/（m2·a）；其次是闊葉林為402.4 g/（m2·a）；草地低于森林均值為288.991 6 g/（m2·a）；混交林為287.9 g/（m2·a）；農(nóng)田耕地為246 g/（m2·a）；灌叢最低僅為145.4821 g/（m2·a），即植被通過光合作用吸收碳的能力由大到小為針葉林、闊葉林、草地、混交林、耕地、灌木。

3.4植被NPP時間序列變化特征

利用BEPS模型模擬得出的區(qū)域NPP采用一元線性回歸分析法進(jìn)行了分析，變化趨勢如圖6所示，在2011—2021年研究區(qū)植被年NPP值在0.77～779.0 g/（m2·a）浮動變化，將一元線性回歸得出結(jié)果按變化斜率進(jìn)行分級研究，定義了8個變化區(qū)間，統(tǒng)計(jì)了各區(qū)間的面積所占比，見表3。由表3可知，塞罕壩地區(qū)NPP整體呈穩(wěn)步增長趨勢；塞罕壩中部闊葉林為主要上升區(qū)，下降區(qū)面積則集中在部分農(nóng)田和草原地帶。

3.5植被NPP對氣候因子的相關(guān)性分析

以年為單位，對2011—2021年的植被NPP與降水量、均溫的相關(guān)系數(shù)進(jìn)行分析，如圖7所示。由圖7（a）可知，塞罕壩機(jī)械林場植被NPP與溫度變化趨勢存在顯著相關(guān)。塞罕壩植被NPP整體上與溫度相關(guān)性較高，相關(guān)系數(shù)為0.2～0.8，因塞罕壩地區(qū)氣候寒冷，冬季時間長，春秋時間短而夏季又不明顯。年均氣溫在-1.2 ℃，一年中冬季時間占全年2/3，所以在該地區(qū)溫度成為了影響植被NPP的主要影響因子。而植被NPP與降水量的相關(guān)性總體而言沒有與溫度的相關(guān)性高，在該地區(qū)植被NPP與降水量的相關(guān)性在-0.3～0.4，由此可知在該地區(qū)降水量對植被NPP的影響較低，而溫度成為該地區(qū)主導(dǎo)因子。

4結(jié)論與討論

4.1討論

本研究采用LAI動態(tài)模型和PROSAIL模型與層狀貝葉斯網(wǎng)絡(luò)算法相耦合的方法構(gòu)建森林LAI高精度時空同化系統(tǒng)，實(shí)現(xiàn)了在20 m尺度上森林LAI時間序列數(shù)據(jù)高精度同化，并以此來驅(qū)動BEPS模型，將模型分辨率提高到20 m，這為研究相對較小尺度區(qū)域的森林生產(chǎn)力提供了準(zhǔn)確依據(jù)。與此同時，本研究模擬了塞罕壩地區(qū)近十年來NPP變化規(guī)律；獲取的高空間分辨率LAI數(shù)據(jù)，極大地提高了LAI數(shù)據(jù)的精度，結(jié)果與JI等[33]研究一致。但本研究未獲取到充足的LAI實(shí)測數(shù)據(jù)樣本，因此在精度驗(yàn)證部分僅與MODIS LAI進(jìn)行比對，在后續(xù)研究中可以補(bǔ)充實(shí)測數(shù)據(jù)，以完善精度評價(jià)。本研究利用2020年二類調(diào)查數(shù)據(jù)和2021年外業(yè)實(shí)測數(shù)據(jù)計(jì)算NPP真實(shí)值，塞罕壩機(jī)械林場為人工林，所以林下植被較少，因此，將灌木和草本忽略不計(jì)，只計(jì)算了活立木的生物量。在氣象因子選擇方面，只選取了溫度與降水量2個主要?dú)庀笠蜃舆M(jìn)行分析研究，在后續(xù)研究中可以選擇加入其他氣象因子，進(jìn)而更加全面地分析氣象因素對該地區(qū)NPP的影響。

4.2結(jié)論

本研究以塞罕壩機(jī)械林場為研究區(qū)域，以2021年MODIS LAI時間序列產(chǎn)品和Sentinel-2反射率數(shù)據(jù)為數(shù)據(jù)源，基于LAI動態(tài)模型、PROSAIL輻射傳輸模型和層狀貝葉斯網(wǎng)絡(luò)構(gòu)建數(shù)據(jù)同化系統(tǒng)，獲得LAI高時空分辨率數(shù)據(jù)，然后將同化的高時空分辨率LAI作為輸入數(shù)據(jù)驅(qū)動BEPS模型，實(shí)現(xiàn)對塞罕壩林場植被凈初級生產(chǎn)力時空變化及影響因素的精準(zhǔn)模擬與分析，研究結(jié)論包括以下4個方面。

1）基于層狀貝葉斯網(wǎng)絡(luò)算法得到的同化LAI，可以獲得高分辨率、低誤差的LAI數(shù)據(jù)。在本研究區(qū)域內(nèi)，春季，HBN_LAI_20mLAI值由0.25緩慢增長至1.85，夏季達(dá)到全年最大值3.66，在秋季LAI值由2.63下降至0.77，冬季出現(xiàn)全年最低值0.18，HBN_LAI_20m符合LAI實(shí)際趨勢。

2）各樣地BEPS模型模擬NPP精度較高（R2=0.77），塞罕壩地區(qū)森林NPP值在0.77～779.0 g/（m2·a）范圍內(nèi)，森林NPP平均值為307.4 g/（m2·a）。 2010—2021年，塞罕壩地區(qū)植被NPP整體呈上升趨勢，塞罕壩中部闊葉林為主要上升區(qū)，下降區(qū)面積則集中在部分農(nóng)田、草原地帶。

3）塞罕壩植被NPP整體上與溫度相關(guān)性較高，相關(guān)系數(shù)在0.2～0.8，與降水的相關(guān)性相對較低，相關(guān)系數(shù)在-0.3～0.4，因此該地區(qū)主導(dǎo)因子為溫度因子。

4）本研究獲知塞罕壩地區(qū)近十年來的植樹造林和積極推動的各項(xiàng)林業(yè)政策、保護(hù)措施在促進(jìn)森林生態(tài)建設(shè)過程中發(fā)揮了重要作用，提升了塞罕壩地區(qū)植被NPP。

【參考文獻(xiàn)】

[1]蘇尚軍，靳東升，張婷.山西省植被凈初級生產(chǎn)力時空變化特征及影響因素[J].山西農(nóng)業(yè)科學(xué)，2022，50（4）：551-558.

SU S J， JIN D S， ZHANG T. Characteristics temporal and spatial change and impact factors of net primary productivity of vegetation in Shanxi Province[J]. Journal of Shanxi Agricultural Sciences， 2022， 50（4）： 551-558.

[2]王金杰，趙安周，胡小楓.京津冀植被凈初級生產(chǎn)力時空分布及自然驅(qū)動因子分析[J].生態(tài)環(huán)境學(xué)報(bào)，2021，30（6）：1158-1167

WANG J J， ZHAO A Z， HU X F. Spatiotemporal distribution of vegetation net primary productivity in Beijing-Tianjin-Hebei and natural driving factors[J]. Ecology and Environment Sciences， 2021， 30（6）： 1158-1167.

[3]LIETH H， WHITTAKER R H. Primary productivity of the biosphere[M]. New York： Springer-Verlag， 1975.

[4]趙穎慧，李思琪，甄貞，等.黑龍江省黑河地區(qū)森林地上生物量和NPP估測及時空格局[J].應(yīng)用生態(tài)學(xué)報(bào)，2016，27（10）：3070-3080

ZHAO Y H， LI S Q， ZHEN Z， et al. Estimation and spatial-temporal pattern of forest aboveground biomass and NPP in Heihe Region， Heilongjiang Province， China[J]. Journal of Applied Ecology， 2016， 27（10）： 3070-3080.

[5]潘竟虎，黃克軍，李真.2001—2010年疏勒河流域植被凈初級生產(chǎn)力時空變化及其與氣候因子的關(guān)系[J].生態(tài)學(xué)報(bào)，2017，37（6）：1888-1899

PAN J H， HUANG K J， LI Z. Spatio -temporal variation in vegetation net primary productivity and its relationship with climatic factors in the Shule River basin from 2001 to 2010[J]. Acta Ecologica Sinica， 2017， 37（6）： 1888-1899.

[6]LIU L， GUAN J Y， HAN W Q， et al. Quantitative assessment of the relative contributions of climate and human factors to net primary productivity in the Ili river basin of China and Kazakhstan[J]. Chinese Geographical Science， 2022， 32（6）： 1069-1082.

[7]潘天石，劉哲，范文義.BEPS模型中光合模塊機(jī)理性參數(shù)的適用性[J].東北林業(yè)大學(xué)學(xué)報(bào)，2018，46（5）：69-73

PAN T S， LIU Z， FAN W Y. Mechanism parameter applicability of photosynthetic modules in BEPS model[J]. Journal of Northeast Forestry University， 2018， 46（5）： 69-73.

[8]FENG X F， LIU G H， ZHOU W Z， et al. Net primary productivity distribution in China from a process model driven by remote sensing[C]. 2005 IEEE International Geoscience and Remote Sensing Symposium （IGARSS 2005）， Seoul， Korea， 2005.

[9]王培娟，孫睿，朱啟疆，等.復(fù)雜地形條件下提高BEPS模型模擬能力的途徑[J].中國圖象圖形學(xué)報(bào)，2006（7）：1017-1025，1044.

WANG P J， SUN R， ZHU Q J， et al. Improvement on the abilities of BEPS under accidented terrain[J]. Journal of Image and Graphics， 2006（7）： 1017-1025， 1044.

[10]MO X G， CHEN J M， JU W M， et al. Optimization of ecosystem model parameters through assimilating eddy covariance flux data with an ensemble Kalman filter[J]. Ecological Modelling， 2008， 217（1/2）： 157-173.

[11]毛學(xué)剛，焦裕欣，張穎.基于BEPS模型的東北三省森林生態(tài)系統(tǒng)NPP模擬[J].森林工程，2017，33（1）：22-27.

MAO X G， JIAO Y X， ZAHNG Y. NPP simulation of forest ecosystem in three provinces of northeast China based on BEPS model[J]. Forest Engineering， 2017， 33（1）： 22-27.

[12]毛學(xué)剛.東北森林碳循環(huán)日步長模型與遙感綜合應(yīng)用研究[D].哈爾濱：東北林業(yè)大學(xué)，2011.

MAO X G. Study on the model of Northeast forest carbon cycle at daily step and the integrated application of remote sensing[D]. Harbin： Northeast Forestry University， 2011.

[13]薛梅，陳蕓芝，閆敏，等.中國東北地區(qū)植被NPP模擬與時空變化分析[J].福州大學(xué)學(xué)報(bào)（自然科學(xué)版），2018，46（6）：821-830.

XUE M， CHEN Y Z， YAN M， et al. Simulation and spatio-temporal variation analysis of net primary productivity in Northeast China[J]. Journal of Fuzhou University （Natural Science Edition）， 2018， 46（6）： 821-830.

[14]XU C， LIU M， AN S， et al. Assessing the impact of urbanization on regional net primary productivity in Jiangyin County， China[J]. Journal of Environmental Management， 2007， 85（3）： 597-606.

[15]ZHANG T L， SUN R， ZHANG R H， et al. Simulation of water and carbon fluxes in Harvard forest area based on data assimilation method[J]. The journal of applied ecology， 2013， 24（10）： 2746-2754.

[16]JONCKHEERE I， FLECK S， NACKAERTS K， et al. Review of methods for in situ leaf area index determination： Part I. Theories， sensors and hemispherical photography[J]. Agricultural & Forest Meteorology， 2004， 121（1/2）： 19-35.

[17]MA H， SONG J L， WANG J D， et al. Improvement of spatially continuous forest LAI retrieval by integration of discrete airborne LiDAR and remote sensing multi-angle optical data[J]. Agricultural and Forest Meteorology， 2014：189-190.

[18] LIU Z L， CHEN J M， JIN G Z， et al. Estimating seasonal variations of leaf area index using litterfall collection and optical methods in four mixed evergreen-deciduous forests[J]. Agricultural and Forest Meteorology， 2015：209-210.

[19]李雪建，毛方杰，杜華強(qiáng)，等.雙集合卡爾曼濾波LAI同化結(jié)合BEPS模型的竹林生態(tài)系統(tǒng)碳通量模擬[J].應(yīng)用生態(tài)學(xué)報(bào)，2016，27（12）：3797-3806

LI X J， MAO F J， DU H Q， et al. Simulating of carbon fluxes in bamboo forest ecosystem using BEPS model based on the LAI assimilated with Dual Ensemble Kalman Filter[J]. Journal of Applied Ecology， 2016， 27（12）： 3797-3806.

[20]XING L Q， LI X J， DU H Q， et al. Assimilating multiresolution leaf area index of moso bamboo forest from MODIS time series data based on a hierarchical bayesian network algorithm[J]. Remote Sensing， 2018， 11（1）：56.

[21]張時煌，彭公炳，黃玫.基于地理信息系統(tǒng)技術(shù)的土壤質(zhì)地分類特征提取與數(shù)據(jù)融合[J].氣候與環(huán)境研究，2004，9（1）：65-78.

ZHANG S H， PENG G B， HUANG M. The feature extraction and data fusion of regional soil textures based on GIS techniques[J]. Climatic and Environmental Research， 2004， 9（1）： 65-78.

[22]BRUAND A， BAJZE D， HARDY M. Prediction of water retention properties of clay soil validity of relationships using a single soil characteristic [J]. Soil Use and Management， 1994， 10（3）： 99-103.

[23]RUNNING S W， COUGHLAN J C. A general model of forest ecosystem processes for regional applications I. Hydrologic balance， canopy gas exchange and primary production processes[J]. Ecological Modelling， 1988， 42（2）： 125-154.

[24]LIU J， CHEN J M， CIHLAR J， et al. A process-based boreal ecosystem productivity simulator using remote sensing inputs[J]. Remote Sensing of Environment， 1997， 62（2）： 158-175.

[25]LIU J， CHEN J M， CIHLAR J， et al. Net primary productivity distribution in the BOREAS region from a process model using satellite and surface data[J]. Journal of Geophysical Research， 1999， 104（D22）：27735-27754.

[26]LIU Q， GU L， DICKINSON R E， et al. Seasonally varying leaf area for climate and carbon models from assimilation of satellite reflectance data into a dynamical leaf model[C]. Agu Fall Meeting Abstracts， 2007.

[27]LI H， CHEN Z X， WU W B， et al. Crop model data assimilation with particle filter for yield prediction using leaf area index of different temporal scales[C]// 2015 the Fourth International Conference on Agro-Geoinformatics， IEEE， 10 September 2015， Istanbul， Turkey， 2015.

[28]李雪建.竹林MODISLAI時間序列同化及在碳通量模擬中的應(yīng)用研究[D].杭州：浙江農(nóng)林大學(xué)，2017.

LI X J. Assimilation of MODIS LAI time series in bamboo forest and its application in carbon flux simulation[D]. Hangzhou： Zhejiang Agriculture & Forestry University， 2017.

[29]DICKINSON R E， TIAN Y， LIU Q， et al. Dynamics of leaf area for climate and weather models[J]. Journal of Geophysical Research Atmospheres， 2008， 113（D16）： D16115.

[30]HARIKA， NARUMANCHI， DISHANT， et al. Estimating bamboo forest aboveground biomass using EnKF-assimilated MODIS LAI spatiotemporal data and machine learning algorithms[C]//2018 IEEE International Conference on Smart Cloud， 2018.

[31]MAO F J， LI X J， DU H Q， et al. Comparison of two data assimilation methods for improving MODIS LAI time series for bamboo forests[J]. Remote Sensing， 2017， 9（5）： 401.

[32]邢璐琪.基于多源遙感數(shù)據(jù)的竹林LAI多尺度同化及在碳循環(huán)模擬中的應(yīng)用[D].杭州：浙江農(nóng)林大學(xué)，2019.

XING L Q. Assimilation of multiresolution lai in bamboo forest and its application in carbon cycle simulation based on multi-source remote sensing data[D]. Hangzhou： Zhejiang Agriculture & Forestry University， 2019.

[33]JI J Y， LI X J ， DU H Q， et al. Multiscale leaf area index assimilation for Moso bamboo forest based on Sentinel-2 and MODIS data[J]. International Journal of Applied Earth Observation and Geoinformation， 2021， 104（C4）： 102519.

[34]JIANG J Y， XIAO Z Q， WANG J D， et al. Multiscale estimation of leaf area index from satellite observations based on an ensemble multiscale filter[J]. Remote Sensing， 2016， 8（3）： 229-229.

[35]BERLINER L M. Hierarchical Bayesian space-time analysis[M]. New York： John Wiley & Sons， Ltd， 2006.

[36]ANG D D， LIANG S L. Using multiresolution tree to integrate MODIS and MISR-L3 LAI products[C]//IEEE International Geoscience and Remote Sensing Symposium， 2010， 1027-1030.

[37]穆少杰，李建龍，周偉，等.2001—2010年內(nèi)蒙古植被凈初級生產(chǎn)力的時空格局及其與氣候的關(guān)系[J].生態(tài)學(xué)報(bào)，2013，33（12）：3752-3764.

MU S J， LI J L， ZHOU W， et al. Spatial-temporal distribution of net primary productivity and its relationship with climate factors in Inner Mongolia from 2001 to 2010[J]， Acta Ecologica Sinica， 2013， 33（12）： 3752-3764.

[38]劉文瑞，李曉婷，李彤，等.基于MODIS和CASA模型的伊春市森林植被NPP變化特征及其影響因子分析[J].生態(tài)學(xué)雜志，2022，41（1）：150-158.

LIU W R， LI X T， LI T， et al. Spatiotemporal variations of forest NPP and related driving factors based on MODIS and CASA models in Yichun[J]， Chinese Journal of Ecology， 2022， 41（1）： 150-158.

[39]穆少杰，李建龍，楊紅飛，等.內(nèi)蒙古草地生態(tài)系統(tǒng)近10年NPP時空變化及其與氣候的關(guān)系[J].草業(yè)學(xué)報(bào)，2013，22（3）：6-15.

MU S J， LI J L， YANG H F， et al. Patio-temporal variation analysis of grassland net primary productivity and its relationship with climate over the past 10 years in Inner Mongolia[J]. Acta Prataculturae Sinica， 2013， 22（3）： 6-15.