1、第二章 资料处理和客观分析2.1 资料处理随着气象观测手段的发展和现代化，得到的资料数量增多，门类拓广，有常规站点观测，有非定点海洋观测，有飞机观测，有气象卫星观测，有定时观测和非定常观测等等，怎样使用这些资料呢？一般来说，对气象资料的要求有两方面：一是可靠性，二是便于使人所周知。气象测站的分布是不规则的，因此我们只能得到这些不规则点上的气象资料，但是数值预报中的网格点是规则的，因而资料无法直接使用。另外，无论是用穿孔纸带或用电信号的形式将气象电报直接输入电子计算机，都要首先按照专门的程序进行译码、检查、整理。因为气象电报的内容，是按照气象电码格式编发的，而它的形式又是按照邮电电码格式编发的，
2、从观测、编码、发报，到传递、转换、接收等，在每个工序和环节上，都存在着出错的可能性。因此，我们所接收到的气象电码，不可避免地存在着一些错误或不妥之处。所以，要正确使用这些气象资料，必须经过必要的处理。比较简单的资料处理可分为以下几个方面：2.1.1记录错情判断 在利用接收到的气象资料之前，首先对资料要作错情判断。一般的做法是根据不同等压面上各种要素值的大小，给出相应略大于其最大值的一个数作为其上限值；也给出相应的略小于其最小值的一个数作为其下限值。例如在我国范围内的各测站，冬季500百帕层上的位势高度最小值不超过500位势什米，我们则取500作为其下限值；最大值不超过600位势什米，我们就取6
3、00作为其上限值。然后利用比较大小的子程序，由计算机对每一组数据进行判断，凡是大于上限值或小于下限值的记录，我们就认为它是错误的，予以舍掉，作为缺测记录。另外还可利用气象要素在时间变化上的连续性和空间分布上的连续性，来判断一个气象要素记录是否错误。对于错误记录，可用下面介绍的补缺测记录的方法，另外补一个值。2.1.2 补缺测或漏传记录 一般可把缺测或漏传记录的测点看作是一个网格点，然后由下节介绍的客观分析方法，利用周围已有的测站记录，插一个值补上。2.1.3 实测风矢量的分解 气象台站观测到的风场资料，是一个既有大小又有方向的风矢量，为便于该资料的利用经常将实测风分解为东西和南北两个分量。分别
4、用u、v表示，并规定：u向东为正，V向北为正。其数量值分别由下式计算： (2.1) 图2.1实测风矢量的分解 这里为实测风速值，为测风报告发布的风向度数。如，实测风为西南风，为240，风速10m/s，算得： 2.1.4 风场订正诊断分析一般都是在有限区域内进行的，多数都采用正方形网格。客观分析后所得到的网格点上的u、v分量值，并不处处与网格区的，轴平行，因此还必须进行风向订正。因为只有在基线上的网格点，其东西、南北方向与，方向一致，其他网格点上东西、南北方向与，方向总有一个偏差角，这显然会给计算带来误差，特别是当计算范围取得较大时，边缘的网格上，这种风向的误差显得更加突出。不进行适当订正是不行
5、的。如图2.2，先考虑在基线以西的某一网格点A，N是北极，NO为基线，NA和CA分别为经过A点的经线和纬线，MA和LA分别和这个正方形网格系统的X轴及Y轴平行。假设A点的风速在经纬的分量分别为u和v。而在网格的X,Y方向的分量分别为、，由于A点所在的经线不与基线相平行，故、分量和、分量彼此也不平行，而是有一夹角，ANO。由图可看出，它们之间有如下的换算关系： (2.2)图2.2风场订正示意图在基线以西的网格AO0，所以a=sin-1(AO/AN)0；在基线以东的网格，AO0，所以a0；在基线上的网格，0，则、。即不必订。 设在某网格点上，10米秒，45由 (2.2)式订正后 2.1.5 平滑和
6、滤波 气象观测资料，总存在着各种各样的误差。比如由气象仪器安装不标准等带来的，非偶然性误差（器差）和由工作人员在观测、编码、发收报等造成的偶然性误差，以及将要素值内插到网格点上时，产生的舍入和插值误差等等。无疑，这些误差都将会影响计算的结果，为了减少误差的影响，通常在计算之前先对原始资料等进行平滑和过滤（滤波），滤掉那些次要的小的天气意义的东西，而保留和突出主要的量，或者，为了研究的需要滤去资料中某些波长的量，而保留与问题有关的量。现分别介绍如下：(1) 一维平滑算子 这是最简单的平滑算子。利用同一直线上三点的资料，又称三点平滑算子 (2.3)这里表示第j点平滑后的值，表示第j点平滑前的值。S
7、为平滑系数（可正可负）。该平滑算子对j点对称，其权重除j+1，j和j-1点外均为0。图2.3 一维平滑算子对函数(x)可展成富氏级数，在点，可写成： (2.4)这里C为常数，A为波动的振幅。为波数，L为波长，为位相。同样在和点，函数(x)可写成 (2.5)将，代入(2.3)式得： (2.6)比较(2.4)、(2.6)式，平滑后的波相未变，改变的只是波的振幅，平滑后的振幅为： (2.7) 令R，称响应函数，即平滑波幅同原波幅之比，表示平滑后的变化。显然，R1，表示平滑后波幅一样，R1，表示平滑后使原波幅被放大。由(2.7)式知： 或者 (2.8)可见对于固定的网格距，响应函数R只与波数k（波长L
8、）以及平滑系数S有关。 由于0sin2(px/L) 1若希望平滑后使原波动衰减，以致滤掉（但不希望出现反位相情况），由只须要有0R1，于是由(2.6)式知：0S1/2。倘若希望平滑后使原波动增幅，即R1，则必有平滑系数S2x的波，可以使其波幅有所恢复。比如取L=6x时： 一次平滑得 R1(6x,1/2)=0.75 二次平滑得 R1-2=R(6x,1/2)R (6x,-1/2)=0.94 二者相比，二次平滑使该波幅恢复了19。表明这种平滑对保留长波是有益的。（2）二维平滑算子 对于平面的问题，须进行二维空间的平滑，把一维推广到二维有两种处理方法： 将计算的场先分别在X方向和Y方向进行平滑，然后取
9、平均，即 (2.12)其中: 这里用到的是，点及其前、后、左、右共五个点的资料，故称为五点平滑格式。 将场先在一个方向平滑，然后再在另一方向上平滑，即 (2.13)其中：2fij同(2.12)，2*fij=fi+1,j+1+fi+1,j-1+ fi-1,j+1+ fi-1,j-1-4fi,j。这里用的是，点及其前后，左右及前后点的左右点（式左右点的前后点）共九个点的资料，故称为九点平滑格式。用Kx，Ky，Lx，Ly分别表示X，Y方向上的波数和波长，其平面波的表示形式可写成 其中为振幅。以之代入(2.12),(2.13)易得其相应的响应函数 (2.14) (2.15)2.1.6 尺度分离 实际的
10、大气运动，包含了各种尺度的天气系统，为着研究的需要，经常要将实际的扰动，分离成不同尺度的波。应用富里叶级数将扰动展成不同尺度的波是常用的分离状的谱分析方法，但是它只能用于沿整个纬圈的半球性分析，不便考察某一特定地区和某些特定尺度系统等之间的关系。而利用上节所讲的平滑滤波技术同时可以满足这些方面的要求，同时又有着明显的天气意义。 为了分离不同尺度的波，须事先设计一个合适的滤波器。 对于任意天气变量A，可以写成： A=+(A-)。表示平滑后的A场，显然：(A- )表示原始场与平滑场的差。如果在作平滑的过程中，滤去高频波（短波），保留低频波（长波），这种平滑就称为一个低通滤波器，它的作用是让低频波通
11、过。表示场的低频部分，即大尺度运动。而(A- )部分，为该场的高颇（短波）部分，即小尺度系统。对(A- )进行运算，也就相当于一个高通滤波器。这里问题的关键，是如何选择滤波器，如何选择平滑系统和平滑次数，要看对具体问题的分析。比如要想分离出短波(波长L为500km左右）可以选用三点平滑算子，取S1/2，连续平滑三次，就得到单波响应函数：R1-3=(1-sin2(px/L)3。 如果计算网格距d取为100km，在场中波长为5d(500km)以下的波，衰减了70以上。保留了30以下，而对波长L=20d(2000km)以上的波，竟保留了95以上。这表明，经过如此过滤之后，在场中主要是波长为2000公
12、里以上的天气波。而在(A- )场中，主要是波长在500km以下的次天气尺度的短波系统。 如果研究的是面（二维）上的问题，即可采用二维平滑，如用九点平滑算子： 其中 (2.16)将实测风场，分解为u、v分量，取d100km，S=1/2，u、v分别用(2.15)式连续平滑三次，然后再合成，即得到相应的流场。 此外，还可以根据这种方法，设计只允许某一范围的波通过，而其它波不能通过的所谓带通滤波器，以分析研究某种波长的波的活动情况，这里变不再介绍了。2.2 客观分析诊断分析一般所需要的资料是网格点上的，而常规的气象观测资料是在固定地点(地面和高空观测站)和固定时间观测到的。为了由这些离散的分布不规则的
13、资料计算出某些物理量，从原则上必须得到每一观测变量在时、空上呈连续分布的场。但实际上并不需要观测值（如位势高度）在(X,Y)平面上连续变化的分布。因为一般可用有限差分方法来计算所需要的导数和梯度，这时只需要把空间上分布不均匀的台站资料内插到规则分布的网格点上就行了。为了得到网格上的资料，可采用两种方法进行内插：一种是主观内插法，即手工分析各种气象要素场的等值线，然后按网格点读取格点数，这种方法叫主观分析；另一种方法是根据直接联系格点值与台站值的方程，从数值上(用计算机)进行内插。这种方法叫客观分析。所用的方程或函数可以是不同的。常用的有限元、多项式、样条等，数值天气预报中还常使用逐步订正法、最
14、优插值法、谱方法、变分法等。另外，曲面拟合方法用于台站稀少的地区的物理量计算，也是比较好的方法。本节我们主要介绍有限元法、多项式法和逐次订正法。 客观分析方案可进行两种类型分析：一种是向量场分析，例如风场，其中所处理的资料不公有量值，还有方向；另一种处理量只有量值的资料，如温度场、湿度场等，这是标量场分析。2.2.1有限元素法图2.5三角形区域单元将我们要进行物理量计算的区域内的所有测站划分成有限个三个角形单元，如图2.5就是其中一个单元，A、B、C表示三个邻近的测站，O是内的一个网格点，若三个测站的要素值在一个平面上，根据有限元素法，它们可按其位置坐标展开以下的线性关系式： (2.17)其中
15、(XA,YA), (XB,YB)，(XC,YC)，分别为A、B、C三测站相对于任意坐标系的位置坐标。SA，SB，SC分别为测站所测得的要素值，这些都是已知的。因此，若有三组观测资料，根据线性代数中的克莱姆法则或高斯消去法等方法，就可求解出方程组的三个系数L1、L2、L3。 对于内的任意一个风格点O，已知它在同一坐标系里的位置坐标(X0,Y0)，则表示O点的要素值S0可由下式求得： S0=L1+L2X0+L3Y0 (2.18) 以上介绍的是有限元素法的一般原理，要用它来内插出所有网格点上的要素值。还要借助于电子计算机。具体做法是：先任取一网格点(i,j)并以此为圆心，以一个网格距为扫瞄半径画圆，
16、然后依次计算各测站到该网格点(i,j)的距离。凡是距离小于或等于扫瞄半径的测站被入选，大于的被舍去，这样一个过程。称为第一次“扫瞄”。这次扫瞄后，若圆内测站数小于3，则可加大扫瞄半径(取1.5个网格距或两个网格距），按上述的过程再扫瞄一次，若圆内测站仍小于3个，可继续加大扫瞄半径，并重复上述过程，一直到圆内测站数大于或等于3个为止。如果是等于3个，就可由这3个测站的资料，构成三元一次线性方程组。用高斯一亚当消去法或用其它的方法求出三个系数L1、L2、L3。，将这三个系数网格点(i,j)的位置坐标代入(2.17)式，便可内插出该网格点上的要素值。如果圆内的测站数大于3个，则再调用一个比较距离的子
17、程序，选取离网格点(i,j)最近的一个测站。利用其要素值，按上述方法，内插出该网格点上相应的要素值。 网格点(i,j)的要素值被内插出后，再进行下一个循环，同样用上述方法内插出下一个网格点的要素值，这样依次循环下去，直至内插出所有网格点上的要素值为止。 实际上，由于观测误差等原因，通常希望有比3个测站还要多的资料，这样求得的系数，可以减少观测随机误差的影响。为了解决这个问题，常采用最小二乘法来确定系数，由此得到（2.18)式： (2.19)式中是以扫瞄半径所画圆内的测站数，n4，解方程组(2.19)求得系数、和，然后再利用公式(2.17)求得圆内格点上的要素值。2.2.2 多项式法 多项式法的
18、原理是寻找一个由多项式表示的的曲面来逼近网格点周围区域各测站实测的气象要素值。如果这个曲面被找到，它即可代表该气象要素在这一网格点附近的空间分布状态，从而可求得网格点上的要素值，具体做法如下：假如某一等压面上的位势高度分布Z(x,y)可以用一个次多项式表示，其表达式为 (2.20)次幂项数p=0P=1P=2P=3P=mam 式中X，Y是该等压面上的两个位置坐标，p，q分别是它们的次数，ak是xp-qyq的系数。如果想用一个二次曲面去逼近高度场，那么可取m2；如果想用一个三次曲面去逼近高度场，那么可取m3；若想用一个更高次曲面去逼近高度场，那么m取更大的数值。现在我们所面临的问题是如何运用网格点
19、周围测站已知的一组实测高度值来确定系数，使得(2.20)式能最佳地逼近网格附近的高度场。对于给定的一组观测资料来说，坐标(xi,yi)(i=1,2,.n)是已知的，实测高度值Ziob也是已知的，m值是按我们的要求确定的，只有ak是未知数，即待确定的。 当(2.19)式中m2时，二次多项式的一般形式为 (2.21) 当m3时，三次多项式的一般形式为 (2.22) (2.21)式中含有六个待定系数，至少要有六组观测值才能求解。(2.2)式中含有10个待定系数，至少需要10组观测资料。m与待定系数个数M的关系为，。实际上，如有限元素一节中所述，为了减少观测随机误差的影响，得到最佳的逼近，常采用最小乘
20、二法来确定系数。如果所考虑的网格点附近有n个测站，n个测站实测高度值为，根据各测站的坐标，由多项式可以再得相应的高度值，某测站i的实测值与计算值的差为，称为残差，最小二乘法就是取残差的平方和为最小。残差平方和为： (2.23)上式为M元二次函数,自变量为(aj),残差平方和为最小的条件是它的一阶偏导数为零，即： (j=1,2,.M) (2.24) 令，()，则，代入(2.23)得： 将此式代入(2.24)得： (j=1,2,.M) 整理得： (j=1,2,.M) (2.25) 令 (2.26)代如(2.24)得： j=1,2,.M (2.27) 上式为求解M个未知数的线性方程组(若m2，则M6
21、；若m3，M10)，形式与(2.19)式相同，解这个方程组，即可确定系数，从而最后确定我们所要寻求的多项式，并用此内插出网格点上的要素值。在电子计算机上实现的方法和步骤，与有限元素法基本相同。2.2.3 逐步订正法 逐步订正法就是将格点周围站记录与终点值进行比较，用格点周围不同半径范围内各测站的观测数值情形与估计值之差的加权平均作修正量逐步对其订正，最终使格点分析值与周围测站记录相比达到完全合理为止。 下面以位势高度的格点分析值的计算为例，简要介绍一下逐步订正法。 以格点为中心，以R为半径划一个圆，设落在圆内的观测值有N个(图2.6)，我们将根据这N个观测值来订正估值。 将初估场(预备场)的格
22、点值 内插到某一个观测点上，例如点，设内插出的高度值为，若点的高度观测值为，则观测值与内插值的差可以算出：对圆内所有个观测值都算出它们与各自内插值的差，.，并计算这些差值的加权平均 (2.28)就是来自高度场的订正值。式中为权重系数。有很多种形式，但最主要的特征是和格点到点的距离的平方成反比，其物理意义是，观测值对格点值的影响程度随着它们之间距离的平方而减小。显然，这种公重系数是各向同性的，它在平面上的分布是一个个同心圆。例如可以取，式中是格点到点的距离，R是扫描半径。也可以取为，是经验系数，随不同层次，不同要素而异。 如果在求取格点的高度分析值时还要使用风场的观测资料，则还需要求出来自风场的
23、订正值。图2.6 圆形区域 设所要求取的格点分析值为 ，则有 假定在图2.5所示的圆内位势高度是线性变化的，利用台劳展式将在点展开： 即： (2.29)将地转风关系 代如得： 对圆内所有有测风报告的站点加权平均，便得到来自风场的订正值 (2.30)M为圆内有测风报告的站点个数，显然MN于是，总订正值为： (2.31)和为加权系数，可根据需要加以调整。比如，对于风场好用的地区，可以取得大一些，反之亦然，但和之间必须满足以下关系： 由此得到格点的分析值为 (2.32) 这是对于影响半径R进行的第一次许订正，然后逐渐缩小R进行重复订正。每次订正后的分析场用作下一次订正的估计值。在连续订正过程中，缩小
24、影响半径可以去掉估计场的大尺度误差，使分析场越来越逼近观测结果。 在逐次订正中还可以采用不同形式的权重系数。2.2.4拉格朗日垂直插值法气象资料在垂直方向上往往是不等距的，为计算方便起见，常常需要将资料（观测结果）进行垂直插值。图2.7线性插值多项式 最简单的插值是线性插值。设对应变量X有函数值Y，对点所对应的函数值为，取Y=P(X)线性插值，就是在点上有函数值，两点作直线于是可得线性插值。 (2.33)也就是一次插值多项式。显然，由此所示得的Y值其误差较大(如图2.7)。常用的是二次插值多项式，为： (2.34)虽然简单，但不够直观。 考虑用任意点的函数本身乘上一个适当的系数，来表示插值多项
25、式，就是拉格朗日插值。 (2.35)式中的系数，又是一个阶多项式。 由；可知应满足条件： (2.36) 考虑，除j=0，=1外，其它点=0，故可将表示为多项式： (2.37)这里C是一具待定的常数，将时，=1代入(2.35)得： 于是： (2.38)同理可求得 的表达式 (2.39)最后拉格朗日插值公式的普遍形式为： (2.40)拉格朗日的二次(n2)插值公式为： (2.41)常被用来求各等压面上的要素值或其他物理量。 由于气压随高度呈对数变化，应用插值公式时常取气压的对数形式lnP。例如，已知探空资料500hPa、700hPa、850hPa的高度值分别为5500、3000、1500（单位：位
26、势米），用拉格朗日二次插值公式求600hPa的高度 2.2.5 观测资料的四维同化 随着观测手段的不断改进和现代化，不仅资料的精度有了提高，资料量扩大。而且新资料不断出现，如雷达探测，卫星云图资料等。充分、合理地使用这些新资料，也是提高天气预报准确率的一个重要方面。 我们知道，原来各个气象台站的观测，都是在统一规定时间内进行的。所以得到的都是定时观测资料。其中最重要的是卫星资料。卫星观测资料是从人造地球卫星上发下来的。人造卫星有两种，一种是极地轨道卫星，它通过极地的上空，绕地球旋转一周，约需90分钟，即隔90分钟就可以获得某地区附近的一次观测资料。另一种是地球同步卫星，它围绕地球旋转的角速度和
27、地球自转的角速度是相等的。它可以观测到的地面的范围，是半径为4000公里多的圆形区域。所以如果在地球赤道上空，均匀地分布着45个地球同步卫星，那么，40N40S这样一个广大的中低纬度地区，都可以被观测到。而且每隔20分钟或更短的时间，就可以得到一次观测资料。 卫星观测资料主要有两种，一种是云图,另一种是辐射资料，可以推算出各个高度上的气温。根据云图，可以反推风和气压的分布。特别是根据云在很短时间（如几分钟）内的变化，可以得出云所在高度上的风。因为在这很短时间内，完全可以认为云的变化是由于风的吹动引起的。此外，还有飞机观测报告，主要是超音速飞机，它的飞行高度一般接近10公里或更高，所以飞机报告的
28、资料多是在250hpa等压面上。在进行250hpa等面分析时，除无线电测风资料外，还可以利用大量的飞机报告和卫星资料。根据卫星资料中云的移动可以推算等压面层上的风向风速。根据云的温度还可推算气压的大小，而云的温度是由红外辐射测量的。特别在低纬度和海洋地区更需要充分利用卫星资料（多数是处于250hpa附近）和飞机报告，才能得到可靠的分析结果。 这样，最近一、二十年来，由于大量的非定时观测资料的出现，客观分析就不只是空间三维的了，还必须把时间这一维也加进去，即把各不同时刻的观测资料纳入统一的分析预报中来，使之自然满足一定的协调条件，此种分析方法称为“四维同化”。 四维同化实际上是一种基于热力-流体力学方程,以数值天气预报(NWP)模式为基础的内插方法.这种方法通常用NWP模式预报值逐日对资料进行内插.在模式中,除了无线电探空资料外,还把诸如飞机探测报告,船测报告,卫星资料等非常规资料包括在一个连续的四维资料同化系统中.如果同化正确,这些附加的资料应该明显改进逐日天气图.通过对基本变量(两个风分量,温度,位势高度和湿度) 逐日客观分析进行平均,就可以得到格点上的环流统计特征量.NWP客观分析的格点资料均匀地覆盖全球,实用起来较方便.22 / 17
