天气诊断分析报告.docx
- 文档编号:28603832
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:24
- 大小:221.22KB
天气诊断分析报告.docx
《天气诊断分析报告.docx》由会员分享,可在线阅读,更多相关《天气诊断分析报告.docx(24页珍藏版)》请在冰豆网上搜索。
天气诊断分析报告
第二章资料处理和客观分析
§2.1资料处理
随着气象观测手段的发展和现代化,得到的资料数量增多,门类拓广,有常规站点观测,有非定点海洋观测,有飞机观测,有气象卫星观测,有定时观测和非定常观测等等,怎样使用这些资料呢?
一般来说,对气象资料的要求有两方面:
一是可靠性,二是便于使人所周知。
气象测站的分布是不规则的,因此我们只能得到这些不规则点上的气象资料,但是数值预报中的网格点是规则的,因而资料无法直接使用。
另外,无论是用穿孔纸带或用电信号的形式将气象电报直接输入电子计算机,都要首先按照专门的程序
进行译码、检查、整理。
因为气象电报的内容,是按照气象电码格式编发的,而它的形式又是按照邮电电码格式编发的,从观测、编码、发报,到传递、转换、接收等,在每个工序和环节上,都存在着出错的可能性。
因此,我们所接收到的气象电码,不可避免地存在着一些错误或不妥之处。
所以,要正确使用这些气象资料,必须经过必要的处理。
比较简单的资料处理可分为以下几个方面:
2.1.1记录错情判断
在利用接收到的气象资料之前,首先对资料要作错情判断。
一般的做法是根据不同等压面上各种要素值的大小,给出相应略大于其最大值的一个数作为其上限值;也给出相应的略小于其最小值的一个数作为其下限值。
例如在我国范围内的各测站,冬季500
百帕层上的位势高度最小值不超过500位势什米,我们则取500作为其下限值;最大值不超过600位势什米,我们就取600作为其上限值。
然后利用比较大小的子程序,由计算机对每一组数据进行判断,凡是大于上限值或小于下限值的记录,我们就认为它是错误的,予以舍掉,作为缺测记录。
另外还可利用气象要素在时间变化上的连续性和空间分布上的连续性,来判断一个气象要素记录是否错误。
对于错误记录,可用下面介绍的补缺测记录的方法,另外补一个值。
2.1.2补缺测或漏传记录
一般可把缺测或漏传记录的测点看作是一个网格点,然后由下节介绍的客观分析方法,禾U用周围已有的测站记录,插一个值补上。
2.1.3实测风矢量的分解
气象台站观测到的风场资料,是一个既有大小又有方向的风矢量,为便于该资料的
利用经常将实测风分解为东西和南北两个分量。
分别用u、v表示,并规定:
u向东为正,
V向北为正。
其数量值分别由下式计算:
(2.1)
uVsin(/180)
vVcos(/180)
这里V为实测风速值,为测风报告发布的风向度数。
如,实测风为西南风,为
240°,风速10m/s,算得:
u10sin(240/180)8.7(m/s)
v10cos(240/180)5(m/s)
2.1.4风场订正
诊断分析一般都是在有限区域内进行的,多数都采用正方形网格。
客观分析后所得到的网格点上的u、v分量值,并不处处与网格区的x,y轴平行,因此还必须进行风向订正。
因为只有在基线上的网格点,其东西、南北方向与x,y方向一致,其他网格点上东西、南北方向与x,y方向总有一个偏差角,这显然会给计算带来误差,特别是当计算范围取得较大时,边缘的网格上,这种风向的误差显得更加突出。
不进行适当订正是不行的。
如图2.2,先考虑在基线以西的某一网格点A,N是北极,NO为基线,NA和CA分别为经过A点的经线和纬线,MA和LA分别和这个正方形网格系统的X轴及丫轴平行。
假设A点的风速在经纬的分量分别为u和V。
而在网格的X,Y方向的分量分别为u'、v',由于A点所在的经线不与基线相平行,故u、v分量和u'、v'分量彼此也不平行,而是有一夹角a,a=ZANO由图可看出,它们之间有如下的换算关系:
(2.2)
u'ucosvcos(90)ucosvsin
v'vcosucos(90)vcosusin
N
图2.2风场订正示意图
C
在基线以西的网格AO>0所以=sin-1(AO/AN)>0;在基线以东的网格,AO<0所以<0;在基线上的网格,a=0,则u'=u、v'=v。
即不必订。
设在某网格点上,u=v=10米/秒,a=45°由(2.2)式订正后
u10cos4510sin4514(米/秒)
v10cos4510sin450(米/秒)
2.1.5平滑和滤波
气象观测资料,总存在着各种各样的误差。
比如由气象仪器安装不标准等带来的,
非偶然性误差(器差)和由工作人员在观测、编码、发收报等造成的偶然性误差,以及将要素值内插到网格点上时,产生的舍入和插值误差等等。
无疑,这些误差都将会影响计算的结果,为了减少误差的影响,通常在计算之前先对原始资料等进行平滑和过滤(滤
波),滤掉那些次要的小的天气意义的东西,而保留和突出主要的量,或者,为了研究的需要滤去资料中某些波长的量,而保留与问题有关的量。
现分别介绍如下:
(1)一维平滑算子
这是最简单的平滑算子。
利用同一直线上三点的资料,又称三点平滑算子
SS
fjfj2fj1fj12fj(1S)fj2fj1fj1(2.3)
这里fj表示第j点平滑后的值,fj表示第j点平滑前的值。
S为平滑系数(可正可负)。
该平滑算子对j点对称,其权重除j+1,j和j-1点外均为0。
对函数f(X)可展成富氏级数,在fj点,可写成:
ik(x.)
fj(x)CAej
和Xj!
点,函数f(X)可写成
比较(2.4)、(2.6)式,平滑后的波相未变,改变的只是波的振幅,平滑后的振幅为:
令R=A/A,称响应函数,即平滑波幅同原波幅之比,表示平滑后的变化。
显然,R=1,表示平滑后波幅一样,R<1,表示平滑使原波幅衰减。
R-0,表示平滑后使原波动被消失(即波动全被滤掉)。
R>1,表示平滑后使原波幅被放大。
由(2.7)式知:
R(k,s)A/A1S(1coskx)
或者R(L,S)12S?
sin2(x/L)(2.8)
可见对于固定的网格距,响应函数R只与波数k(波长L)以及平滑系数S有关。
由于0Wsin2(△x/L)<1若希望平滑后使原波动衰减,以致滤掉(但不希望出现反位相情况),由只须要有0WR<1,于是由(2.6)式知:
0WSW1/2。
倘若希望平滑后使原波动增幅,即R>1,则必有平滑系数S<0b
作为一特例,取S=1/2,即最大平滑系数,此时响应函数为
R(L,1/2)1sin2(x/L)cos2(x/L)(2.9)
对于L=2Ax的波,R=0,表明波长为2倍网格距的波,通过这种平滑,可认完全被滤掉。
对于L2x的波,平滑使波幅有不同程度的衰减,但由于余弦函数,在0/2之
间是减函数,随角度增加,余弦函数减小,L越大,x/L越小,R越大,表示平滑后波幅随波长的增大而减衰得越来越小。
取L=3Ax时,由(2.7)式知,R(3Ax,1/2)=0.25,原波幅衰减了75%。
取L=6Ax时,R(6Ax,1/2)=0.75原波幅衰减了25%。
取L=10Ax时,R(10Ax,1/2)=0.905,原波幅衰减得更少,不足10%。
可见,取S=1/2滤波时,虽然可以滤去高频波,但同时也削弱了天气波,不甚理想。
理想滤波应该是保留需要的波,滤去所不需要的波,从响应函数曲线(图2.4)上看图
形最好近似为矩形。
图2.4S=1/2的响应函数曲线
倘若为去掉短波,并且尽可能少的改变长波,可以采用不同平滑系数,仍用同一平
滑算子,函数进行多次平滑的办法。
可以证明,
取平滑系数S,S,……S
n,作n次平滑后的响应函数R为:
R1n(k,S)
n
F^R2?
?
?
Rn1Si(1
i1
coskx)或者
R1n(k,L)
n
2
12Sisin(kx)
i1
(2.10)
作为一个例子,这里举一个二次(n=2)平滑的情况,并且将平滑系数分别取为
S1=1/2,S2=1/2,平滑后的响应函数,由(2.7)知
Ri2R(L,1/2)R2(L,1/2)1sin2(x/L)?
1sin2(x/L)
1sin4(x/L)(2.11)
同(2.7)式,即取S=1/2的一次平滑情况比较。
取上述二次症滑后,对波长L>2Ax
的波,可以使其波幅有所恢复。
比如取L=6Ax时:
一次平滑得R1(6△x,1/2)=0.75
二次平滑得R1-2=R(6Ax,1/2)?
R(6△x,-1/2)=0.94
二者相比,二次平滑使该波幅恢复了19%。
表明这种平滑对保留长波是有益的。
(2)二维平滑算子
对于平面的问题,须进行二维空间的平滑,把一维推广到二维有两种处理方法:
①将计算的场先分别在X方向和丫方向进行平滑,然后取平均,即
玷
1fi
fi,j
fj
i,j
S2
f2f..
i,ji,j
4
(2.12)
2
其中:
2fi.i,j
f“
i1,j
fi1,jfi,j1
fi,j14fi,j
这里用到的是i,j点及其前、后、左、右共五个点的资料,故称为五点平滑格式。
②将场先在一个方向平滑,然后再在另一方向上平滑,即
2
jss
fi,jjfi,jfi,jS(1S)2fi,j宁2fi,j(2.13)
24
其中:
▽2fij同(2.12),▽*fij=fi+1,j+1+fi+1,j-1+fi-1,j+1+fi-1,j-1-4fi,j。
这里用的是i,j点及其前后,左右及前后点的左右点(式左右点的前后点)共九个点的资料,故称为九点平滑格式。
用Kx,Ky,Lx,Ly分别表示X,丫方向上的波数和波长,其平面波的表示形式可写成
CAei(kxXkY丫)或CAe(2X/Lx2丫/Ly)
其中A为振幅。
以之代入(2.12),(2.佝易得其相应的响应函数
1
R五占一RXRy1Ssin2(kXX/2)sin2(kYY/2)
2
1Ssin2(X/LX)sin2(Y/LY)(2.14)
R九点RX?
Ry12Ssin2(kXX/2)?
12Ssin2(kYY/2)
12Ssirf(X/Lx)?
12Ssirf(Y/Ly)(2.15)
2.1.6尺度分离
实际的大气运动,包含了各种尺度的天气系统,为着研究的需要,经常要将实际的扰动,分离成不同尺度的波。
应用富里叶级数将扰动展成不同尺度的波是常用的分离状的谱分析方法,但是它只能用于沿整个纬圈的半球性分析,不便考察某一特定地区和某些特定尺度系统等之间的关系。
而利用上节所讲的平滑滤波技术同时可以满足这些方面的要求,同时又有着明显的天气意义。
为了分离不同尺度的波,须事先设计一个合适的滤波器。
对于任意天气变量A,可以写成:
A=?
+(A-?
)。
?
表示平滑后的A场,显然:
(A-?
)表示原始场与平滑场的差。
如果在作平滑的过程中,滤去高频波(短波),保留低频波
(长波),这种平滑就称为一个低通滤波器,它的作用是让低频波通过。
表示场的低频部分,即大尺度运动。
而(A-?
)部分,为该场的高颇(短波)部分,即小尺度系统。
对(A-?
)进行运算,也就相当于一个高通滤波器。
这里问题的关键,是如何选择滤波器,如何选择平滑系统和平滑次数,要看对具体问题的分析。
比如要想分离出短波(波长L为500km左右)可以选用三点平滑算子,取S=1/2,连续平滑三次,就得到单波响应函数:
R-3=(1-sin2(△x/L))3。
如果计算网格距d取为100km在?
场中波长为5d(500km)以下的波,衰减了70%以上。
保留了30%以下,而对波长L=20d(2000km)以上的波,竟保留了95%以上。
这表明,经过如此过滤之后,在?
场中主要是波长为2000公里以上的天气波。
而在(A-?
)场中,主要是波长在500km以下的次天气尺度的短波系统。
如果研究的是面(二维)上的问题,即可采用二维平滑,如用九点平滑算子:
将实测风场,分解为u、v分量,取d=100kmS=1/2,u、v分别用(2.15)式连续平滑
三次,然后再合成,即得到相应的流场。
此外,还可以根据这种方法,设计只允许某一范围的波通过,而其它波不能通过的
所谓带通滤波器,以分析研究某种波长的波的活动情况,这里变不再介绍了。
§2.2客观分析
诊断分析一般所需要的资料是网格点上的,而常规的气象观测资料是在固定地点
(地面和高空观测站)和固定时间观测到的。
为了由这些离散的分布不规则的资料计算出某些物理量,从原则上必须得到每一观测变量在时、空上呈连续分布的场。
但实际上并不需要观测值(如位势高度)在(X,Y)平面上连续变化的分布。
因为一般可用有限差分方法来计算所需要的导数和梯度,这时只需要把空间上分布不均匀的台站资料内插到规则分布的网格点上就行了。
为了得到网格上的资料,可采用两种方法进行内插:
一种是主观内插法,即手工分析各种气象要素场的等值线,然后按网格点读取格点数,这种方法叫主观分析;另一种方法是根据直接联系格点值与台站值的方程,从数值上(用计算
机)进行内插。
这种方法叫客观分析。
所用的方程或函数可以是不同的。
常用的有限元、多项式、样条等,数值天气预报中还常使用逐步订正法、最优插值法、谱方法、变分法等。
另外,曲面拟合方法用于台站稀少的地区的物理量计算,也是比较好的方法。
本节我们主要介绍有限元法、多项式法和逐次订正法。
客观分析方案可进行两种类型分析:
一种是向量场分析,例如风场,其中所处理的资料不公有量值,还有方向;另一种处理量只有量值的资料,如温度场、湿度场等,这是标量场分析。
2.2.1有限元素法
将我们要进行物理量计算的区域内的所有测站划分成有限个三个角形单元,如图
2.5就是其中一个单元,AB、C表示三个邻近的测站,0是ABC内的一个网格点,若
三个测站的要素值在一个平面上,根据有限元素法,它们可按其位置坐标展开以下的线
性关系式:
其中(XA,YA),(XB,YB),(XC,YC),分别为A、BC三测站相对于任意坐标系的位
置坐标。
SASBSC分别为测站所测得的要素值,这些都是已知的。
因此,若有三组观测资料,根据线性代数中的克莱姆法则或高斯消去法等方法,就可求解出方程组的三个系数L1、L2、L3。
对于ABC内的任意一个风格点0,已知它在同一坐标系里的位置坐标(X°,Yo),则表示0点的要素值So可由下式求得:
So=Li+L2Xo+LsY)(2.18)
以上介绍的是有限元素法的一般原理,要用它来内插出所有网格点上的要素值。
还要借助于电子计算机。
具体做法是:
先任取一网格点(i,j)并以此为圆心,以一个网格
距为扫瞄半径画圆,然后依次计算各测站到该网格点(i,j)的距离。
凡是距离小于或等
于扫瞄半径的测站被入选,大于的被舍去,这样一个过程。
称为第一次“扫瞄”。
这次
扫瞄后,若圆内测站数小于3,则可加大扫瞄半径(取1.5个网格距或两个网格距),按上述的过程再扫瞄一次,若圆内测站仍小于3个,可继续加大扫瞄半径,并重复上述过
程,一直到圆内测站数大于或等于3个为止。
如果是等于3个,就可由这3个测站的资料,构成三元一次线性方程组。
用高斯一亚当消去法或用其它的方法求出三个系数Li、
L2、Lao,将这三个系数网格点(i,j)的位置坐标代入(2.17)式,便可内插出该网格点上的要素值。
如果圆内的测站数大于3个,则再调用一个比较距离的子程序,选取离网格点(i,j)最近的一个测站。
利用其要素值,按上述方法,内插出该网格点上相应的要素值。
网格点(i,j)的要素值被内插出后,再进行下一个循环,同样用上述方法内插出下一个网格点的要素值,这样依次循环下去,直至内插出所有网格点上的要素值为止。
实际上,由于观测误差等原因,通常希望有比3个测站还要多的资料,这样求得的
系数,可以减少观测随机误差的影响。
为了解决这个问题,常采用最小二乘法来确定系数,由此得到(2.18)式:
nnn
L2XiL3YS
i1i1i1
nnn
L2Xi2L3XiYSXi
i1i1i1
nnn
L2XiYL3Y2SY
i1i1i1
n>4,解方程组(2.19)求得系数L1、L2和L3,
然后再利用公式(2.17)求得圆内格点上的要素值。
2.2.2多项式法
多项式法的原理是寻找一个由多项式表示的的曲面来逼近网格点周围区域各测站
实测的气象要素值。
如果这个曲面被找到,它即可代表该气象要素在这一网格点附近的空间分布状态,从而可求得网格点上的要素值,具体做法如下:
假如某一等压面上的位势高度分布Z(x,y)可以用一个m次多项式表示,其表达式为
mpi
Z(x,y)akxpqyqk-p(p1)q1(2.20)
p0q02
次幂
项数
p=0
①
P=1
②③
P=2
④⑤⑥
P=3「
⑦⑧⑨⑩
?
P=m
?
?
?
?
?
?
?
?
?
?
?
?
am
式中X,丫是该等压面上的两个位置坐标,p,q分别是它们的次数,ak是xp-qyq的系数。
如果想用一个二次曲面去逼近高度场,那么可取m=2;如果想用一个三次曲面去逼
近高度场,那么可取m=3;若想用一个更高次曲面去逼近高度场,那么m取更大的数值。
现在我们所面临的问题是如何运用网格点周围测站已知的一组实测高度值来确定
系数ak,使得(2.20)式能最佳地逼近网格附近的高度场。
对于给定的一组观测资料来说,坐标(Xi,yJ(i=1,2,..…n)是已知的,实测高度值
Zob也是已知的,m值是按我们的要求确定的,只有ak是未知数,即待确定的。
当(2.19)式中2时,二次多项式的一般形式为
2
(2.21)
2
Z(x,y)aia?
xa3yaqXa§xyaey
当m=3时,三次多项式的一般形式为
(2.22)
Z(x,y)aia?
xa3yaqX2a§xyaey2
3223
a?
xagXya9xya^y
(2.21)式中含有六个待定系数,至少要有六组观测值才能求解。
(2.2)式中含有10
个待定系数,至少需要10组观测资料。
m与待定系数个数M的关系为,
1
M^(m1)?
(m2)。
实际上,如有限元素一节中所述,为了减少观测随机误差的影响,得到最佳的逼近,常采用最小乘二法来确定系数ak。
如果所考虑的网格点附近有n
个测站,n个测站实测高度值为Zi0b(i1,2,……n),根据各测站的坐标(Xi,yj,由多项
式可以再得相应的高度值,某测站i的实测值与计算值的差为Z°bZi,称为残差,最
小二乘法就是取残差的平方和为最小。
残差平方和为:
n
(2.23)
E(aj)(Zi0bZi)2
i1
上式为M元二次函数,自变量为(aj),残差平方和为最小的条件是它的一阶偏导数为零,
即:
E⑻)0
aj
(j=1,2,....
...M)
(2.24)
令Z0bBm1i,
pqq
XiyiBki
,(
1
k—P(P1)q
2
1
mp
M
p
乙akxi
yiak
Bki代入(2.23)
>得:
p0q0
k1
n
nM
2
E(ai)Z
3Zi2
■Ii
BM1i
akBki
将此式代入(2.24)
得:
i1
i1k1
n
M
2B…
j,i
BM1iakBki0
(j=1,2,...
...M)
i1
k1
整理得:
Mn
n
akBj
?
Bl-
1iki
Bji?
BM1i
(j=1,2,
......M)(2.25)
k1i1
i1
n
令Cjk
B..B,.jiki
j
1,2,……
M
’.(2.26)
i1
1
k
1,2,……
.M1
代如(2.24)得:
M
Cjkak
k1
CjM1
j=1,2,……M
(2.27)
上式为求解
M个未知数
ai的线性方程组(若m=2,则M=6;若m=3,
M=10),形式
),贝U,
与(2.19)式相同,解这个方程组,即可确定系数ai,从而最后确定我们所要寻求的多项
式,并用此内插出网格点上的要素值。
在电子计算机上实现的方法和步骤,与有限元素法基本相同。
2.2.3逐步订正法
逐步订正法就是将格点周围站记录与终点值进行比较,用格点周围不同半径范围内
各测站的观测数值情形与估计值之差的加权平均作修正量逐步对其订正,最终使格点分
析值与周围测站记录相比达到完全合理为止。
下面以位势高度的格点分析值的计算为例,简要介绍一下逐步订正法。
以格点为中心,以R为半径划一个圆,设落在圆内的观测值有N个(图2.6),我们将根据这N个观测值来订正估值。
将初估场(预备场)的格点值d内插到某一个观测点上,例如i点,设内插出的高度值为i,若i点的高度观测值为i,则观测值与内插值的差可以算出:
iii
对圆内所有N个观测值都算出它们与各自内插值的差1,2,……N,并计
算这些差值的加权平均
NN
CiWi(2.28)
i1i1
C就是来自高度场的订正值。
式中叫为权重系数。
Wi有很多种形式,但最主要的特征是和格点到i点的距离的平方成反比,其物理意义是,观测值对格点值的影响程度随着它们之间距离的平方而减小。
显然,这种公重系数是各向同性的,它在平面上的分布是一个个同心圆。
例如可以取W(R2ri2)(R2ri2),式中n是格点到i点的距离,R是扫描半径。
也可以取为W1(1『),是经验系数,随不同层次,不同要素而
异。
如果在求取格点的高度分析值时还要使用风场的观测资料,则还需要求出来自风场
的订正值CV0。
图2.6圆形区域
设所要求取的格点分析值为d,则有
假定在图2.5所示的圆内位势高度是线性变化的,利用台劳展式将d在1点展开:
x
对圆内所有有测风报告的站点加权平均,便得到来自风场的订正值
M
Wif(vxuy)i
i1
Md
W
i1
M为圆内有测风报告的站点个数,显然M〈N于是,总订正值为:
CaCavCv(2.31)
a和aV为加权系数,可根据需要加以调整。
比如,对于风场好用的地区,aV可以取得
大一些,反之亦然,但a和aV之间必须满足以下关系:
aaV1
由此得到格点的分析值为
(2.32)
这是对于影响半径R进行的第一次许订正,然后逐渐缩小R进行重复订正。
每次订正后的分析场用作下一次订正的估计值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 天气 诊断 分析 报告