基于SiB2模型的低丘红壤区花生根域土壤水分动态模拟与数据同化

张　乾，景元书

(1.气象灾害预报预警与评估协同创新中心/江苏省农业气象重点实验室，南京 210044；2.南京信息工程大学，南京 210044)

0　引　言

土壤水分是研究作物生长与地气相互作用之间联系的重要指标，由于其在空间和时间上变化非常大，因此，建立准确的陆面过程模型评估陆面条件尤其是土壤水分变化在陆地大气相互作用中的表现，研究大尺度和流域尺度的水分和能量循环，获取完备的地球表层系统的时空信息[1]，对提升区域尺度天气预报精确性有着积极的意义。根据基本气象资料对若干段时间内作物根域土壤水分的变化进行准确的预估，有利于更深层次地研究和理解地球表层生物物理过程。同时，陆面模型在进行模拟时，模拟误差会随着时间不断积累，因此可以建立数据同化方案，在观测时间段内对模型的模拟进行校正以提高模拟精度。

Sellers 等人于1996年建立的SiB2(Simple Biosphere Model version 2，简单生物圈模型)，是目前研究单层植被中水、能量以及CO2通量变化较成熟的陆面模型之一。SiB2在全球尺度下使用卫星图像，获得了多种类型植被的物候学资料作为模型运行的重要参数，应用能量、动量和质量守恒定律, 能够方便有效地模拟出土壤、大气和生物圈诸多参数以及陆面(植被和裸地)和大气之间的相互作用，该模式已被嵌入其他区域和全球模式中[2]。

南方低丘红壤区是我国水土流失较为严重的地区之一，由于红壤性质上的酸、瘦、黏等弱点，红壤分布区域降水时空分布的不均匀，高温伏旱频繁，以及不合理开发利用造成的水土流失、土壤退化、土壤污染等，导致红壤地区的生态环境恶化，红壤资源潜在的生产能力得不到应有的发挥，使整个地区农业及经济持续发展受到严重影响[1]。花生是我国南方红壤地区重要的经济作物之一，利用模型以及相应的数据同化方案准确地模拟花生根域土壤水分变化，有利于指导当地农业生产，趋利避害。

1　研究地区与研究方法

1.1　研究区概况

研究区位于江西省鹰潭市余江县下辖刘垦场三分场的红壤生态试验站( 116°55′E，28°15′N) 。该地区地形以低丘岗地为主，南北有少量丘陵，平均海拔37.56 m，坡度3°～8°，土壤多为红砂岩或红黏土发育的普通红壤[4]，是典型的我国南方低丘红壤区之一。属于亚热带湿润季风气候，四季分明、气温偏高。年平均气温17.6 ℃，年平均降水量1 788.8 mm，主要集中在春末夏初时期，年平均日照时数1 739.4 h。研究区域内，经济作物以水稻、花生、柑橘、板栗为主，植被分布主要以樟科、茶科等常绿阔叶林以及部分竹类、藤本、蕨类等构成。

1.2　数据来源及处理

土壤湿度和地表温度数据来源于鹰潭市余江县刘垦场三分场的FDR水分传感器观测的2015年花生地的水分数据。进行统计分析前，对土壤水分数据进行单样本k-s检测，其中6-7月的数据采用平方根反正切函数进行转化，最终使所有数据均符合正态分布，可以进行统计分析[5]。根据对花生地实地测量，得到适用于试验区土壤特性和植被条件的 SiB2所需要的全部参数。

1.3　试验设计

在试验区代表性地段，选择上坡、下坡建造的20 m×5 m(100 m2)的试验小区，种植的花生品种为粤油551[6]，小区由水泥板围成，水泥板高出地面30 cm，以阻挡小区外径流流入[4]。FDR水分传感器每30 min测定一次10～100 cm深度的土壤温度以及水分含量。分别选取5月10-19日、6月8-17日、8月17-26日3个均为期10 d的时间段作为花生幼苗期、开花期、饱果成熟期生长阶段的代表时间段，将3个时间段的逐小时深度分别为10、15、20 cm的土壤湿度作为观测值，分别在SiB2中代入测定好的3个时间段的花生形态学参数、模型主程序运行参数等，选择5月10日、6月8日、8月17日3个日期0时刻仪器测量得到的植被温度、地表温度、深土温度、冠层空气温度、参照高度温度、土壤表层水分含量、植被根域水分含量、深层土壤水分含量数据分别作为花生幼苗期、开花期、饱果成熟期SiB2运行的初始值，最后代入相应的自动气象站气象资料。其中，修正SiB2长波辐射的计算，得到土壤水分模拟值后，引入集合卡尔曼滤波进行数据同化以调高模拟精度。

1.4　误差统计

对于模型模拟和同化结果，采用均方根误差(RMSE, Root Mean Square Error)、平均误差(AVE, Average Error)和相对误差(RE, Relative Error)三种误差分析方法[7]。

(1)

(2)

(3)

式中：N为整个观测期的时长；Obst为t时刻的观测值；Xt为t时刻的模型模拟和同化结果[7]。

2　SiB2和数据同化方案

2.1　SiB2简介

陆面过程作为气候-陆地系统间重要的反馈机制之一，是天气、气候、水文和环境预报的关键因素。针对不同的陆地下垫面，具有的地表反照率、比辐射率、地面粗糙度等理化性质以及相应的植被类型各不相同，局地气候受其影响的敏感性也不同的陆面过程难点，1996年，Seller等人在SiB1的基础上建立了SiB2，融入了更科学的植被冠层光合作用传导模型，分别描述CO2和水汽在叶片内部和外部的传输过程，克服了SiB1单纯使用经验公式的缺陷。Sellers还对水文子模型进行了改进，以确切地描述基本径流，并对土壤剖面内的层间交换进行了更可靠的计算，模式定义了9种植被类型和7种土壤类型，在模拟地表能量收支、径流传输等取得了好的效果。

SiB2 利用气象驱动数据，可以模拟一定时间一定区域内的净辐射、感热、潜热、地表热通量、地表有效辐射温度、土壤湿度变化。在模拟地面以下土壤的热传导和水热输送过程中，SiB2定义了表层、根域和深层3 层土壤。公式(4)～(6)为每层土壤的水分控制方程[7]。

(4)

(5)

(6)

式中：θi(i=1,2,3)分别为表层、根区、深层土壤中液态水的体积含水量，m3/m3；ρw为液态水的密度，kg/m3；Di(i=1, 2, 3)为每层土壤厚度，m；Qi,i+1(i=1,2)为第i和i+1 层之间的水流，m/s；Q1为上边界进入土壤表层的水流，m/s；Q3为重力排水，m/s；Eg为土壤表层的蒸发速率，m/s；Ect为通过气孔的植被蒸腾速率(m/s)[7]。

2.2　参数设置

根据试验地条件，将SiB2中的植被类型设定为9(农田/C3草地)，土壤类型设定为3(壤土)[9]，在实地测量的基础上确定了在花生幼苗期、开花期和饱果成熟期内低丘红壤区农田小流域适于SiB2运行的所有参数，见表1～表3。

2.3　长波辐射修正

运行SiB2所需要的驱动数据包括向下短波辐射、向下长波辐射、水汽压、气温、风速、降水6个基本量，其中鹰潭市余江县红壤生态试验站自动气象站提供了短波辐射、相对湿度、气温、风速、降水量的逐小时观测资料。

水汽压使用相对湿度换算得到，T为气温。

(7)

向下长波辐射是来自整层大气的辐射，即大气逆辐射。SiB2中可以根据已有驱动数据换算向下长波辐射，但我们选择利用江西省鹰潭市气象局自动气象站的逐小时云量数据，修正brunt公式从而计算得到地面有效长波辐射。

表1　花生不同生长期形态学参数

表2　花生和土壤光学特性参数

表3　试验区作物生理参数及土壤特性参数

(8)

则向下长波辐射为：

(9)

式中：ε为地面比辐射率，5、6、8月的值分别为0.97、0.97、0.975[10]；σ为斯蒂芬玻尔兹曼常数；T、d为地面温度；C为春夏季江西地区云对地面有效辐射的减弱系数(0.59)[13]；N为逐小时云量。

2.4　数据同化方案

数据同化的基本思想是假设模型的模拟结果和观测都有一定的不确定性，通过一定的数据同化方法融合这两种来源的数据可以得到一个更优的模型状态估计结果[16]。传统的卡尔曼滤波算法是一种利用观测数据和预测值(两者都含有噪声和误差)来逼近隐藏变量真实值的算法，模型预测值被视为先验值，观测值和模型预测值一起用于获得隐藏变量的更新估计值，即为滤波过程。由于在利用气象模型进行数值模拟的过程中，模拟误差会不断积累，最终导致时间序列越长，结果越不可靠，解决误差积累的核心问题是如何在模型运行过程中引入更准确的初值作为校正。所以我们引入集合卡尔曼滤波算法，在整个模拟过程中设置多个同化点，每次同化引入观测值进行滤波，得到的更逼近真实值的土壤水分数据作为下一阶段的模型初值。

集合卡尔曼滤波是顺序数据同化算法的典型代表之一，用状态样本的集合来代表模型状态后验概率密度函数，通过这些样本的向前积分计算状态的均值和方差，克服了卡尔曼滤波线性化的缺点[17]，从而可以应用到陆面过程等非线性系统模型中，其计算步骤如下：

(10)

ωi,k～N(0,Qk)

(3)计算k+1时刻的卡尔曼增益矩阵Kk+1：

(11)

(12)

(13)

(14)

(15)

vi,k～N(0,Rk)

(16)

(17)

(18)

(5)进入下一时刻返回步骤(2)。

集合卡尔曼滤波主要包括预报和分析两个过程。首先进行模式状态的预报，接着引入观测资料，然后根据观测资料对模式状态进行重新分析，公式(10)为变量的集合预报，公式(13)、(14)为预测误差的分析，公式(16)为对变量的分析，公式(18)为对分析变量的误差分析。设定集合数为48，同化间隔为6 h，SiB2模型误差方差矩阵(Q)和观测误差方差矩阵(R)分别为[7,18]：

3　结果与分析

3.1　模拟与同化结果

各生育期模拟与同化结果见图1-图3。

图1　幼苗期根域土壤水分模拟与同化结果

图2　开花期根域土壤水分模拟与同化结果

图3　饱果成熟期根域土壤水分模拟与同化结果

幼苗期降雨量多且集中，对花生根域土壤湿度影响较大，短期的大量突发性降水使SiB2对根域土壤湿度估计过高，在最初阶段就积累了大量误差，导致模拟结果整体偏高。同化结果虽然在有降水时也产生了过高估计，但随着时间推移可以不断消减误差累积，在无降水或降水偏少时，同化结果接近观测值。

开花期由于降水减少，同时高温伏旱增多，SiB2对根域土壤湿度估计过低，误差随着时间逐渐增大。同化则取得了很好的效果，整体略低于观测值。

饱果成熟期降水等气象因素对模拟结果的干扰明显减小，模拟和同化结果与观测值都保持了较好的一致性，同化结果更接近观测值。

3.2　误差统计

对试验区幼苗期、开花期、饱果成熟期花生地分别为10、15、20 cm深度的根域土壤水分的模拟结果显示(见表4)，SiB2对饱果成熟期的花生根域土壤湿度模拟精度最好，整体偏低；对开花期的模拟精度次之，整体偏低，其离散程度低于饱果成熟期模拟结果；对幼苗期的模拟结果较差，其整体偏高，离散程度较高。在利用集合卡尔曼滤波算法进行数据同化后，均方根误差分别降至原先的1/4、1/4、1/2；平均误差分别降至原先的1/5、1/5、1/2；相对误差分别降至原先的1/10、1/5、1/2。

表4　花生不同生长期模拟结果误差比较

4　结　语

本文模拟了低丘红壤区花生幼苗期、开花期、饱果成熟期3个生长阶段花生根域深度土壤湿度10 d内的变化情况，并发展了基于集合卡尔曼滤波数据同化方案，得出如下结论：

(1)针对江西低丘红壤区花生地土壤性质、作物生物学特性等进行的实地测量，以及对SiB2长波辐射计算的改进，使得SiB2可以初步适用于对当地幼苗、开花、饱果成熟等关键生长期花生作物根域土壤水分的短期动态变化模拟，在经过数据同化后，3个生长期根域土壤湿度模拟结果均方根误差分别降至原先的1/4、1/4、1/2；平均误差分别降至原先的1/5、1/5、1/2；相对误差分别降至原先的1/10、1/5、1/2。

(2)模型的模拟和同化结果对降水等气象因素都具有较好的响应能力，短期强降水会使得模拟结果偏高，长期的无降水或伏旱会使得模拟结果逐渐变低，类似的极端天气现象如果发生在早期或者持续较长时间，则会使得模拟结果积累大量误差，在进行数据同化后则可以有效地避免这种情况。

(3)在花生不同生长期内，根域土壤水分动态变化都表现了不同的特点：幼苗期花生生长发育迅速，生物学参数不稳定影响模型计算，且由于花生枝干发育不完全，保持水土的能力不强，受气象因素干扰较大，春夏之交降水多发，因而模拟结果精度和一致性较差，进行长时间的模拟不现实；开花期花生生长发育速度变慢，植株枝干基本发育完全，有一定的水土保持能力，受气象因素影响较小，模拟精度以及一致性均好于幼苗期，但夏季多发伏旱，降水缺乏使得模拟结果虽然一致性较高，但整体结果偏低；饱果成熟期花生发育完全，植株枝干生长茂密、叶面积指数最大，具有较强的水土保持能力，不易受到气象因素的影响，模拟精度最佳，一致性次于开花期。所以对幼苗期和开花期的花生根域土壤湿度进行模拟时，有必要进行数据同化。

试验中造成误差的原因可能有，一是实际测量不同生长期花生的形态参数如叶面积指数、冠层覆盖比、平均叶长叶宽等误差较大，尤其是在花生幼苗期内，花生植株发育迅速，叶片、枝干、根系等生长快，最初设定的形态参数随着时间的推移会越来越不准确导致误差偏大；二是小区周围有较浓密高大的树林遮蔽了部分阳光，阻碍了风的影响，引起了误差，突发性降水、伏旱等极端天气对模拟结果的影响显著；三是SiB2模型存在很多局限性，集合卡尔曼滤波的同化性能受模型误差的制约较大，同时，如何确定观测误差、集合的大小、同化频率等对同化结果至关重要。江西低丘红壤区地处亚热带季风区水文环境复杂、气候多变、作物生长发育迅速，因此，在应用SiB2对低丘红壤区农田小流域进行长期的模拟试验时仍然有必要对其进行改进以及发展更加精确高效的区域尺度陆面同化方案。

致谢:南京信息工程大学应用气象学院朱博钧博士生对研究给予了帮助，谨致谢意！

参考文献：

[1]李新,黄春林,车涛,等.中国陆面数据同化系统研究的进展与前瞻[J].自然科学进展,2007,17(2):163-173.

[2]高志球,卞林根,程彦杰,等.利用生物圈模型(SiB2)模拟青藏高原那曲草原近地面层能量收支[J].应用气象学报,2002,13(2):129-141.

[3]赵其国,黄国勤,马艳芹.中国南方红壤生态系统面临的问题及对策[J].生态学报,2013,33(24):7 615-7 622.

[4]景元书,张斌.亚热带花生旱坡地的水分动态特征分析[J].中国农学通报,2006,22(11):467-471.

[5]张川,张伟,陈洪松,等.喀斯特典型坡地旱季表层土壤水分时空变异性[J].生态学报,2015,35(19):6 326-6 334.

[6]景元书,张斌,孙秀芬.坡地对花生气候生产潜力利用率的影响[J].花生学报,2003,32(1):17-20.

[7]黄春林,李新.基于集合卡尔曼滤波的土壤水分同化试验[J].高原气象,2006,25(4):665-671.

[8]Sellers PJ, Randall DA, Collatz GJ, et al. A revised land surface parameterization (sib2) for atmospheric gcms.part i: model formulation[J]. Journal of Climate, 1996,9(4):676-705.

[9]Sellers J P,Tucker J C,Collatz J G, et al. A revised land surface parameterization (sib2) for atmospheric gcms. part ii: the generation of global fields of terrestrial biophysical parameters from satellite data[J]. Journal of Climate, 1996,9(4):706-737.

[10]姚镇海,邱新法,施国萍,等.中国比辐射率空间分布特征分析[J].干旱区地理,2015,38(6):1 234-1 240.

[11]王艺, 朱彬, 刘煜,等. 中国地区近10年地表反照率变化趋势[J]. 气象科技, 2011,39(2):147-155.

[12]姚彤, 张强. 我国北方不同类型下垫面地表反照率特征[J]. 物理学报, 2014,63(8):462-470.

[13]高国栋,陆渝蓉,吴必文,等.地面长波辐射量的计算及云对其影响的初步分析[J].气象科学,1985,(2):77-84.

[14]郭元喜, 龚道溢, 汪文珊,等. 中国东部夏季云量与日气温统计关系[J]. 地理科学, 2013,33(1):104-109.

[15]Nakajima T, Higurashi A, Kawamoto K, et al. Nakajima, T. et al. A possible correlation between satellite-derived cloud and aerosol microphysical parameters. Geophysical. Research. Letters. 28, 1171-1174[J]. Geophysical Research Letters, 2001,28(7):1 171-1 174.

[16]韩旭军, 李新. 非线性滤波方法与陆面数据同化[J]. 地球科学进展, 2008,23(8):813-820.

[17]Anderson B D O, Moore J B, Eslami M. Optimal Filtering[J]. IEEE Transactions on Systems Man & Cybernetics, 1979,12(2):235-236.

[18]李超, 惠建忠, 唐千红,等. EnKF优化土壤湿度方程中参数的性能研究[J]. 自然资源学报, 2015,(12):2 131-2 140.