《兰州大学资料同化课件讲解》由会员分享，可在线阅读，更多相关《兰州大学资料同化课件讲解（53页珍藏版）》请在金锄头文库上搜索。
1、第7讲 大气资料的四维变分同化方法,1 四维变分同化基本原理四维同化的概念-利用模式消化吸收多时刻观测，不断改进预报，优化大气状态的估计。变分方法是一种实现四维同化的有力工具。 在进行大气资料分析时，我们有两种基本的可用信息：（1）观测；（2）大气遵循的物理规律。前面我们在作资料分析中用到过一些简化的物理约束，四维变分同化利用完整的大气模式来作为物理约束。 四维变分同化的基本思想是调整初始场，使由此产生的预报在一定时间区间（同化窗口）内与观测场距离最小,4DVAR示意图：,按照这样的思想，四维同化变分同化可以表述为极小化下面的目标泛函：这里x0=x(0), xt=x(t). xt是由下面的预报
2、模式产生的解:离散形式（n=0 成为三维同化）,4DVAR是微分方程反问题将已知微分方程和定解条件（初条件，边条件）求方程的解的问题作为正问题，那末，已知方程的解（部分解）或解的某种函数反求定解条件或者方程的一些未知项的问题被称之为微分方程的反问题。因此，四维变分同化也是一类微分方程的反问题。求反问题的解的过程称为反演。我们可将观测y近似看作预报模式（方程）的解的某种函数，那末上面表述的四维变分同化就是由观测反演初值的问题。四维变分同化的一个显著特点是利用了过去时间的观测资料，而且同化后的场是模式的一个预报场，不会出现不协调的问题。四维变分同化方法还有能力从一部分观测变量去反演另外的变量。比如
3、，由高度的观测反演风场。,关于反问题的进一步讨论如果将由“原因”推得“结果”的问题称为正问题，则由“结果”推求“原因”的问题可称为反问题。正问题 z=R(u) 这里算子R 已知，由u求z为正问题。反问题是原来的已知条件未知（或部分未知），而原问题的解已知，即由z求u的问题，形式上可写为u=R-1z 这里R-1是R 的逆算子。我们要研究的反问题一般是指R-1的显式表达式不可知的情况，只能由u的“表现”z间接推求u。,一个简单的例子扩散输送问题定解条件：几类反问题：待定微分方程中的未知参数的反问题算子识别；待定初始条件的反问题逆时间过程问题；待定边界条件的反问题边界控制问题；待定边界形状的反问题几
4、何反问题。还有的反问题是几类相混合。,解反问题的主要困难不适定性解的存在性、唯一性和稳定性不满足吉洪诺夫的论著不适定问题的解法首先引入“条件适定”的概念，基本思想是：放弃求精确解转而求近似解解决了解不存在的困难；近似解总存在，但不唯一，此时再加适当约束条件，找出具有稳定性的解来。,广义解，目标泛函反问题Au =z 的广义解：uU， zZ对于给定的zZ在集U上使 取极小值的u*U, 称为方程Au=z在U上的广义解。若zAU，广义解等于经典解。(AU为U的映象), 为距离。 经典解解不存在: z不属于AU,假定算子方程 Au =z 的逆算子A-1存在但不连续依赖于z，由u=A-1z计算u不再现实。
5、正则化的思想是构造一个连续的算子去逼近A-1，从而得到稳定的（但是近似的）解。具体而言，他将求稳定的反问题的解归结为求下面泛函的极小值： 非负泛函,:观测误差.正则系数，比如罚函数：,如何定? 给出让 迭代。,2, 最优控制理论,四维变分同化将问题提为一个最优控制问题。 最优控制问题的一般提法：在系统工程中，一个系统通常可以用n个变量完全描述清楚。我们以向量X 记这n个变量，并称之为状态向量。系统的运动方程可以用时间间隔T1,T2上的状态方程来表示：,其中是n 维状态向量；u是r维（r n）控制向量，它是从系统之外按一定要求施加到系统上来的；f是n 维向量函数。由方程可以看出，只要f 满足一定
6、的条件，在确定的初始状态下，如果在时间间隔t0, TT1, T2上给定了一个控制规律uu（t）, 那末状态方程在t0,T上将有唯一解。这个解表示了系统的相点在n维状态空间的一个运动,控制规律u不同时，相应的系统运动也是不同的所谓最优控制问题就是要选择适当的控制规律 u(t) 使相应的品质指标(cost function, 代价函数或目标函数)达到最小这里第一项代表了对终端条件的要求，第二项代表了对整个变化过程的要求，而总的品质指标是这两方面要求的综合,最优控制问题的一般提法应为：给定了系统的状态方程以及品质指标和初始条件（t0, x0），目标集，控制域后，要寻找一个容许控制u*（t），使系统在
7、时间间隔t0, T上由初始状态X0转移到目标集上的某一点Xr，且使相应的品质指标J（u）为极小这里终端时间T可以是固定的，也可以是自由的 (容许控制的作用下，系统由初始状态X0转移到终端状态X有些情况，我们要求终端状态应该是状态空间中某一个点集中的一个点集合常称为目标集当只是状态空间中的一个确定的点时，称终端是固定的当是整个状态空间时，称终端是自由的),根据这些要求所求得的控制函数称为最优控制或极值控制，记为u*（t）而相应的状态向量称为最优轨迹或极值值轨迹，记为X*(t）,从本质上讲，所求得的最优控制u*是时间的函数但从形式上来看，它可以直接表示为时间的函数u*（t），也可以表示为状态的函数
8、u*（x）如果求得的最优控制被直接表示为时间的函数，即u*u*（t），则称所求得的是开环最优控制若求出的最优控制被表示为状态的函数，即 u*u*(x）则称所求得的是闭环最优控制从控制理论的一般知识可以知道，闭环控制比开环控制有不少优越之处，所以在实际应用中，总是力图寻找闭环控制规律但是，在实际计算中，寻找闭环最优控制规律要比寻找开环最优控制规律困难得多,这里所谓“最优”的标准就是品质指标J 取极小。但J 的形式的确定完全取决于对本身的具体要求以及工作人员的工作经验，并无普遍适用的法则。可以说，许多最优控制问题的最困难之点就在于无法提出一个恰当的品质指标泛函。,3, 用变分法解最优控制问题,显然
9、最优控制问题是一个求带约束的泛函的极小问题，前面讲述的变分方法同样适用于解决这一问题。以下用一个例子来说明一般做法。,反演初值和k(x), 定义目标泛函：,无约束极小：,变分：,分部积分得到：,利用了边界：,要求终端和边界：,要求变分为零得到伴随方程（它的定解条件前面已经给出）两点边值问题 ： （迭代求解）,用下降算法,梯度：,这里x0=x(0), xt=x(t)，xt是由下面的预报模式产生的解:控制变量也可以是方程的参数，比如s(x),回到四维变分同化来，也是一个最优控制问题极小化下面的目标泛函：,从另一角度来推导： 先对伴随算子做些说明。在希尔伯特（Hilbert）空间（H空间），称A*为
10、线性算子A 的伴随算子，只要对于属于H的任意元素x，y均有 =, 为内积. 在Hilbert空间，内积的定义： (在N 维向量空间内积的定义显然，N 维向量空间的矩阵算子A的伴随算子就是其转置AT ，因为xTAy= yT ATx),预报模式方程：线性化得到切线性模式: 伴随模式:,（7.2.1),(7.2.2),（7.2.3),F*是 的伴随算子(转置）：将u*和(7.2.2)作内积，减去u和(7.2.3)作内积,然后在(0,)区间对时间积分得到（空间边界变分为零条件下）:,算子 的伴随算子,如果目标函数是：让得到表明J对初值的梯度是 。,4 四维变分同化的计算过程 实际问题中，微分方程都被离
11、散化，这时目标函数为控制方程： 这里n是有观测的总时间层数，将目标函数)写为,（7.4.1),为了获得x0的最优估计x*0, 需计算J 相对x0的梯度：其中的计算需要对伴随模式积分。由于x0的扰动x0引起的J0的扰动由（7.4.1）略去高阶小量,由上面两式得到,（9.4.6）,对预报模式 作扰动得到切线性模式： 代入(7.4.6)中的 ,(7.4.6) 取极限 ，(7.4.6)意味着(7.4.7) 记为,即下面 方程的解：,(7.4.9),(7.4.8),即每一个 是由时间 出发，以 为初值向后积分伴随方程到t=t0得到。 由于方程（7.4.8-9）都是线性的，实际上只须对伴随模式由tr到t0
12、作一次积分，初值为零，而在到达观测时刻tr时，伴随变量加上如(7.4.9)右端表示的一个强迫项。,实际过程,举例： 预报,X(1,0)=a， X(2,0)=b，X(1,1)= X(1,0)* X(2,0)+c X(1,0) X(2,1)= X(1,0) +d X(2,0)X(1,2)= X(1,1)* X(2,1)+c X(1,0) X(2,2)= X(1,1) +d X(2,1)目标函数： J=0.5（X(1,2)-y1）*2+( X(2,2) y2)*2,线性方程,X(1,0)= a， X(2,0)= b， X(1,1)= X(2,0)* X(1,0)+ X(1,0)* X(2,0)+c
13、X(1,0) X(2,1)= X(1,0)+ dX(2,0）,需要求,先计算,然后,类似：,伴随码（Adjoint code)生成技术,为了得到目标函数对初值的梯度我们必须积分伴随模式。伴随模式方程实际是原模式的切线性方程的伴随方程。伴随模式的解析形式只能作为理论推导用，实际问题是离散化的，预报程序中还有些是不能写成解析公式的。要保证相应的伴随模式严格成立，通常的作法是先根据原模式计算程序写出切线性模式程序，再直接根据切线性模式程序一一对应地写出伴随程序。一个预报模式的程序有上万条语句，首先写出他的切线性模式程序，然后根据切线性模式程序写出伴随程序，工作量是巨大的。,这里我们用一些简单例子介绍
14、一般规则。 在预报模式中，一个循环过程一般可以表示 为矩阵与向量的乘 YAX 这里A与X无关（线性），则切线性程序和原来相同， YAX 伴随过程为 X*ATY* 如果A与X有关（非线性），那末切线性程序为 YAX 伴随过程为 X*ATY*,例如DO I=1,N-1X(I)= a * Y(I+1)*2END DO切线性程序为:DO I=1,N-1X(I)= a * 2*Y0(I+1)* Y (I+1)END DO 为了以下方便，记a * 2*Y0(I+1)B（I）显然，Y是输入变量，X是输出变量。 注意，这里Y0是“基本变量”，X和Y是扰动变量。基本变量是由原模式产生的。,需要区分以下两种情况：
15、 （1）只有X是输出变量，Y在以后不在用到（Y不是输出变量）； （2）X和Y都是输出变量，就是说Y在以后还要到。对第一种情况，上述循环的矩阵是,(N -1)N,对第二种情况，矩阵是 (2N -1)N,切线性程序，第一种情况是第二种情况实际是相应的伴随程序分别为,写出程序是 （1） DO I=1,N-1Y(I+1)= a * 2*Y0(I+1)*X(I) END DO (在切线性程序中Y后来不再用到的情况) 或者 （2） DO I=1,N-1Y(I+1)= a * 2*Y0(I+1)*X(I)Y(I+1) END DO （在切线性程序中Y后来还要用到的情况）,分析上面的式子我们不难找到规律。 一般如果切线性程序是 X=aYbZ, 这里a, b是系数，Y和Z是输入，那末伴随程序可以容易写出： YaX, z=bX。或者YaXY, Z=bX+Z (Y, Z在切线性程序中后来还要用到的情况)。,下面的求和过程（线性程序）DO K=1, KKDO I=1, IIY(I,K)=0.0DO L=1,KK Y(I,K)=Y(I,K)+X(I,L)*FLOAT(K)end doend doend do假定y是输出，x是输入，X在第3重循环后不再运用。但是要注意X（I,L）在第3重循环内实际上被重新运用了(对不同的K, X被重新运用),所以伴随程序应该写为：,
