基于集合卡尔曼滤波和生物圈模型的陆面数据同化

基于集合卡尔曼滤波和生物圈模型的陆面数据同化

1对土壤水分的数据同化64数据的同源性特征（cdda，foys，简称数据同源性）是指在考虑数据的时空分布、观测场地和背景场误差的基础上，在数值模拟的动态执行过程中整合新的观测数据的方法。陆面数据同化系统(LDAS,LandDataAssimilationSystem)是近年来将四维同化方法应用到地球表层科学和水文学中而迅速发展起来的,它作为独立的领域出现是在1995年之后。当前,陆面数据同化的研究主要对给定的陆面模型和水文模型,采用不同的数据同化算法同化地表观测资料、卫星和雷达数据,优化地表和根区土壤水分、温度、地表能量通量等的估算。李新等\,Pathmathevan等利用模拟退火算法对陆面过程模型(SiB2)进行了土壤水分和地表温度的同化实验;Houser等、Schuurmans等把数据同化算法与水文模型结合,进行了土壤水分和潜热通量的同化试验;Galanlowicz等\,Hoeben等、Reichle等、Walker等也对土壤水分的同化方法进行了研究和比较;Crow等、Kumar、Lakshmi进行了地表温度的同化方法的研究。集合卡尔曼滤波是Evensen等在1994年根据Epstein的随机动态预报理论提出的顺序数据同化算法。它将模型状态预报看成近似随机动态预报,用一个状态总体(设数目为N)去代表随机动态预报中的概率密度函数,通过向前积分,状态总体很容易计算不同时间的概率密度函数所对应的统计特性(如均值与方差)。集合卡尔曼滤波的最大特点是它克服了卡尔曼滤波要求线性化的模型算子和观测算子的缺点。经过近十年的研究和发展,集合卡尔曼滤波算法逐渐发展成熟,也逐步受到数据同化研究人员的重视,不同作者[22,23,24,25,26,27,28]对集合卡尔曼滤波算法的理论和应用进行了广泛探讨。土壤水分是地气相互作用、水文循环等研究的关键变量,它影响地表能量通量、径流、辐射平衡、物质迁移等。土壤水分的准确估计对于研究和理解地球表层生物物理过程起着重要作用。土壤水分在空间和时间上变化非常大,建立稠密的地球表层变量的观测网几乎不现实。此外,陆面过程或水文模型虽然可以模拟土壤水分的连续变化,但随着时间的变化,误差不断累积,导致模拟结果较差。通过建立陆面数据同化系统可以把相关的多源信息(地表观测、卫星、雷达等)融入陆面过程模型或水文模型,在有观测的时间段校正模型对土壤水分的预报值,降低误差的累积,从而可以提高对土壤水分的估算精度。2方法陆面数据同化系统主要由数据同化算法、陆面过程模型、数据(驱动数据、参数、观测数据、输出数据)等构成(图1)。2.1层土壤的水分控制方程陆面过程是指发生在地表和大气之间水分、热量和动量交换的作用过程,包括地面上的热力过程、水文过程和生物过程,地气间的能量和物质交换以及地面以下土壤中的热传导和水热输送过程等。本文采用了发展比较成熟的SiB2模型模拟土壤水分的变化。在SiB2中定义了表层、根区和深层3层土壤。公式(1)~(3)为每层土壤的水分控制方程。∂θ1∂t=1D1[Q1-Q12-Egρw],(1)∂θ2∂t=1D2[Q12-Q23-Ectρw],(2)∂θ3∂t=1D3[Q23-Q3],(3)∂θ1∂t=1D1[Q1−Q12−Egρw],(1)∂θ2∂t=1D2[Q12−Q23−Ectρw],(2)∂θ3∂t=1D3[Q23−Q3],(3)式中θ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,Ect分别是土壤表层的蒸发率(m/s)和通过气孔的植被蒸腾率(m/s)。2.2数据同化算法2.2.1卡尔曼滤波的计算步骤1960年,Kalman等针对随机过程状态估计提出Kalman滤波的思想,它包括预报和分析两个步骤。在预报阶段,根据前一时刻的模式状态生成当前时刻模式状态的预报值。在分析阶段,引入观测数据,利用最小方差估计方法对模式状态进行重新分析。随着模式状态预报的继续进行和新的观测数据的陆续输入,这个过程不断向前推进。在卡尔曼滤波算法中,认为模型误差方差矩阵(Q)和观测误差方差矩阵(R)是已知的。卡尔曼滤波的计算步骤如下:(1)设定初始时刻t0的状态变量Xa0a0和它的误差方差矩阵Pa0a0。(2)根据第k时刻(k=0,1,2,…,)的分析场的状态变量Xakak和背景场误差方差矩阵Pakak计算第k+1时刻预报场的状态变量Xfk+1fk+1和背景场误差方差矩阵Pfk+1fk+1。Xfk+1=ΜXak,(4)Ρfk+1=ΜΡakΜΤ+Q,(5)Xfk+1=MXak,(4)Pfk+1=MPakMT+Q,(5)式中M为线性的模型算子,MT为模型算子的转置,Q为模型误差方差矩阵。(3)计算k+1时刻的卡尔曼增益矩阵Kk+1。Κk+1=Ρfk+1ΗΤ(ΗΡfk+1ΗΤ+R)-1‚(6)式中H为线性的观测算子,HT为观测算子的转置,R为观测误差方差矩阵。(4)计算k+1时刻分析场的状态变量Xak+1和它的误差方差矩阵Pak+1。Xak+1=Xfk+1+Κk+1(Ζk+1-ΗXfk+1),(7)Ρak+1=(Ι-Κk+1Η)Ρfk+1,(8)式中Zk+1为k+1时刻的观测值,I为单位矩阵。(5)进入下一时刻,再返回步骤(2)。如此循环往复地执行预报和分析步骤。从卡尔曼滤波的计算公式我们可以看出,在计算过程中该算法需要线性的模型算子和观测算子。对于复杂的非线性系统(特别是地球系统科学),给出线性的模型算子是非常困难的。模型算子和观测算子线性化问题成为制约卡尔曼滤波广泛应用的瓶颈。2.2.2模型状态变量的计算集合卡尔曼滤波是Evensen在1994年根据Epstein的随机动态预报理论提出的顺序数据同化算法。它利用了蒙特卡罗方法的思想,用符合高斯分布的一组随机变量(设数目为N)去代表随机动态预报中的概率密度函数,通过向前积分,计算下一时刻状态总体的概率密度函数,并得到该时刻的统计特性(如均值与协方差)。在本研究中,使用Houtekamer提出的算法。集合卡尔曼滤波的计算步骤如下:(1)初始化背景场。给定N个符合高斯分布的随机变量Xi(i=1,…,N),在本文中,状态变量为地表、根区、深层土壤体积含水量,Xi=[θ1,θ2,θ3]Τi。(2)计算每个随机变量在第k+1时刻的预报值Xfi,k+1。Xfi,k+1=Μk(Xai,k)+wi,kwi,k~Ν(0,Qk),(9)式中Xfi,k+1为随机变量i在k+1时刻的预报值,Xai,k为随机变量i在k时刻的分析值,Mk(·)为非线性的模型算子,在本文中为SiB2模型;Qk为模型误差方差矩阵,wi,k为期望为0,方差为Qk的高斯白噪声。(3)计算k+1时刻的卡尔曼增益矩阵Kk+1。Kk+1=Pfk+1HT(HPfk+1HT+Rk)-1,(10)Ρfk+1=1Ν-1Ν∑i=1(Xfi,k+1-¯Xfk+1)⋅(Xfi,k+1-¯Xfk+1)Τ(11)Ρfk+1ΗΤ=1Ν-1Ν∑i=1[Xfi,k+1-¯Xfk+1][Η(Xfi,k+1)-Η(¯Xfk+1)]Τ,(12)ΗΡfk+1ΗΤ=1Ν-1Ν∑i=1[Η(Xfi,k+1)-Η(¯Xfk+1]⋅[Η(Xfi,k+1)-Η(¯Xfk+1)]Τ,(13)¯Xfk+1=1ΝΝ∑i=1Xfi,k+1,(14)式中Xfi,k+1为第i个随机变量在k+1时刻的预报值,¯Xfk+1为k+1时刻状态变量预报值的平均值,Rk为观测误差方差矩阵,Hk(·)为观测算子,本文中的状态变量和观测均为3层土壤体积含水量,所以Hk(·)为3×3的单位方阵。(4)计算k+1时刻分析场的状态变量平均值Xak+1和背景场误差方差矩阵Pak+1。Xai,k+1=Xfi,k+1+Kk+1[yk+1-H(xfi,k+1)+vi,k+1]vi,k~N(0,Rk),(15)¯Xak+1=1ΝΝ∑i=1Xai,k+1,(16)Ρak+1=1Ν-1Ν∑i=1(Xai,k+1-¯Xak+1)(Xai,k+1-¯Xak+1)Τ,(17)式中Xai,k+1为第i个随机变量在k+1时刻的分析值,¯Xak+1为k+1时刻状态变量分析值的平均值,vi,k为期望为0,方差为Rk的高斯白噪声。(5)进入下一时刻,返回步骤(2)。从公式(10)~(15)中可以看出,集合卡尔曼滤波在计算卡尔曼增益矩阵和k+1时刻分析场的状态变量时,利用蒙特卡罗方法有效地避免了卡尔曼滤波需要线性化的模型算子和观测算子的难题。2.3驱动数据和观测数据“全球能量与水循环亚洲季风之青藏高原试验”(GAME-Tibet)已经积累了大量观测资料并且已被用于青藏高原的水热平衡模拟研究,本文所需的驱动数据和观测数据均取自GAME-Tibet试验区MS3608观测站1998年7月6日至8月9日的观测数据,时间分辨率为1h。该站点位于31°13.6′N,91°47.0′E,海拔高度4610m,地表覆盖类型为矮草。2.3.1辐射和降水条件SiB2模型所需的驱动数据包括短波辐射、气温、风速、水汽压、长波辐射和降水。其中短波辐射、气温、风速资料直接来自于自动气象站(AWS)观测,长波辐射使用SiB2中的标准方法-Brunt公式计算,水汽压由湿度换算,降水资料取自自动记录的雨量筒。2.3.2errs等和李新等的定义青藏高原MS3608站点的土壤、植被等参数主要参照Sellers等和李新等的定义。不同土壤的饱和水力传导系数变化范围很大(1.2×10-4~6×10-6m/s),对于模拟结果影响很大,本文采用6×10-5m/s。2.3.3深度4.10,5.1土壤水分数据由时域反射仪(TDR)测量,深度分别为4,20,60,100,160和196cm。在同化实验中只利用了4,20,100cm的土壤水分观测数据。2.4otmoanstraft安氏原理对于模拟、同化、观测结果,我们采用了两种误差分析方法:(1)均方根误差(RMSE,RootMeanSquareError);(2)平均误差(AVE,AverageError)。RΜSE=√1ΝΝ∑t=1(Οbst-Xt)2,(18)AVE=1ΝΝ∑t=1(Οbst-Xt),(19)式中N为整个观测期,Obst为t时刻的观测值,Xt为t时刻的SiB2模拟值或同化结果。3结果分析3.1亚风场模拟结果有6个月的情况,自然集合数取48,利用1998年7月6日、7日逐时4,20,100cm的土壤水分数据作为初始状态变量的集合,同化间隔为6h(利用00,06,12和18时的4,20,100cm土壤体积含水量观测数据同化SiB2模型输出的表层、根区、深层土壤水分),模型误差方差矩阵(Q)和观测场误差方差矩阵(R)分别为Qi=[2.5×10-5,0,00,1.44×10-5,00,0,2.5×10-5],Rk=[1×10-4,0,00,4×10-4,00,0,1×10-4].图2a~c是1998年7月6日至8月10日亚洲季风实验青藏高原实验区MS3608站点表层、根区、深层土壤水分同化结果。图2a是表层土壤水分的同化结果。由图可看到,降水对于表层土壤水分的影响很大。当有降水时,SiB2模拟结果与4cm观测值相比有些偏大;无降水时,SiB2模拟结果接近观测值。对于同化结果,有降水时,同化降低了SiB2对土壤水分过高的估计,使结果接近观测值;无降水时,同化结果与观测相比略偏低,在部分时段同化效果比SiB2直接模拟的效果还要差。在整体上,同化还是取得了较好的效果,但扰动较大。图2b是根区土壤水分的同化结果。SiB2对根区土壤水分的模拟结果整体偏大,降水对根区土壤水分的影响仍然比较大。有大的降水时,模拟的根区土壤水分急剧增大。在根区,同化对土壤水分的估计有比较明显的效果。图2c是深层土壤水分的同化结果。降水对深层土壤水分影响很小,土壤水分的变化幅度不大。经历一次较强的降水后,深层土壤水分增大的幅度较小。SiB2模拟结果略低于100cm的观测值,同化效果显著,同化结果几乎和观测值重合。3.2同化土壤水分模型图3为表层、根区、深层土壤水分同化后误差方差随时间的变化。表层、根区、深层土壤水分方差都比较稳定,在有观测时,三层土壤水分的误差明显降低,所以同化观测可以消除误差的累积效应,抑制误差的增长。表层、根区、深层误差方差分别为<5×10-5,1.5×10-4,2.5×10-5。由表1可以看到,在SiB2模型运行过程中同化土壤水分的观测值降低了土壤水分模拟值的误差。表层、根区、深层土壤水分模拟值的均方根误差约分别降低为原来的1/2及1/9,平均误差也分别降低为原来的1/4及1/17。在同化试验中,考虑到土壤水分观测数据的代表性问题,观测误差给定的较大,所以同化前后相比,表层、根区的误差减少幅度不大。4应用结果分析我们发展了基于集合卡尔曼滤波和SiB2的试验型单点陆面数据同化方案,利用GAME-Tibet土壤水分观测数据进行了同化试验,得到以下结论:(1)基于集合卡尔曼滤波的数据同化方法可以降低模拟的误差,同化后和同化前相比,表层、根区、深层土壤水分模拟值的均方根误差分别降低为原来的1