同化理论的发展及其在遥感方面的应用

同化理论发展及其在遥感方面的应用
摘要：在为数值预报模式提供准确、合理初值问题上,资料同化是一种行之有效的方法。

其基本含义是根据一定的优化标准和方法,将不同空间、不同时间、采用不同观测手段获得的观测数据与数学模型有机结合,纳入统一的分析与预报系统,建立模型与数据相互协调的优化关系,使分析结果的误差达到最小。

其中应用最为广泛的同化方法是变分法和卡尔曼滤波法。

另外遥感信息与作物生长模型结合来进行作物监测和产量预测，已经逐渐成为一种接受度较大、应用较为广泛的方法之一。

关键字：同化 变分 卡尔曼滤波 遥感
1．同化的概念
在为数值预报模式提供准确、合理初值问题上,资料同化是一种行之有效的方法[1]。

它是由早期气象学中的分析技术发展起来的[2-3]。

其基本含义是根据一定的优化标准和方法,将不同空间、不同时间、采用不同观测手段获得的观测数据与数学模型有机结合,纳入统一的分析与预报系统,建立模型与数据相互协调的优化关系,使分析结果的误差达到最小。

一般一个资料同化系统包括观测数据集、动力模型和数据同化方案三部分。

以模式的一种初估状态或其他一些不重要的状态为初始场,由模式得到的解常称之为背景场;结合观测数据集,通过同化过程产生能够相对准确反映真实状态的一种最优估计,称之为分析场。

一般而言,分析场是背景场和观测场的加权平均,其方差始终比观测场和背景场的方差小。

2．同化方法的发展
2.1逐步订正法
1954年,Gilchrist 等提出了理想的逐步订正法。

其原理是从每一个观测中减去背景场得到观测增量,通过分析观测增量得到分析增量,然后将分析增量加到背景场上得到最终的分析场。

每一个分析格点上的分析增量通过其周围影响区域内观测增量的线性组合而加权,观测权重与观测位置和格点之间的距离成反比。

该方法的表达式可写为:
11(,)[()()]
()()(,
)n b i a b n i w i j y i X i X j X j w i j ==-=+∑∑
()b X i 为插值到观测点i 上的背景场信息；y(i)为相应的观测值；()a X i 为格点j 的订正值；
权函数22,22
,(,)max 0,i j i j R d w i j R d ⎛⎫-= ⎪ ⎪+⎝⎭
，其中,i j d 为格点i ，j 之间的距离，R 是给定的影响半径，w(i,j)随距离的增加而递减，即分析场是观测场和背景场加权平均的结果。

2.2最优插值法
1963年Gandin 提出了最优插值法(简称OI),资料同化才有了基于统计估计理论的基础。

从统计意义上来说,它是一种均方差最小的线性插值方法。

相比逐步订正法而言,这里最大的改进就是权重考虑了背景场和观测误差的统计特征,故又称之为统计插值法[4]。

这种方法在20世纪80—90年代的数值天气预报业务中得到了广泛应用。

它是通过简化算法给定最佳线性无偏估计方程中的权矩阵。

其计算公式为:
([])a b b X X K y H X =+-
a X 为模型的分析场，
b X 为模型的背景场，y 为观测向量，H 为观测算子。

分析场的增益
矩阵1()T T K BH HBH O -=+，通常也称为权矩阵。

对于最优插值法,关键就是如何求K 。

此方法的基本假设是,对于一个模型变量,在确定它的增量时,只有几个观测值是重要的。

如果在观测值的选择方面给出正确假设,最优插值法就可以很简单地实施,需要计算量相对较小。

它的缺点是不能适用于观测算子较复杂的系统,当对模型的不同部分采用不同的观测数据集时,会使分析场产生伪噪声,并且无法确保大尺度与小尺度分析之间的一致[5]。

2.3变分法
2.3.1三维变分
三维变分(简称3DV AR)不仅可以用来处理观测算子为非线性的情况，而且也避免了求OI 中的增益矩阵K 。

它是通过迭代一个代价函数，产生一个分析时刻的真实状态的最优估计，通常目标泛函可表示为:
1111()()()[()]()[()]22
T T b o b b o o J x J J X X B X X y H x E F y H x --=+=--+-+- 式中，x 为模式状态变量构成的分析变量；b X 为x 的先验估计，常称背景场向量；B 为
背景误差协方差矩阵；E 、F 分别代表观测误差和代表性误差，代表性误差是一个观测算子在从模式空间向观测空间转换时产生的不精确性的一种估计，将观测算子近似为线性算子就可产生这种误差；H 是观测算子，它将模式格点上的分析变量通过空间插值和物理变换到观测空间与观测相当的量，方便与观测量比较；o y 为观测值；T 代表矩阵的转置；-1代表矩阵的
逆。

变分问题就是求上式取得极小时的分析状态的迭代解，这个解代表了真实状态的一个最大似然估计(最小方差)。

3DV AR 是在某一时刻进行的分析，前一时刻的同化结果可作为后一时刻模式运行的初始场。

在实现过程中，一个主要困难是如何合理设计背景场误差协方差矩阵B 的模型，不仅需要正确定义所有模型变量之间的背景协方差,而且要求各变量的误差合理，因为这会决定在分析场中观测场和背景场的权重，另外背景误差协方差矩阵必须是正定的。

该方法有很多优点。

首先，不同OI 中的局部分析，3DV AR 进行的是三维空间的全局分析，避免了分析不是全局最优的问题；其次,观测算子可以是非线性的；第三，它可以将一些额外的项，加到目标泛函里来实现一些外部的弱约束，如平衡项或平滑项；第四，计算中只需要求切线性算子及其伴随算子，等等。

但在使用3DA VR 时，无法用后面时刻的资料来订正前面的结果，同化的解在时间上不连续，这是它的不足之处。

2.3.2四维变分
为了在整个时间序列的分析过程中保持变量的动力协调，四维变分资料同化的基本思想被提出。

四维变分(简称4DV AR)是3DV AR 的推广，它可以使用时间上分布的观测资料，是在一个时间段上做同化，它和3DV AR 所采用的模型状态方程和观测投影算子都是相同的。

其目标函数可定义为：
111
11()()()([])([])22n T T b o b b i i i i i i i i J x J J X X B X X y H x O y H x --==+=--+--∑ 其中，下标i 表示给定的观测时刻，分布于时段[0,]N T ；i O 为第i 时刻的观测误差协方差
矩阵；B 为背景误差协方差矩阵,它由初始时刻指定；i H 为第i 时刻的观测算子；i y 、i X 分别为第i 时刻的观测和状态预报量，且i X 满足方程0(),(0,1,)i i X M x i n →==，0i M →为从起始时刻0
预报到i 时刻的模式预报算子。

类似3DA VR,4DA VR 的最终目的仍是求使得目标函数达到最小时所对应的a X ，就是要求的分析场，使得以它为初值的模式运行结果与这段时间内对应观
测值时刻的观测的偏差之和达到最小。

4DV AR 的核心部分就是伴随模型的建立，通常可有两种途径实现：第一种是伴随的差分，即由连续的正模型得到连续的伴随模型方程和连续的目标泛函的梯度表达式，再对连续的伴随模型和目标泛函的梯度进行差分离散；第二种途径是差分的伴随,它是由离散的正模型直接导出离散的伴随模型及目标泛函的梯度，但国际上一直认可采用的大多是后一种途径。

4DV AR 有很多优点：可同化多时刻的资料(当前观测可对过去的分析结果产生影响)；B 矩阵在同化窗口是隐式发展的，可以考虑一点flow-dependent;另外可在代价函数中加上其他的弱约束项。

其不足之处在于，目前为了减小代价虽然有了一些降阶和“自动伴随编译器”等简化办法，但对于复杂模式而言，由正模型的离散形式直接得到伴随方程的离散形式，仍然是一个比较艰巨的任务；另外，在计算方面，4DV AR 比3DV AR 代价大得多。

2.4卡尔曼滤波法
2.4.1线性卡尔曼滤波
KF 以分析误差的最小方差为最优标准,在假定系统是线性的、噪音是白色、高斯型条件下的一种递归资料处理方法。

与4DV AR 寻求整个同化时段的最优解不同,它只着眼于求解观测时刻的最优分析值,在给出这一时刻的最优分析值的同时,也给出了分析误差的分布,这是变分方法所不具备的。

其基本思想是,首先进行模式状态的预报,然后引入观测数据,根据观测数据对模式状态进行重新分析,接着再进行预报,进而完成预报、分析、再预报的循环过程。

若将背景场误差协方差矩阵B 和分析场误差协方差矩阵A(即a t X X -)分别以f P 和a P 代替,令M 为模型预
报算子,其模型误差的协方差矩阵为Q,在认为模型误差与观测误差彼此不相关时,对于线性系统而言,H 、M 都是线性的,假设误差符合正态分布且增长是线性的,该方法的计算可分为预报和分析两个部分。

其计算流程如下图所示。

在分析部分中, ()K i 的计算需要前步误差()f P i ,分析值()a X i 是由()K i 和()f X i 计算得到, ()a P i 的计算是为了计算下一步误差(1)f P i +。

它与4DV AR 相比优点在于预报误差随模式动力发展而发展,f P 矩阵flow-dependence; KF
显式发展了背景误差协方差矩阵,而4DV AR 只能隐式发展背景误差协方差矩阵,且KF 不用写模式的伴随,因而理论上它比依赖于伴随模式的反向积分的4DV AR 能更方便地应用于真实大气中。

但具体实施起来却有相当的难度[6],因为背景误差协方差矩阵的及时传播、计算需要2N 倍的模式积分时间,而且还需要存储N ×N 维B 矩阵的空间,Q 的估计也很难给出。

这些都意味着要付出昂贵的计算时间与占用大量的机器内存。

2.4.2 集合卡尔曼滤波
实际的大气和海洋预报模式中大多是高维的非线性系统,因此要直接应用广义卡尔曼滤波仍然有难度,主要困难在于背景误差协方差矩阵的预报。

但随着集合预报思想的发展,并逐渐和卡尔曼滤波的研究相结合,形成了集合卡尔曼滤波的方法(简称EnKF),很好地克服了上述困难。

它是用集合样本来估计背景误差协方差矩阵,解决了误差协方差矩阵预报中遇到的问题。

其最
大的优点就是,用集合的思想解决了背景误差协方差矩阵的估计和预报的困难,同时也解决了
卡尔曼滤波中非线性近似的应用问题,也避免了4DV AR中伴随模式的使用。

值得注意的是,这种方法是在假设背景误差的概率密度函数服从高斯分布的前提下,而实际中往往不满足,所以
这是应用中存在的一个问题,同时它也会存在滤波发散和不平衡问题。

3．同化在遥感上的应用
遥感技术结合地理信息系统和全球定位系统等其它现代高新技术，可以实现信息收集和分析的定时、定量、定位，它具有客观性强，不受人为干扰的特点。

我国运用这一技术开展农业监测，可以提高农业决策的科学化水平，进而为农业生产提供高质量的服务。

其中农作物遥感监测是遥感应用的一个重要主题，而卫星遥感技术具有快速、宏观、准确、客观、及时、动态等特点，用于大范围作物长势监测和产量预测具有得天独厚的优势[7]。

但是，由于遥感监测常常受卫星遥感数据空间分辨率、时间分辨率等因素的影响，且遥感信息大多反映的只是地表或作物冠层瞬间的物理状况，而对于作物生长发育、产量形成、以及个体与群体之间生长状况的关系等，目前通过遥感手段还难以揭示。

作物生长模型能够利用环境因素模拟作物生长过程，揭示作物生长和发育的原因与本质，已成为农业研究最有力的工具之一。

因此，遥感信息与作物生长模型结合来进行作物监测和产量预测，已经逐渐成为一种接受度较大、应用较为广泛的方法之一。

利用二者各有的优势，如何将二者结合，发挥各自的优势是最近几年来国内外的研究热点。

作物生长模型与遥感数据的同化是指通过调整作物生长模型中与作物生长发育和产量形成密切相关的、其它方法难以获得的初始条件或参数的值来缩小某些遥感“观测值”与相应模型模拟值之间的差距，从而达到估计这些初始值或参数值的目的。

在遥感数据与作物生长模型的同化过程中，既可以利用某个/某些作物生长状态变量（如LAI）作为同化变量，也可以利用作物冠层反射率等作为同化变量。

4．总结与展望
综上所述，我们回顾了同化理论的发展历史并同时列举了几种应用最为广泛的同化方法，最后略述了同化方法在遥感方面的应用。

我们可以看到同化方法在经过了长时间的发展之后已经日渐成熟，被应用于多个领域，并且我相信，随着进一步的理论结合实际，它将会被不断改进、不断优化，最终达到人们所期望的程度。

参考文献：
1.官元红,周广庆；资料同化方法的理论发展及应用综述；气象与减灾研究；2007年12月第30卷
第4期。

2.丑纪范；四维同化的理论和新方法[M]∥数值天气预报中的若干新技术.北京:气象出版
社,1995:262-294.
3.Bouttier F, Courtier P. Data assimilation concepts and methods [M]∥ECMWF. Meteorological
training course lecture Series.2002:1-59.
4.韩桂军；海洋数据同化方法研究[R].北京:中国科学院大气物理研究所,2004.
5.Lorenc A C.A global three-dimensional multivariate statistical interpolation scheme
[J].Mon Wea Rev,1981,109 (4):701-721.
6.高山红,吴增茂,谢红琴；Kalman滤波在气象数据同化中的发展与应用[J].地球科学进
展,2000,15(5):571-582.
7.杨邦杰, 裴志远, 周清波, 等. 我国农情遥感监测关键技术研究进展[J]. 农业工程学报, 2002,
18(3): 191-194。