兰州大学资料同化课件讲解

第7讲 大气资料的四维变分同化方法 1 四维变分同化基本原理 四维同化的概念-利用模式消化吸收多时刻观测， 不断改进预报，优化大气状态的估计。 变分方法是一种实现四维同化的有力工具。 在进行大气资料分析时，我们有两种基本的可用 信息：（1）观测；（2）大气遵循的物理规律。 前面我们在作资料分析中用到过一些简化的物理 约束，四维变分同化利用完整的大气模式来作为 物理约束。 四维变分同化的基本思想是调整初始场，使由此 产生的预报在一定时间区间（同化窗口）内与 观测场距离最小 4DVAR示意图： 按照这样的思想，四维同化变分同化可以表述为极 小化下面的目标泛函： 这里x0=x(0), xt=x(t). xt是由下面的预报模式产生的解 : 离散形式 （n=0 成为三维同化） 4DVAR是微分方程反问题 将已知微分方程和定解条件（初条件，边条件） 求方程的解的问题作为正问题，那末，已知方程的 解（部分解）或解的某种函数反求定解条件或者方 程的一些未知项的问题被称之为微分方程的反问题 。因此，四维变分同化也是一类微分方程的反问题 。 求反问题的解的过程称为反演。我们可将观测y 近似看作预报模式（方程）的解的某种函数，那末 上面表述的四维变分同化就是由观测反演初值的问 题。四维变分同化的一个显著特点是利用了过去时 间的观测资料，而且同化后的场是模式的一个预报 场，不会出现不协调的问题。四维变分同化方法还 有能力从一部分观测变量去反演另外的变量。比如 ，由高度的观测反演风场。 关于反问题的进一步讨论 如果将由“原因”推得“结果”的问题称为正问 题，则由“结果”推求“原因”的问题可称为反问 题。 正问题 z=R(u) 这里算子R 已知，由u求z为正问题。反问题是原 来的已知条件未知（或部分未知），而原问题的解 已知，即由z求u的问题，形式上可写为 u=R-1z 这里R-1是R 的逆算子。我们要研究的反问题一般 是指R-1的显式表达式不可知的情况，只能由u的“ 表现”z间接推求u。 一个简单的例子 扩散输送问题 定解条件： 几类反问题： 待定微分方程中的未知参数的反问题算子识别； 待定初始条件的反问题逆时间过程问题； 待定边界条件的反问题边界控制问题； 待定边界形状的反问题几何反问题。 还有的反问题是几类相混合。 解反问题的主要困难不适定性 解的存在性、唯一性和稳定性不满足 吉洪诺夫的论著不适定问题的解法首 先引入“条件适定”的概念，基本思想是： 放弃求精确解转而求近似解解决了解不存在 的困难；近似解总存在，但不唯一，此时再 加适当约束条件，找出具有稳定性的解来。 广义解，目标泛函 反问题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不再现实。正则化的思想是构造 一个连续的算子去逼近A-1，从而得到稳定的（但是近 似的）解。具体而言，他将求稳定的反问题的解归结 为求下面泛函的极小值： 非负泛函,:观测误差.正则系数，比如罚函 数： 如何定? 给出让 迭代。 2, 最优控制理论 四维变分同化将问题提为一个最优控制问题。 最优控制问题的一般提法： 在系统工程中，一个系统通常可以用n个变量完 全描述清楚。我们以向量X 记这n个变量，并称之 为状态向量。系统的运动方程可以用时间间隔 T1,T2上的状态方程来表示： 其中是n 维状态向量；u是r维（r =, 为内积. 在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 四维变分同化的计算过程 实际问题中，微分方程都被离散化，这时目标函数 为 控制方程： 这里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作一次积分，初值为零，而在到 达观测时刻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 X(1,0) X(2,1)= X(1,0)+ dX(2,0） 需要求 先计算 然后 类似： 伴随码（Adjoint code)生成技术 为了得到目标函数对初值的梯度我们必须积 分伴随模式。伴随模式方程实际是原模式的切线 性方程的伴随方程。伴随模式的解析形式只能作 为理论推导用，实际问题是离散化的，预报程序 中还有些是不能写成解析公式的。要保证相应的 伴随模式严格成立，通常的作法是先根据原模式 计算程序写出切线性模式程序，再直接根据切线 性模式程序一一对应地写出伴随程序。一个预报 模式的程序有上万条语句，首先写出他的切线性 模式程序，然后根据切线性模式程序写出伴随程 序，工作量是巨大的。 这里我们用一些简单例子介绍一般规则。 在预报模式中，一个循环过程一般可以表示 为矩阵与向量的乘 YAX 这里A与X无关（线性），则切线性程序和原来相同， YAX 伴随过程为 X*ATY* 如果A与X有关（非线性），那末切线性程序为 YAX 伴随过程为 X*ATY* 例如 DO I=1,N-1 X(I)= a * Y(I+1)*2 END DO 切线性程序为: DO I=1,N-1 X(I)= a * 2*Y0(I+1)* Y (I+1) END DO 为了以下方便，记a * 2*Y0(I+1)B（I） 显然，Y是输入变量，X是输出变量。 注意，这里Y0是“基本变量”，X和Y是扰动变量 。基本变量是由原模式产生的。 需要区分以下两种情况： （1）只有X是输出变量，Y在以后不在用到（Y不是 输出变量）； （2）X和Y都是输出变量，就是说Y在以后还要到。 对第一种情况，上述循环的矩阵是 (N -1)N 对第二种情况，矩阵是 (2N -1)N 切线性程序，第一种情况是 第二种情况实际是 相应的伴随程序分别为 写出程序是 （1） DO I=1,N-1 Y(I+1)= a * 2*Y0(I+1)*X(I) END DO (在切线性程序中Y后来不再用到的情况) 或者 （2） DO I=1,N-1 Y(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, KK DO I=1, II Y(I,K)=0.0 DO L=1,KK Y(I,K)=Y(I,K)+X(I,L)*FLOAT(K) end do end do end do 假定y是输出，x是输入，X在第3重循环后不再运用 。但是要注意X（I,L）在第3重循环内实际上被 重新运用了(对不同的K, X被重新运用),所以伴 随程序应该写为： DO K=1, KK DO I=1, II X(I,K)=0.0 End do End do （将