基于叶面积指数改进的直角双曲线模型在玉米农田生态系统中的应用

孙敬松，周广胜

(1.中国科学院植物研究所植被与环境变化国家重点实验室，北京 100093;2.中国气象科学研究院，北京100081;3.中国科学院研究生院，北京 100049)

陆地植被的总初级生产力(GPP)(光合CO2的同化与光呼吸平衡)导致每年90—130 Tg的碳被同化，每年15%的大气碳被陆地植被光合作用固定［1］。太阳辐射是植物光合的主要驱动因子，是影响陆地生态系统生产力与碳收支的重要环境因子。通常，能被植物光合作用直接吸收利用的那部分太阳辐射，波长在0.4—0.7μm范围内，称为光合有效辐射(PAR)［2-3］。PAR常被用作植被净初级生产量模型、气候模型、水文模型、生物地球化学模型等的关键变量［4］，弄清GPP与PAR之间的关系是准确评估陆地生态系统的生产力及其碳收支的关键。研究表明，植被的光合作用与光合有效辐射之间的关系可以用直角双曲线Michaelis-Menten模型很好地描述(式 3)［5-6］。

以往研究中，利用直角双曲线模型对生态系统GPP进行模拟或插补时，模型的两个参数Amax和α值是不变的［5-6］。然而，有研究表明在不同的生长季进行模拟，这两个参数值是不同的［7-8］。为了消除模型参数季节性波动对模型模拟结果的影响，许多研究通常是在较短的时间尺度上利用直角双曲线模型进行模拟(如半个月，10d或5d)［7-9］。尽管如此，这种方法过于繁琐而且缺乏机理性，尤其是在短的时间尺度上进行模拟，当涡相关数据量缺失过多时，少量的数据不足以准确拟合直角双曲线模型。

农田生态系统碳吸收大于或等于自然生态系统［6，10］，对全球碳循环有着极为重要的作用。因此，以玉米(Zea mays L.)农田生态系统为例，探讨态系统的环境因子和生物因子在不同生长季对直角双曲线模型中两个参数Amax和α值的影响。

本研究基于中国气象局沈阳大气环境研究所锦州农田生态系统野外观测站连续5a的观测资料(2004—2008年)，分析了反映GPP与PAR关系的Michaelis-Menten直角双曲线模型参数Amax和α值的季节变化，进而阐明环境因子和生物因子(温度、水分、VPD、辐射、叶面积指数(LAI))对模型参数季节变动的影响。在此基础上，本研究找出影响模型参数季节变化的主导因子及其过程，并将影响过程整合到直角双曲线模型中改进模型，并对玉米农田生态系统GPP进行模拟和比较。

1 材料和方法

1.1 研究地点

研究地点位于辽宁省锦州市中国气象局锦州农田生态系统野外观测站(41°49'N，121°12'E，海拔17.0 m)，属温带大陆性季风气候，年极端最高气温35℃，极端最低气温－29.9℃，全年平均气温为7.8—11.0℃，年平均降雨量 550—630 mm［11］。

玉米农田观测场内安装有涡度相关系统与微气象梯度仪，可以对水热、碳通量及其太阳辐射、温度、湿度等环境因子进行连续观测，涡相关和微气象梯度仪安装于一个25 m×25 m的围栏内。围栏内以及周围地势平坦开阔，无任何高大建筑物遮蔽，试验田的土壤为典型棕壤，生长季种植作物均为春玉米(中晚熟、稀植品种，华单4号和新科10号)，植株平均密度为4.3株/m2，经实地观测，植株高度最高可达3 m左右。大范围均一的下垫面满足试验要求盛行风向的风浪区长度。在非生长季农田下垫面无作物，地表为裸土覆盖;涡度相关仪安装于距地3.5 m处，不再随作物生长调节安装高度。自2004年6月开始自动观测至今。

1.2 观测系统

1.2.1 涡度相关通量观测与数据处理

通量观测采用开路涡度相关通量观测系统。该系统由CR5000数据采集器、CSAT3超声风速仪(Campbell公司)、Li-7500 CO2/H2O分析仪(Li-Cor公司)和1个1 G的PC卡组成。CR5000数据采集器主要是控制测量、运算及数据存储;CSAT3超声风速仪用以测量三维风速和超声虚温;Li-7500测量CO2绝对密度、H2O绝对密度和大气压力。整个系统安装于距地面3.5 m高度处。仪器采样频率为10 Hz，由数据采集器自动存储10 Hz原始数据。

CR5000数据采集器有2 M内存，可以保存一段时间的观测数据。系统工作时由涡度相关原理计算在线通量，并存储30 min的CO2通量(Fc)、潜热通量(LE)与显热通量(Hs)和时间系列数据，系统同时对在线通量自动进行虚温订正与空气密度变化订正。观测数据由于受到天气因素、地形条件以及仪器本身故障等影响，需要先进行数据质量分析与控制，才能用于研究分析。数据质量控制以涡度相关法通量观测原理以及涡度相关仪器本身局限性为依据，对30 min数据进行剔除［12］。

在通量观测研究中，通常将CO2进入生态系统标记为负值，反之为正值，因此当NEE为负值时，表示生态系统吸收CO2，反之为释放CO2。

1.2.2 太阳辐射与土壤温度测定

观测样地安装有微气象梯度观测塔，太阳辐射来自气象梯度观测塔3.5 m高度的CNR1净辐射仪和LI190SB光量子传感器，土壤温度来自观测塔107热敏电阻温度探头观测的5 cm土壤温度。根据太阳辐射量区分白天和夜间(太阳总辐射＜20 wm－2)［12］。

1.2.3 叶面积指数测定

每隔20 d左右一次性采集完整的5株标准株玉米鲜样，量取每株叶片的叶面积，叶面积指数(LAI)为单位面积上所有叶片表面积的总和与相应土地面积之比。

1.3 涡相关通量数据的拆分和总初级生产力(GPP)的估算

生态系统碳交换(NEE)受光合与呼吸作用调控，由于涡度相关系统直接测量的NEE需要利用公式(1)进行组分拆分以估算GPP［13-14］:

本研究利用夜间数据通过Lloyd-Taylor模型建立R e与土壤5 cm温度的关系(式2)［15］估算白天生态系统呼吸(自养呼吸和异养呼吸)Reco:

式中，T0为常数，一般设为 －46.02℃，E0为活化能参数 决定了呼吸的温度敏感性，通过模型的拟合得出，Reco，ref为参考温度下的呼吸值，通过模型拟合得出，Tref为参考温度，一般设为10℃，T为5 cm土壤温度。

1.4 数据分析

玉米一般在4月下旬至5月上旬播种，5月下旬至6月初进入三叶期，9月中旬成熟。为方便分析直角双曲线模型在不同季节的模拟情况，选取5年(2004—2008)同一时段的(6月1日—9月15日)的观测资料和直角双曲线模型(式3)进行分析:

式中，GPP(mg CO2m－2s－1)为生态系统总初级生产力，α是当PAR接近于0时的生态系统光能利用率［16］，称为生态系统的表观初始光能利用率，PAR为入射到冠层上的光合有效辐射(μmol photon m－2s－1)，Amax为生态系统光合作用最大时的碳交换(mg CO2m－2s－1)。统计分析软件和函数绘图分析软件则主要利用SPSS 15.0(SPSS，Chicago，IL，USA)和 Origin8.0(Microca，OriginLab，USA)。

2 结果分析

2.1 数据分析

2.1.1 模型参数的季节动态

传统的光合作用模拟一般是利用直角双曲线模型(式3)在整个生长季或15d的时间尺度上进行，在模拟的时间窗口中，直角双曲线的参数通常设为常数［5，13-14，17］。为此，首先对玉米农田生态系统的GPP进行整个生长季和每隔半月的直角双曲线模拟(表1)。

表1 直角双曲线模型的参数表观初始光能利用率(α)和生态系统光合作用最大时的碳交换(A max)、相关系数(R2)、样本数(n)和均方根误差分析(RMSE)Table1 Simulated results of the rectangular hyperbola model，includes parameters of the apparent quantum yield(α)and the value of NEE at light saturation(A max)，coefficient of determination(R2)and root-mean-square error(RMSE)

不同生长阶段的玉米农田生态系统GPP模拟中，直角双曲线模型参数α和Amax的值发生了波动，其中α值的波动较小(变异系数=28.85%)，说明Amax除了受PAR的影响以外还可能受生物量、温度、水分等因子季节变化的影响，而α值较为稳定。这种影响也使得直角双曲线模型在整个生长季对玉米农田生态系统GPP的模拟结果较差(表1)。

2.1.2 模型参数α和Amax季节动态的影响因子分析

为了弄清影响玉米农田生态系统GPP-PAR之间直角双曲线关系的主导因子，假定直角双曲线模型中的参数α和Amax分别受不同环境或生物因子的影响而产生波动。为了验证这个假设，选取有叶面积指数观测的5d之内的数据(有叶面积观测的和前后2d共5d的数据)，首先对双曲线模型中的参数α、Amax与气温、土壤温度、空气湿度、土壤水分、太阳辐射、饱和水汽压差(VPD)、LAI等因子的均值进行pearson相关分析发现，Amax和LAI具有显著的相关关系(P=0.002)(图1)，而α与其它因子的相关关系均不显著，这进一步说明α值的保守性。

图1 A max与LAI线性回归分析Fig.1 The linear regression between A max and LAI

2.1.3 考虑叶面积指数的直角双曲线模型

Amax与LAI的线性关系表明，直角双曲线模型的参数Amax季节动态可能由LAI引起。因此，可以利用它们的线性关系来预测模型的参数变化，以准确地估算生态系统的GPP。基于这一假设，利用Amax与LAI的线性关系改进直角双曲线模型(式4):

式中，GPP(mgCO2m－2s－1)为生态系统的总初级生产力，PAR为入射到冠层上的光合有效辐射(μmolphoton m－2s－1)，LAI为叶面积指数，其中，生态系统的表观初始光能利用率α，Amax=a LAI+b为生态系统光合作用最大时的碳交换，a、b为方程的参数。

模型的变量PAR可以通过常规气象观测得到半小时、小时、天尺度的值，而LAI在每个生长季的观测次数有限，会在一定程度上限制了模型的准确度。所以，引入一个Logistic模型(式5)［18］，将整个生长季LAI的观测频度拓展到天的时间尺度。

式中，LAI为每天的叶面积指数，t为玉米播种后天数

利用玉米农田5a的LAI观测资料对方程进行拟合。结果表明，Logistic模型对LAI的模拟具有较高的准确性(图2)。

利用PAR(半小时尺度的数据)和由Logistic模型获得的LAI(时间尺度为天)等数据，对改进的直角双曲线模型(式5)进行拟合分析，详见表2。

表2 改进直角双曲线模型的参数a、b、α、相关系数(R2))和均方根误差分析(RMSE)Table2 Simulated results of the improved rectangular hyperbola model，includes parameters of model(α，a，b)，coefficient of determination(R2)and root-mean-square error(RMSE)

在整个生长季上进行模拟，未改进的直角双曲线模型仅能解释玉米农田生态系统27%的GPP波动(表1)，而改进的直角双曲线模型对玉米农田生态系统GPP的模拟结果可以解释玉米农田生态系统55%的GPP波动(表2)，并能反映模型参数随LAI变化的过程与机理。

图2 玉米农田生态系统LAI的生长曲线模拟Fig.2 The simulation of daily LAI

3 讨论

利用模型模拟以外的数据，也就是没有相应的叶面积指数观测的那一时段的碳通量数据，对改进的直角双曲线模型和未改进的直角双曲线模型模拟结果进行了比较。生长季的叶面积指数变化引起了玉米农田生态系统冠层光合作用对PAR响应的变化，从而影响直角双曲线模型的参数Amax值。因此，将直角双曲线模型的Amax参数作为常数模拟玉米农田生态系统总初级生产力将产生较大误差(表1)。虽然通过减小模型模拟的时间尺度(例如半个月)，可以提高直角双曲线模型模拟的准确度(图3C)，但这种模拟方法不仅比较繁琐，而且缺乏机理性，难以说明相邻的两个时间段内模型参数的突变。改进的直角双曲线模型耦合了动态的LAI对玉米农田生态系统光合作用的影响过程，阐明了直角双曲线模型的重要参数Amax对生物因子的响应，使得玉米农田生态系统光合速率对PAR的响应曲线在日时间尺度产生变异。通过线性回归模型比较发现，改进的直角双曲线模型(图3C)和半月时间尺度的模拟结果(图3B)较接近，均好于直角双曲线模型在整个生长季尺度的模拟(图3A)。

图3 观测值和模型模拟结果的比较Fig.3 Relationship between measured and simulated GPP

尽管本研究中，改进的直角双曲线模型与半个月时间尺度上未改进的直角双曲线模型模拟的准确性相近，但是改进的直角双曲线模型不仅简捷，而且可解释参数Amax的连续变化，尤其是当缺失的通量数据过多时，直角双曲线模型不能在半个月或更短的时间尺度上进行拟合时，改进的直角双曲线模型可以较为准确地对缺失数据进行插补。

4 结论

为了消除直角双曲线模型参数季节波动对模型模拟结果的影响，通常在较短的时间尺度上对GPP进行模拟和插补［7-9］。然而，这种方法过于繁琐缺乏机理性，尤其是在短的时间尺度上，利用少量的碳通量数据难以准确拟合直角双曲线模型。

LAI是控制GPP-PAR关系的关键因子，它的变化可以解释Amax值的变化。改进的模型引入了LAI和参数Amax的关系，提高了模型在整个生长季尺度上对GPP模拟的准确度，从机理上阐明了直角双曲线模型在不同的时间尺度上模拟的准确度问题，为准确评估陆地生态系统碳循环提供了借鉴。

生态系统碳收支动态评估是全球变化研究的关键之一［18］，改进的直角双曲线模型可以方便、准确地评估生态系统的GPP。同时，结合卫星遥感获得的植被指数数据，通过与叶面积指数的关系，可以较方便地进行区域或全球尺度上陆地生态系统的碳收支评估。本研究中，环境因子如土壤、湿度和VPD等没有影响模型参数的变化，这些因子可能是通过其它机理影响生态系统的碳收支，因此，将来的研究将着重考虑这些因子对生态系统碳收支的影响过程和机理。

［1］ Williams M，Rastetter E B，Fernandes D N，Goulden M L，Shaver G R，Johnson L C.Predicting gross primary productivity in terrestrial ecosystems.Ecological Applications，1997，7(3):882-894.

［2］ Udo SO，Aro T O.Global PAR related to global solar radiation for central Nigeria.Agricultural and Forest Meteorology，1999，97(1):21-31.

［3］ Jacovides CP，Timvios F S，Papaioannou G，Asimakopoulos D N，Theofilou CM.Ratio of PAR to broadband solar radiation measured in Cyprus.Agricultural and Forest Meteorology，2004，121(3/4):135-140.

［4］ Thomas V，Finch D A，McCaughey J H，Noland T，Rich L，Treitz P.Spatial modelling of the fraction of photosynthetically active radiation absorbed by a boreal mixedwood forest using a lidar-hyperspectral approach.Agricultural and Forest Meteorology，2006，140(1/4):287-307.

［5］ Ruimy A，Jarvis P G，Baldocchi D D，Saugier B.CO2fluxes over plant canopies and solar radiation:a review.Advances in Ecological Research，1995，26:1-68.

［6］ Goulden M L，Daube B C，Fan SM，Sutton D J，Bazzaz A，Munger JW，Wofsy SC.Physiological responses of a black spruce forest to weather.Journal of Geophysical Research，1997，102(D24):28987-28996.

［7］ Gu L，Baldocchi D，Verma SB，Black T A，Vesala T，Falge E，Dowty P.Advantages of diffuse radiation for terrestrial ecosystem productivity.Journal of Geophysical Research，2002，107(4050):1155-1168.

［8］ Polley W H，Emmerich W，Bradford JA，Sims P L，Johnson D A，Saliendra N Z，Svejcar T，Angell R，Frank A B，Phillips R L，Snyder K A，Morgan J A.Physiological and environmental regulation of interannual variability in CO2exchange on rangelands in the western United States.Global Change Biology，2010，16(3):990-1002.

［9］ Zhang L M，Yu G R，Sun X M，Wen XF，Ren C Y，Fu Y L，Li Q K，Li Z Q，Liu Y F，Guan D X，Yan JH.Seasonal variations of ecosystem apparent quantum yield(α)and maximum photosynthesis rate(P-max)of different forest ecosystems in China.Agricultural and Forest Meteorology，2006，137(3):176-187.

［10］ Law B E，Falge E，Gu L，Baldocchi D D，Bakwin P，Berbigier P，Davis K，Dolman A J，Falk M，Fuentes JD，Goldstein A，Granier A，Grelle A，Hollinger D，Janssens I A，Jarvis P，Jensen NO，Katul G，Mahli Y，Matteucci G，Meyers T，Monson R，Munger W，Oechel W，Olson R，Pilegaard K，Paw K T，Thorgeirsson U H，Valentini R，Verma S，Vesala T，Wilson K，Wofsy S.Environmental controls over carbon dioxide and water vapor exchange of terrestrial vegetation.Agricultural and Forest Meteorology，2002，113(1/4):97-120.

［11］ Han G X，Zhou GS，Xu Z Z，Yang Y，Liu JL，Shi K Q.Biotic and abiotic factors controlling the spatial and temporal variation of soil respiration in an agricultural ecosystem.Soil Biology and Biochemistry，2007，39(2):418-425.

［12］ Papale D，Reichstein M，Aubinet M，Canfora1 E，Bernhofer C，Kutsch W，Longdoz B，Rambal S，Valentini R，Vesala T，Yakir D.Towards a standardized processing of Net Ecosystem Exchange measured with eddy covariance technique:algorithms and uncertainty estimation.Biogeosciences，2006，3(4):571-583.

［13］ Stoy P C，Katul GG，Siqueira M B S，Juang JY，Novick K A，Uebelherr JM，Oren R.An evaluation of models for partitioning eddy covariancemeasured net ecosystem exchange into photosynthesis and respiration.Agricultural and Forest Meteorology，2006，141(1):2-18.

［14］ Reichstein M，Falge E，Baldocchi D，Papale D，Aubinet M，Berbigier P，Bernhofer C，Buchmann N，Gilmanov T，Granier A，Grünwald T，Havrankova K，Ilvesniemi H，Janous D，Knohl A，Laurila T，Lohila A，Loustau D，Matteucci G，Meyers T，Miglietta F，Ourcival JM，Pumpanen J，Rambal S，Rotenberg E，Sanz M，Tenhunen J，Seufert G，Vaccari F，Vesala T，Yakir D，Valentini R.On the separation of net ecosystem exchange into assimilation and ecosystem respiration:review and improved algorithm.Global Change Biology，2005，11(9):1424-1439.

［15］ Lloyd J，Taylor J A.On the temperature-dependence of soil respiration.Functional Ecology，1994，8(3):315-323.

［16］ Wofsy SC，Goulden M L，Munger J W，Fan S M，Bakwin P S，Daube B C，Bassow S L，Bazzaz F A.Net exchange of CO2in a midlatitude forest.Science，1993，260(5112):1314-1317.

［17］ Goulden M L，Daube B C，Fan SM，Sutton D J，Bazzaz A，Munger JW，Wofsy SC.Physiological responses of a black spruce forest to weather.Journal of Geophysical Research，1997，102(D24):28987-28996.

［18］ Nightingale JM，Coops N C，Waring R H，Hargrove WW.Comparison of MODISgross primary production estimates for forests across the U.S.A.with those generated by a simple process model，3-PGS.Remote Sensing of Environment，2007，109(4):500-509.

参考文献:

［18］ 王信理.在作物干物质积累的动态模拟中如何合理运用Logistic方程.中国农业气象，1986，7(1):14-19.