自动生成四维变分同化共轭码的软件方法

[收稿日期]2007 05 18

[作者简介]潘纪明(1960 ),男,江苏南京人,南京审计学院信息科学学院工程师,主要研究方向为网络工程和软件开发;吴令云(1979 ),女,江苏南京人,南京审计学院信息科学学院教师,主要研究方向为计算机网络。

第4卷 第4期2007年11月

南京审计学院学报

Journal o f N anjing A udit U niversity

V ol.4,No.4Nov.,2007

计算机科学与技术

自动生成四维变分同化共轭码的软件方法

潘纪明,吴令云

(南京审计学院信息科学学院,江苏南京 210029)

[摘 要]四维同化(4 D assimilatio n)是当前气象(meteo ro log y)科学研究的重要方向,计算四维同化需要编写大量共轭码(adjo int code),手工编写共轭码非常困难,本文介绍了笔者设计的自动编写共轭码语句的软件代码及其原理和算法,实验证明,该软件有很好的效果。

[关键词]气象;四维同化;共轭码;软件

[中图分类号]T P 311.52 [文献标识码]A [文章编号]1672 8750(2007)04 0076 03

四维变分同化简称为四维同化,是气象预报的重要方法,它能把观测数据与大气动力学模式结合起来,

目前已有一些国家开始运用这一方法,它也成为目前气象研究的重要课题。四维同化问题必须用到预报方程,同时要由预报方程产生共轭码以实现最优化。气象科学一般都用Fo rtran 语言编写软件。四维同化用到许多预报方程,它们是微分方程,四维同化软件中把这些微分方程化为差分形式,一般用Fortran 语言写成;它的共轭码也是用Fortran 语句写成的。共轭码是四维同化软件的重要部分,而共轭码行数很多,数以万计,稍有不慎就会出错,这就有必要编写一个由预报方程写出共轭码程序的软件。现将笔者开发的写共轭码的软件介绍如下:

一、四维同化基本原理

本文仅简介四维变分同化的基本原理,实际情况可能更复杂,但基本原理是一致的。从数学观点来看,四维变分同化就是求代价函数的条件极值,以下简介代价函数及其条件极值。

设所考虑的气象问题涉及p 个气象要素,诸如横向风速、纵向风速、温度、湿度等,称它们为大气控制变量。在所考虑区域D 上取一些格点x k ,建立代价函数:

J =1/2

N i=0

m

k=1

[U(t i

,

x k )-U O

(t i ,x k )]T

W ik [U(t i ,x k )-U O

(t i ,x k )]+1/2 m

k=1

[U(t 0,x k )-U kb ]T W bk [U(t 0,x k )-U kb ]

(1.1)

其中U(t,x k )是在区域D 的格点x k 上大气控制变量所成p 维向量,t i 是观测时间,U O

(t i ,x k )是U(t i ,x k )的观测值,W ik ,W bk 是权矩阵,U kb 表示对初始时刻U(t 0,x k )的背景估计。

约束条件为动力学方程:

U t

=N (U)(1.2)

!

76!

(1.2)一般化为差分形式用Fo rtr an语句写成。

(1.2)扰动得切线方程:

U

t=N∀(U) U(1.3)其中N∀(U)是N(U)的导算子。

四维同化就是在约束条件(1.2)下使(1.1)最小而求得初始场U(t0,x k)的估计值。为求该条件极值,总采用优化方法,一般都要求出(1.1)中J的梯度。由于J的自变量太多,往往数以千计,求梯度的一般方法不适用。常用的一个重要方法是共轭码方法,其数学原理如下:

设D是a#b矩形区域网格,观测时间是等间格的,t i+1-t i=⊿t,用向量U i表示在时刻t i,所有格点上控制变量所成向量:

U i=U(t i,x1)∃

U(t i,x m)

,

其中m=ab,U i是abp维向量,再令U O i表示U i在时刻t i,观测值i=0,1,2,...N,

其中,U O i=U O(t i,x1)

∃

U O(t i,x m)

。

(1.2)差分为:

(U i+1-U i)/t=f(U i),

上式化为:

U i+1=U i+tf i(U i) i=0,1,2,......N-1(1.4)将(1.4)扰动得:

U i+1=(I+M i) U i i=0,1,...N-1(1.5)其中M i是矩阵:

M i=t u i+1,1

u i,1∃∃

u i+1,1

u i,a bp ∃∃∃∃

∃∃∃∃ u i+1,a bp

u i,1∃∃

u i+1,abp

u i,a bp

(1.5)迭代可得:

U i+1=M i M i-1∃M0 U0(1.6) (1.1)可写为:

J=1/2

N

i=0

(U i-U O i)T W i(U i-U O i)+1/2(U0-U b)T W b(U O-U b)

其中U b表示背景估计。U b=u1b ∃∃u mb

W i,W b是权矩阵,将它扰动得

J=( U0)T W b(U0-U b)+

N

i=0

( U i)T W i(U i-U O i)

将(1.5)代入得

J=( U0)T W b(U0-U b)+

N

i=0

( U0)T M T0∃M T i-1W i(U i-U O i)

!

77

!