包志意，范文義

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

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

中圖分類號：S718.55文獻標識碼：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)植被利用光合作用所固定的有機物總量與自養(yǎng)呼吸消耗的差值，通過植被與外界環(huán)境因子相互作用下所產(chǎn)生[1]。NPP是反映植被的生產(chǎn)能力和生態(tài)過程的指標，用以評價生態(tài)系統(tǒng)植被質(zhì)量狀況和區(qū)域生態(tài)系統(tǒng)功能[2-3]。

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

在以往的森林生產(chǎn)力模擬研究中，由于用于驅(qū)動模型的LAI數(shù)據(jù)缺乏可信度，導致模型模擬結果存在較大誤差[14]。葉面積指數(shù)（Leaf Area Index，LAI）定義為單位地表面積綠葉總面積的一半，對植物生理過程以及生態(tài)系統(tǒng)生產(chǎn)力的形成有著重要影響[15-16]，是碳循環(huán)模擬的重要驅(qū)動數(shù)據(jù)[17-18]，獲得高分辨率時間序列LAI數(shù)據(jù)是準確模擬區(qū)域碳循環(huán)的關鍵。但從單一時相遙感數(shù)據(jù)中反演LAI，難以獲取長時間序列LAI[17]，MODIS衛(wèi)星遙感提供LAI產(chǎn)品又存在噪聲大、在時間序列上波動較大等缺陷。因此，獲取高精度的LAI數(shù)據(jù)對于區(qū)域生態(tài)系統(tǒng)生產(chǎn)力的精確模擬變得尤為重要。越來越多的學者采用數(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模型的一個關鍵輸入?yún)?shù)，在以往研究中大多采用MODIS LAI產(chǎn)品，其存在空間分辨率低、混合像元較多以及區(qū)域應用中出現(xiàn)明顯高估或低估現(xiàn)象等問題，如何獲取可信度和時空分辨率更高的LAI數(shù)據(jù)變得尤為重要。因此，本研究以塞罕壩機械林場為研究對象，以2021年MODIS LAI時間序列產(chǎn)品和Sentinel-2反射率數(shù)據(jù)為數(shù)據(jù)源，基于LAI動態(tài)模型、PROSAIL輻射傳輸模型和層狀貝葉斯網(wǎng)絡（Hierarchical Bayesian Network， HBN）構建數(shù)據(jù)同化系統(tǒng)，獲得LAI高時空分辨率數(shù)據(jù)，然后將同化的高時空分辨率LAI作為輸入數(shù)據(jù)驅(qū)動BEPS模型，實現(xiàn)對塞罕壩地區(qū)植被凈初級生產(chǎn)力時空變化及影響因素的精準模擬與分析，為植被生態(tài)系統(tǒng)碳循環(huán)的時空模擬提供更精確的葉面積指數(shù)數(shù)據(jù)，為今后塞罕壩森林合理經(jīng)營管理、科學規(guī)劃提供技術支撐。

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

1.1研究區(qū)概況

塞罕壩機械林場位于河北省承德市圍場滿族蒙古族自治縣最北部區(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ù)點。在NASA網(wǎng)站（https：//ladsweb.modaps.eosdis.nasa.gov/）下載了2021年全年的MODIS LAI產(chǎn)品。通過使用MODIS重投影工具（MRT）軟件，將其重投影至WGS-84坐標系，空間分辨率為500 m×500 m。使用ENVI 5.3對研究區(qū)域進行裁剪，并提取了LAI值。最后，使用Savitzky-Golay（SG）濾波算法，以減少數(shù)據(jù)噪聲和消除異常。SG平滑后的LAI（MODIS_SG_LAI）被用作貝葉斯層狀網(wǎng)絡（HBN）的初始輸入數(shù)據(jù)。

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

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

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

本研究所使用的氣象數(shù)據(jù)為逐日降水量、最高溫度、最低溫度、太陽輻射以及相對濕度等。以上數(shù)據(jù)由“國家青藏高原科學數(shù)據(jù)中心”下載。對氣象數(shù)據(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ù)中國土壤類型分類編碼與美國土壤分類系統(tǒng)編碼的對應關系轉(zhuǎn)換為美國土壤分類系統(tǒng)分類標準的土壤質(zhì)地矢量圖，其次，將轉(zhuǎn)化后的土壤質(zhì)地類型圖根據(jù)土壤有效持水量與土壤質(zhì)地間的聯(lián)系把土壤質(zhì)地圖轉(zhuǎn)換為土壤有效持水量圖，最后將土壤有效持水量矢量數(shù)據(jù)轉(zhuǎn)換成空間分辨率為1 km 柵格數(shù)據(jù)。

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

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

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

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

2研究方法

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

研究中LAI數(shù)據(jù)的同化是以空間分辨率為500 m的MODIS LAI 和空間分辨率為20 m的sentinel-2反射率數(shù)據(jù)為基礎數(shù)據(jù)，由LAI動態(tài)模型，PROSAIL模型與層狀貝葉斯網(wǎng)絡算法組成的同化系統(tǒng)實現(xiàn)高空間分辨率LAI數(shù)據(jù)的同化。然后將同化后的LAI數(shù)據(jù)結合氣象數(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)驗模型，模型得出結果用于模擬樹冠反射率并參與多尺度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模型是先進的冠層反射率模型，該模型將PROSPECT模型與SAIL模型相結合。模型能精確地模擬冠層反射率[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構建貝葉斯同化方法

貝葉斯同化方法包括數(shù)據(jù)模型、過程模型和參數(shù)模型組成。利用條件獨立、思想復雜的問題化解為若干個相對簡單的模型，并通過條件概率將其聯(lián)系起來。復雜的聯(lián)合概率推理問題也會被轉(zhuǎ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倍的關系，一個父節(jié)點對應25個子節(jié)點。過程模型定義了真實過程在3個尺度上的條件依賴關系，從而得到所有尺度狀態(tài)的后驗概率。參數(shù)模型定義了所有參數(shù)的先驗分布，采用分辨率相關似然推斷（ resolution-specific likelihood inference，RESL）和分辨率相關限制似然推斷（ resolution-specific restricted-likelihood inference，RESREL）實現(xiàn)參數(shù)極大似然推斷。

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

2.1.4向上濾波

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

頂層數(shù)據(jù)初始概率分布和所有層、所有節(jié)點的概率分布均服從正態(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）表示定標器中節(jié)點i 觀測數(shù)據(jù)的真值；Ych（i，r）表示對應于節(jié)點Ych（i，r）的子節(jié)點（i ，j）的真實值；式（8）中的1表示所有元素都為1的列向量，-1為逆矩陣，1′表示轉(zhuǎn)置；W（i，r）從參數(shù)模型中獲取；Q′（i，r）是一個正交矩陣。

2.1.5向下平滑

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

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

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é)點（i，j）的觀測數(shù)據(jù)[30-31]。具體細節(jié)請參考文獻[32-34]。

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

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

2.3模型驗證方法

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

Lf=Bf×Tf×Cf 。 （11）

Lfr=Lf×e? 。 （12）

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

2.4變化趨勢分析

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

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

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

2.5植被NPP與氣因子相關性分析

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

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

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

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

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

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

3結果與分析

3.1LAI同化

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

3.2BEPS模型精度驗證

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

3.3植被NPP時空變化分析

研究區(qū)十年間植被NPP均值空間分布如圖5所示。在空間上，由土地覆蓋類型圖結合圖5可知，塞罕壩機械林場區(qū)域由西到中，植被NPP值呈現(xiàn)依次降低的情況，這是由于該區(qū)域受到氣候、土壤質(zhì)地及土地覆蓋類型的影響所呈現(xiàn)的獨特分布格局?；贐EPS模型對2011—2021年進行模擬，結果表明，塞罕壩地區(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采用一元線性回歸分析法進行了分析，變化趨勢如圖6所示，在2011—2021年研究區(qū)植被年NPP值在0.77～779.0 g/（m2·a）浮動變化，將一元線性回歸得出結果按變化斜率進行分級研究，定義了8個變化區(qū)間，統(tǒng)計了各區(qū)間的面積所占比，見表3。由表3可知，塞罕壩地區(qū)NPP整體呈穩(wěn)步增長趨勢；塞罕壩中部闊葉林為主要上升區(qū)，下降區(qū)面積則集中在部分農(nóng)田和草原地帶。

3.5植被NPP對氣候因子的相關性分析

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

4結論與討論

4.1討論

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

4.2結論

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

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

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整體上與溫度相關性較高，相關系數(shù)在0.2～0.8，與降水的相關性相對較低，相關系數(shù)在-0.3～0.4，因此該地區(qū)主導因子為溫度因子。

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

【參考文獻】

[1]蘇尚軍，靳東升，張婷.山西省植被凈初級生產(chǎn)力時空變化特征及影響因素[J].山西農(nóng)業(yè)科學，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)境學報，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].應用生態(tài)學報，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)力時空變化及其與氣候因子的關系[J].生態(tài)學報，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模型中光合模塊機理性參數(shù)的適用性[J].東北林業(yè)大學學報，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]王培娟，孫睿，朱啟疆，等.復雜地形條件下提高BEPS模型模擬能力的途徑[J].中國圖象圖形學報，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]毛學剛，焦裕欣，張穎.基于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]毛學剛.東北森林碳循環(huán)日步長模型與遙感綜合應用研究[D].哈爾濱：東北林業(yè)大學，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].福州大學學報（自然科學版），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]李雪建，毛方杰，杜華強，等.雙集合卡爾曼濾波LAI同化結合BEPS模型的竹林生態(tài)系統(tǒng)碳通量模擬[J].應用生態(tài)學報，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)技術的土壤質(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時間序列同化及在碳通量模擬中的應用研究[D].杭州：浙江農(nóng)林大學，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)模擬中的應用[D].杭州：浙江農(nóng)林大學，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)力的時空格局及其與氣候的關系[J].生態(tài)學報，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)學雜志，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時空變化及其與氣候的關系[J].草業(yè)學報，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.