華北平原農田生態系統的Noah-MP模型參數敏感性分析

王 茜，查元源

（武漢大學水資源工程與調度全國重點實驗室，武漢 430072）

0 引 言

農田生態系統與人類生命、生活、生產息息相關，它為我們提供賴以生存的糧食，同時也消耗著全球淡水抽取量的70%和總耗水量的90%[1]。作為水-糧食-能源三者相互聯系的復雜系統的縮影，準確刻畫農田生態系統中涉及的水、能量及CO2通量變化過程，為區域水資源配置、糧食生產、生態環境保護提供理論依據，對保障國家水安全、糧食安全、生態安全有著重要意義。陸面模式作為模擬陸面上發生的各種物理、化學、生物和水文過程，以及這些過程與大氣的相互作用的有效工具，歷經簡單模式階段與生物大氣模式階段，進入了新一代多模式階段[2]。從簡單地考慮水量平衡，到將陸氣間水分與熱量的交換過程進行耦合；從不考慮植被的作用，到顯式引入植被的生物物理作用，再到通過耦合葉片的光合作用與蒸騰作用實現陸氣交換過程的水碳耦合，模式的物理機制逐漸完善，對陸面過程與陸氣交換的刻畫也逐漸細化。新一代陸面模式中，Noah-MP[3]模型集合多種參數化方案選項，提供多達4 608 種參數化方案組合，方便使用者依據適用性靈活選擇，具有廣闊的應用前景。

隨著陸面模式的不斷完善，參數不確定性成為了限制模型模擬性能的主要因素。進入大數據時代，衛星、無人機等多種遙感技術快速發展，使得基于天-空-地多源觀測數據建立高分辨率的土壤屬性與植被屬性等地表參數庫成為可能[4]。因此，對陸面模式進行參數敏感性分析，識別出影響模型模擬效果的敏感參數，不僅能夠指導參數優化以提升模型模擬性能，也能為地表參數庫的構建指引方向。

近年來，許多學者展開了Noah-MP 模型的參數敏感性分析研究。在密西西比河流域評估Noah-MP 陸面模式模擬水文變量的表現時，CAI 等[5]簡單地對3 個徑流模擬關鍵參數進行了敏感性測試，并在此基礎上以納什效率系數（NSE）為性能指標進行參數優化，從而提升模型徑流模擬能力。ARSENAULT 等[6]利用10 個國際FluxNet 站點的數據對動態植被參數、靜態植被參數及土壤參數進行了全局敏感性分析，找出了12 個敏感參數，并指出針對特定區域的具體問題需要對參數分布范圍進行進一步約束。而HUO等[7]重點關注草及常綠闊葉林2種典型植被，提出了一種空間抽樣與參數抽樣相結合的新方法，在美國大陸的網格單元中篩選出植被總初級生產力（GPP）和潛熱通量（LE）模擬的敏感參數。SHU 等[8]在黃土高原的研究顯示，相比于直接同化葉面積指數（LAI），基于參數敏感性分析的參數優化能更大幅度地降低葉面積指數的模擬值與觀測值的均方根誤差（RMSE）。上述研究均表明，受不同的下墊面及氣候區的影響，需開展針對性的Noah-MP 模型參數敏感性分析試驗，以助挖掘參數優化對模型模擬性能提升的巨大潛力。目前，已有大量有關Noah-MP 模型在華北平原農田生態系統的應用型研究[9-11]，但在該區域開展的參數敏感性分析研究未見報道。華北平原是確保全國糧食安全的重要基地[12]，基于Noah-MP 模型開展參數敏感性分析，有助于準確量化這一典型農業生態系統模式與大氣間的輻射傳輸、物質和能量交換，對推動區域經濟、社會和生態可持續性協調發展有著促進作用。

1 數據與方法

1.1 觀測站點與數據

中國國家生態系統觀測研究網絡（CERN）的禹城農業綜合試驗站（116.570 2°E，36.829 0°N，平均海拔22.0 m）與欒城農業生態系統試驗站（114.699 3°E，37.893 0°N，平均海拔50.1 m）均位于華北平原，屬暖溫帶半濕潤季風氣候區，以冬小麥-夏玉米為主要作物，一年2 季?？紤]到華北平原輪作特點，即每年10月至次年6月為冬小麥生長季，6月至10月為夏玉米生長季，結合數據質量與可獲取性等原因，本研究以2020年6月10日至2021年6月9日的輪作年為研究時段。獲取的觀測數據包括基本氣象數據、多層土壤水（0.1、0.2、0.3、0.4和0.5 m處）及水熱碳通量。研究時段內的數據存在大量缺失，需要進行插補處理。插補過程參考R 語言的fluxnetLSM[13]程序包，該程序包主要用于插補遵循Fluxnet 2015 協議的數據。處理后的氣象數據仍存在大量連續缺省值難以作為模型所需輸入數據，故不加以利用，而土壤水與通量觀測數據缺失率低于15%，可以用于模型敏感性分析檢驗。

1.2 Noah-MP模型

Noah-MP 模型（The Community Noah Land Surface Model with Multi-Parameterization Options）是在Noah 3.0 版本的基礎上進行改進和完善的新一代陸面模式。Noah 模型將陸面簡單視為0～0.1、0.1～0.4、0.4～1.0 和1.0～2.0 m 的4 層土壤，能量平衡計算在土壤與大氣的交界面展開，Noah-MP 模型在此基礎上調整模型架構，將冠層與地表分開，引入3 層雪層，同時對積雪/融雪、產匯流、植被生長等關鍵過程進行優化，并為表1 所示的12 個陸面過程提供多種參數化方案選項。Noah-MP 模型能很好地模擬地表能量收支平衡及陸地水文循環過程，已被作為陸面參數化方案嵌入GFS、CFS 及WRF 等氣候/天氣預測模式中。為簡化輸入輸出文件讀取過程，修改網格化模型的輸入輸出端口，得到用于站點敏感性分析的單點模型（https://github.com/wangxi1998/NoahMP-Sensitivity-NCP/tree/main/POINT_NOAHMP）。

表1 模型參數化方案設置Tab.1 Model parameterization scheme

1.3 Noah-MP模型參數敏感性分析設計

1.3.1 模型模擬試驗的設計

模型所需的氣象數據來自GLDAS-2.1[14,15]，空間分辨率為0.25°，時間分辨率為3 h。該數據為Noah-MP 等陸面模式的“在線”模擬集合，故常用作Noah-MP 模型的氣象輸入數據。利用GEE 平臺提取2 個站點在研究時段的數據（2 m 氣溫、2 m 比濕、10 m 風速、氣壓、下行短波輻射、下行長波輻射及降水），將其插值得到0.5 h分辨率的數據作為輸入。模型初始場參考模型指導手冊確定，分層土壤溫度與地表溫度取全球均值（274.0、278.0、280.0、284.0 和287.0 K），土壤含水率為飽和含水率的40%（取為0.24 m3/m3）。

選用的參數化方案及參數參考ARSENAULT 等[6]和HUO等[7]的研究（見表1 和表2）。選取的42 個參數包含30 個植被相關參數與12 個土壤相關參數；植被相關參數中有12 個與輻射傳輸過程密切相關，其余參數為動態植被模擬模塊涉及參數。模型模擬過程中，模型將土地利用類型與土壤類型作為索引條件，再在模型自帶的查找表中匹配對應參數。因此，大部分參數的取值上下限通過查找表確定。少量參數間具有一定大小關系，如冠層頂部高度HVT 與冠層底部高度HVB、田間持水量SMCREF 與孔隙度SMCMAX、凋萎系數SMCWLT與田間持水量SMCREF 及蒸發停止的土壤含水率SMCDRY，需要對它們的取值范圍進行約束。

表2 敏感性分析選用的參數及其取值范圍Tab.2 Selected parameters and the value ranges for sensitivity analysis

模型運行步長為0.5 h，并利用研究時段內的氣象數據進行多輪預熱。為了對模型參數敏感性進行系統性的分析，本研究選取潛熱通量（LE）、顯熱通量（H）、凈生態系統碳交換量（NEE）、第1 層土壤含水率（SM1）及第2 層土壤含水率（SM2）進行分析。ARSENAULT 等[6]通過數值試驗發現當樣本量接近1 500 時，全局敏感度基本收斂。更重要的是，那些在小樣本量下表現出較高靈敏度的參數在較大樣本量下也表現出較高靈敏度，反之亦然。為了減少計算成本，先利用ARSENAULT 等[6]研究中的算例對42 個參數抽樣分別抽樣250次、500 次、1 000 次與1 500 次，探究合理的抽樣次數。有研究表明，動態植被模塊開啟時，氣孔阻抗的土壤水分限制因子的參數化方案對模型參數的取值有較大影響。因此，除對2個站點進行參數敏感性分析外，以CLM 和Noah 2 種氣孔阻抗的土壤水分限制因子方案為代表，探索基于土壤水勢與基于土壤含水率的參數化方案對參數敏感性的影響。

1.3.2 參數敏感性分析評估方法

本研究采用Matlab 工具箱中的拉丁超立方抽樣法（LHS）對參數抽樣，抽樣得到的參數代入模型運行后，基于Sobol 方差分析計算參數的敏感性指標（實現代碼獲取網址為https://github.com/wangxi1998/NoahMP-Sensitivity-NCP/tree/main）。全局敏感度Ti（Ti越大代表該模型輸出變量對第i個參數越敏感）的計算方法如下[6,16]：

式中：x 為進行敏感性分析的N 個參數的向量（N=42）；xi為第i個參數，x～i為除xi外的N-1個其他參數的向量。蒙特卡洛抽樣的集合平均計算方法如下：

式中：M 為拉丁超立方抽樣次數；xm為第m 次抽樣過程得到的參數向量；是拉丁超立方抽樣的2個平行抽樣與的參數；f(x)為性能指標。

性能指標選用衡量預測值與真實值之間偏差的均方根誤差，計算方法如下：

式中：K 為模型輸出變量的時間維度；Rk為第k 個模型運行步長時刻的輸出結果；Ok為對應時刻的觀測值。

2 結果與分析

2.1 模型預熱輪數與抽樣次數M的確定

模型預熱以相鄰2輪的平均相對誤差達0.001為結束條件，判斷標準見式(4)。2站點的預熱過程類似，此處僅展示禹城站結果（見圖1），圖1中結果均為各輪預熱后模型的最終結果。

圖1 禹城站預熱過程Fig.1 Spin-up process of Yucheng station

式中：Rk為第k 個模型運行步長時刻的輸出結果；K 為模型輸出變量的時間維度；n為變量達到穩定所需預熱輪數。

不同輸出變量達到穩定所需預熱時間的結果如表3所示。結果表明，土壤含水率（SM1和SM2）的模擬結果達到穩定所需的時間最短，熱通量（LE、H和土壤熱通量G）相較需要更長的時間，其中顯熱通量的預熱結果不穩定，預熱多輪后其相對誤差成周期性循環，取進入循環前一年作為最終所需預熱輪數；NEE 模擬穩定所需時間需要100 a 以上，與其值較?。?0-4量級），達到數值穩定更為困難不無關系。此結果與密西西比流域的研究結果存在差異，CAI 等[5]認為土壤水預熱需要8 a，而顯熱與潛熱通量預熱的時間更短，大約需要4 a。

圖2 為不同抽樣次數的輸出變量敏感參數分析結果。結果表明：對于不同抽樣次數而言，隨著抽樣次數的減少，大部分參數的全局敏感度大小發生變化，尤其是植被相關參數（如QE25、SLA、LTOVRC 等）變化較為顯著，全局敏感度較高的參數變化較明顯。盡管如此，無論抽樣次數如何變化，對輸出變量不敏感的參數依舊只有非常小的全局敏感度。因此，在不考慮敏感性參數精確排序的條件下，若只為確定哪些參數更重要、更為敏感，拉丁超立方抽樣法（LHS）的抽樣次數M 可以取參數數量N 的5 倍左右，這與HUO 等[7]的研究結論具有一致性。下文中M 取為250 以節約計算成本。

圖2 不同抽樣次數的輸出變量參數敏感參數分析Fig.2 Parameter sensitivity analysis of output variables for different sampling number

2.2 站點參數敏感性結果分析

為了對參數的敏感性進行區分，簡單地將全局敏感度高于0.10 的參數劃分為高敏感性參數，將全局敏感度高于0.05且不高于0.10 的參數劃分為中敏感性參數，禹城站與欒城站的敏感性參數劃分結果見表4。

表4 禹城站與欒城站敏感性參數匯總（btr = 2 : CLM）Tab.4 Summary of sensitivity parameters for two stations（btr = 2 : CLM）

綜合2 個站點的結果，由于LE 與H、SM1 與SM2 涉及的物理過程一致，因此其敏感性參數也大致相同。對LE 與H 來說，模擬過程不僅包含水量平衡還涉及到冠層能量平衡過程，因此它們對植被動態生長過程中的植物器官周轉率LTOVRC及25 ℃時葉綠素吸收的光量子效率QE25也很敏感。LE、H與土壤有關的敏感參數主要為BEXP、SMCMAX、SMCREF、SMCWLT，這些參數都與植被作用下的土壤水狀態密切相關。故而它們也是SM1 的敏感參數。與SM1 不同的是，由于SM2所處位置更深，土壤水與植被的相互作用更小，因此對凋萎系數SMCWLT 及田間持水量SMCREF 這類與植物相關的土壤參數敏感性較弱，而對土壤孔徑分布指數BEXP 及土壤孔隙率SMCMAX這類土壤物理性質參數更敏感。

NEE 的敏感性參數為微生物呼吸參數MRP 和葉組織轉化速率LTOVRC，但NEE 對MRP 的敏感性顯著高于其他敏感參數，這一方面可以證實農田生態系統的固碳潛力主要集中在農田土壤[17]，另一方面與模型對農田生態系統的作物碳同化過程模擬能力較弱有關，Noah-MP 模型的模擬結果低估了冬小麥與夏玉米生長旺季（4 月和8 月）時農田生態系統的碳匯作用（見圖3）。同時，LTOVRC 是一個與葉片物質量更新周轉進入土壤有機碳有關的參數，這同樣能說明在模型模擬農田生態系統NEE 過程中，農田土壤碳庫的變化過程占主導，而與光合作用碳同化過程有關的作物碳庫變化過程較弱。與光合作用相關的VCMX25、QE25 及SLA，雖然也是重要參數，但NEE 對它們敏感性較低。

圖3 禹城站與欒城站NEE模擬值與真實值的比較Fig.3 Comparison of simulated and observed NEE for two stations

整體上，2個站點的敏感參數基本相同（見圖4）。結合所有輸出變量來看，2 個站點的區別體現在土壤孔隙率SMCMAX、田間持水量SMCREF、凋萎系數SMCWLT 和土壤孔徑分布指數BEXP 這4 個參數上。相比禹城站，欒城站模擬的LE、H 對SMCMAX 的變化更不敏感，而SM1、SM2 受SMCREF 和SMCWLT 的影響更大。SMCMAX 用于確定土壤水達到靜水平衡時的地下水位，其值為劃分土壤水與地下水的閾值，體現了土壤水與地下水間的轉換關系。禹城站SMCMAX 更為敏感，說明禹城站的地下水位更淺，地下水對土壤水分狀態的影響較大。站點資料顯示禹城站地下水位埋深一般為1.5～4.0 m，欒城站地下水位埋深大部分在15 m以下。SMCWLT 用來刻畫水分脅迫對植物生長、器官凋亡等生理過程的影響，SMCREF則決定了降水與灌溉水的下滲與土壤的疏干過程，BEXP 為土壤水分特征方程的重要參數，與土壤含水狀態密切相關。欒城站的淺層土壤水（SM1、SM2） 對SMCWLT、SMCREF、BEXP 更為敏感，說明欒城站的氣候更為干旱（禹城站該輪作年干旱指數為1.929，欒城站為2.540），土壤含水率受降水與灌溉的影響大，作物生長更可能受到水分脅迫，由氣象數據根據FAO 推薦的Penman-Monteith 方法計算的日參考作物需水量與降水見圖5。

圖4 禹城站與欒城站各輸出變量參數敏感參數分析（btr = 2 : CLM）Fig.4 Parameter sensitivity analysis of each output variable for two stations

圖5 禹城站與欒城站的日參考作物需水量與降水Fig.5 Daily reference crop evapotranspiration and precipitation for two stations

2.3 氣孔阻抗的土壤水分限制因子參數化方案的比較

在植被覆蓋情況下，氣孔阻抗不僅決定了植物蒸騰過程中水汽進入大氣的能力，而且也決定了植物的光合作用中二氧化碳進入葉肉組織的能力。植物的氣孔開閉受到土壤水分狀態的調控，在土壤含水率較低時，植物能感應根區土壤水

分狀態并對氣孔進行調節，防止自身過分失水而死亡。在Noah-MP 模型中，土壤水分對氣孔開閉影響的描述存在于氣孔阻抗的土壤水分限制因子參數化方案之中。2.2 節中選用的是基于土壤水勢的CLM 方案，本節將結合2.2 節的計算結果，比較CLM 方案與基于土壤含水率的Noah 方案在參數敏感性上的異同。與2.2節的劃分原則一致，禹城站與欒城站的中/高敏感參數見表5。考慮參數抽樣引起的變異性，在比較不同氣孔阻力參數化方案時，與比較站點間差異一樣，重點關注全局敏感度的相對大小變化，分析那些相對變化明顯的參數。

表5 禹城站與欒城站敏感性參數匯總（btr = 1 : Noah）Tab.5 Summary of sensitivity parameters for two stations（btr = 1 : Noah）

結合2 個站點來看（圖6、圖7），使用不同氣孔阻力參數化方案對SM1、SM2、NPP 的參數敏感性幾乎無影響。而將參數化方案從CLM 更換為Noah 2 個站點都出現了以下現象：①對LE 來 說LTOVRC、QE25的敏感性降低，SMCREF、SMCWLT 的敏感度升高；②對H 來說QE25 的敏感性降低。由于控制氣孔阻力的土壤水分因子參數化方案的本質區別在于Noah 方案基于土壤含水率而CLM 方案基于水勢，在Noah 方案中氣孔阻抗由SMCREF 和SMCWLT 直接決定，因此Noah 方案中LE 對這2 個參數更敏感。模型中的水碳耦合的關鍵樞紐在于葉片的氣孔，氣孔開放程度不僅影響植物蒸騰還影響光合過程，進而間接影響作物的生長。LTOVRC、QE25 的敏感性降低能在一定程度上反映出LE 對土壤水勢的響應靈敏度高于土壤含水率。H 與水通量的聯系沒有LE 那么緊密，因此與土壤水直接相關的SMCREF、SMCWLT的敏感度變化不明顯。而H 對QE25 的敏感性降低但LTOVRC 的敏感度變化不明顯，這反映出氣孔開閉對H 的影響主要是通過影響植物的光合過程實現。對比圖6 與圖7 可以看出，在氣候更為干旱的欒城站，不同參數化方案下BEXP 的敏感度發生了非常明顯的變化，原因在于BEXP直接參與CLM方案的水勢計算過程，因此在氣候干旱地區使用該方案時需要重視對BEXP的校準。

圖6 禹城站不同參數化方案的各輸出變量參數敏感參數分析Fig.6 Parameter sensitivity analysis of each output variable for different parameterization schemes at Yucheng station

圖7 欒城站不同參數化方案的各輸出變量參數敏感參數分析Fig.7 Parameter sensitivity analysis of each output variable for different parameterization schemes at Luancheng station

3 討 論

華北平原主要的土地利用類型為耕地，是一個大型的農田生態系統，水熱碳通量的變化受到種植、管理等人類活動的影響。與其他植被覆蓋類型相比，陸面模式在模擬農田時存在著巨大的挑戰。為了準確刻畫華北平原農田生態系統中的水、能量及CO2通量變化過程，降低參數不確定性對陸面模式的影響，需要開展模型參數敏感性分析。

本研究在ARSENAULT 等[6]的基礎上對用于敏感性分析的單站點模型進行修改，完善預熱前后淺層土壤中的速效碳庫FASTCP、深層土壤中的遲效碳庫STBLCP 等模型狀態變量的繼承關系，代碼修改前后欒城站的敏感性分析結果見圖8。結果顯示，模型代碼修改前后，輸出變量LE、H、SM1 與SM2對不同參數的敏感程度基本一致，NEE 的敏感性參數發生巨大變化。模型代碼修改后，LE、H 對土壤參數的敏感性得以凸顯，揭示了熱通量與土壤水的潛在聯系，參數BEXP 成為最敏感參數，與HUO 等[7]的研究結果一致。BEXP 是土壤的形狀參數，與土壤的孔隙結構有關，同時其作為土壤水分特征方程的重要參數，能夠反映土壤的持水性，與陸氣間水分交換息息相關。對于NEE 來說，模型代碼修改后，敏感性參數不再是植物生理相關參數，微生物呼吸參數MRP 的全局敏感性指標顯著高于其他參數，這恰能體現農田生態系統碳匯功能主要依賴土壤固碳而不是農作物生物量固碳。同時，本研究發現為使凈生態系統碳交換量（NEE）達到穩定值需要經歷100 a 以上的預熱過程。修改代碼并對模型進行充分預熱后，為土壤微生物呼吸提供能量的淺層土壤中的速效碳庫FASTCP不斷積累，這也導致MRP 的敏感性增強。值得注意的是，NEE 預熱穩定后FASTCP 的值為102量(t/hm2)，遠高于第2 次全國土壤普查數據中農田表層土壤有機碳密度（26.6～32.5 t/hm2）[17]，這是Noah-MP 陸面模式對NEE 等碳循環物理量模擬能力欠缺導致的，開發者應重視對碳循環模塊的優化。在對比2種氣孔阻抗的土壤水分限制因子參數化方案對模型敏感性的影響時發現，基于土壤水勢的CLM 方案能更好地將干旱地區植物根區的缺水狀態反映到LE 的變化上；而YANG 等[18]也得出了類似的結論，認為基于土壤含水率的Noah 方案在降水充沛的西雙版納站能獲得最接近觀測值的模擬結果。ARSENAULT[6]等同樣發現在半干旱的草原與農田生態站CLM方案對BEXP 的敏感度高于Noah 方案，并認為PSISAT 在不同方案間的敏感度差異也源自同樣的機理。本研究并未發現PSISAT 在不同方案間的敏感度差異，這可能是因為選擇的下墊面單一，只聚焦于華北平原生態系統。

圖8 代碼修改前后欒城站參數敏感性分析結果Fig.8 Parameter sensitivity analysis results at Luancheng station before and after code modification

本研究揭示了Noah-MP 模型模擬華北平原農田生態系統水熱碳通量及土壤含水率時的敏感參數，并對比了不同氣孔阻抗的土壤水分限制因子參數化方案下參數敏感度，對參數優化與參數化方案選擇有著指導意義。但本研究仍存在一些不足，例如，只考慮用戶自定義參數的敏感性，而未考慮模型嵌入參數，CUNTZ 等[19]發現對潛熱通量、總徑流及各徑流分量而言，許多模型嵌入參數甚至比用戶自定義參數更敏感，因此未來可以將參數敏感性分析的范圍拓展至模型嵌入參數。

4 結 論

（1）若不需要對參數按敏感性強弱進行嚴格排序，在使用拉丁超立方抽樣法（LHS）篩選模型的敏感性參數時，為節約計算成本，抽樣次數可以取為參數的5倍左右。

（2）用Noah-MP 模型模擬華北平原農田生態系統典型站點時，LE、H 的敏感參數有：LTOVRC、QE25、BEXP、SMCMAX、SMCREF、SMCWLT。其中，BEXP、SMCMAX、SMCREF、SMCWLT 也是SM1 和SM2 的敏感參數。NEE 對MRP、LTOVRC 的敏感性顯著高于VCMX25、QE25 及SLA 等與作物光合作用相關的參數，體現農田生態系統碳匯功能主要依賴農田土壤固碳。

（3）在氣候干旱的地區需要重視對土壤參數SMCMAX、SMCREF、SMCWLT、BEXP 的校準，尤其在使用CLM 氣孔阻抗的土壤水分限制因子參數化方案時，BEXP 的校準需要特別關注。

致謝：感謝禹城農業綜合試驗站與欒城農業生態系統試驗站為本研究提供的數據支持，感謝Arsenault Kristi R.在GitHub 平臺貢獻的開源代碼，使本研究能在此基礎上提出改進。