引言

统计预报中的Kalman滤波方法是继传统MOS、PP方法后，被气象业务运用得较好的数值预报产品释用新方法^[1]。陆如华等(1994)^[2]用了Kalman滤波做连续量(地面气温、风场等)为对象的要素预报优于回归方法，且用Kalman滤波建立统计预报模型，既能适应自然天气季节变化，又能适应数值预报模式的更新改进。梁钰等(2006)^[3]用卡尔曼滤波做沙尘天气预报，但又采用“判别指标”；赵玉广等(2004)^[4]仍用PP方法做雾的预报；胡春梅等(2006)^[5]用逐步回归方法，采用气候持续因子，做西北太平洋热带气旋强度预报；林健玲等(2006)^[6]用条件数计算选取回归因子，做区域日降水量预报，且好于逐步回归方法；曾黎明等(2005)^[7]用逐步相似法对热带气旋进行路径相似预报；张建海等(2005)^[8]用多时刻因子代替单时刻因子，以改善统计预报方程质量；高洁等(2005)^[9]用事件概率回归和选用相关系数较高的因子，做(小概率)辐射雾雾消预报。另外，近年多有“支持向量机回归方法”在气象预报中的应用^[10]。

本文研究表明，Kalman滤波与回归一样可以归结成为MOS预报方法，从而可以用于对降水做MOS预报；通过设计并选用对于预报对象是线性化的预报因子，和采用多时刻因子，用相似-Kalman滤波做降水定性(概率)、定量预报；和采用单时刻的线性化预报因子，用相似-Kalman滤波做地面气温和风场预报。因是统计预报，这里又将历史样本的预报准确率当作预报对象的发生概率，从而完成“概率天气预报”。最终建立“湖北分县降水(概率)、地面气温和风场实时统计MOS预报系统”(以下简称“系统”)。

1 资料、方法与预报产品

“系统”用我国业务数值模式T213L31预报产品，取1000、925、850、700、600、500、400、300、250、200hPa共10层24、30、36、42、48小时预报模式大气资料。可将设计计算的网格点预报因子场，做球面坐标系上双三次曲面拟合^[11]，并按各个预报站点的经纬度进行插值。从而取得对应站点降水、地面气温和风场观测资料，并按预报时次建立各个站点的时序样本。

“系统”用Kalman滤波，且用了“条件”相似Kalman滤波，以做定性定量降水、地面气温和风场统计MOS预报。“条件”相似-Kalman滤波是指：当预报有较大天气变化的“条件”时，就做相似-Kalman滤波，这时要在历史样本中找到与当前预报最相似样本，否则，仍做单纯Kalman滤波，即无“条件”出现时，默认最相似样本就是当天时序最后样本。另外，降水概率预报也是使用统计方法所做客观预报，以区别主观预报与经验概率。

在物理学上，气温要素为一个自由度，而风和降水都是两个自由度，且降水还是不连续量(场)。故地面风场是分别做其两个水平分量预报。又降水分为定性与定量预报：定性预报是将样本中有关预报因子与预报对象都“+1、-1”化；而定量预报是在选样本时，为保证降水过程的时空连续性，限定了样本的降水(过程)发生站点数与平均雨量两个方面条件。

“系统”提供湖北全省(部分临省)定点(分县)、定时(6小时分辨率)、定量降水与概率、地面气温和风场MOS预报产品。其中风场预报包括两个内容：①风向和②风速(率)。而降水与概率预报包括三个内容：①有无降水发生；②若有降水发生则给出预报降水量，否则不给出预报降水量；③有无降水都给出预报降水发生概率。这里是用历史样本预报雨/晴准确率，代替预报有/无降水发生概率，并将概率分为九个等级：0~15~25~35~45~55~65~75~85~100。

2 预报因子选取 2.1 降水定性、定量预报因子

降水定量预报方程选取T213L31各个预报时次模式大气中4个预报因子：

R₁：水汽通量散度降水(mm/h)^[12-13]；

R₂：凝结函数降水(mm/h)^[12-13]；

R₃：气块不稳定降水^[12-13](R_A)和对流不稳定降水^[12-13](R_B)之线性组合(mm/h)；

R₄：T213L31产品降水预报场(mm/h)；

R₅：MOS降水预报场(mm/h)。

即有：

$ {R_1} = \int_{Ps}^{PT} {\mathop - \limits^\nabla } V{q_s}\frac{{{\rm{d}}p}}{g}\;\;\;\left( {RH \ge 80\% } \right) $

(1)

式(1)中V、RH分别为模式大气中的全风速和相对湿度，q_s为RH≥80%时的饱和比湿，p_s是地面气压，p_T为模式大气顶层，并取p_T=200hPa，g为地球重力常数。这里当R₁>0时，有水汽通量辐合与降水；而当R₁＜0时，有水汽通量辐散与“反降水”。显然，当整个气柱RH均小于80%时，R₁=0。

$ {R_2} = \int_{Ps}^{PT} - F\omega \frac{{{\rm{d}}p}}{g}\;\;\;\left( {RH \ge 80\% } \right) $

(2)

式(2)中F(≥0)为凝结函数^[12]，w为模式大气中的P坐标垂直速度^[12]。因ω＜0时R₂>0，有降水发生；而ω>0时R₂＜0，也有“反降水”物理意义，故R₂与R₁一样有正负分布，同时也有当整个气柱RH均小于80%时，R₂=0。

大气中还存在条件不稳定气团，所谓的“气块干稳定湿不稳定”^[12]：即气块湿绝热上升时的不稳定能量。积分因大气垂直不稳定能量释放而造成的气块不稳定降水：

$ {R_a} = \int_{Ps}^{{P_T}} {\left( {{q_{so}} - {q_s}} \right)} \frac{{{\rm{d}}p}}{g}\; $

(3)

式(3)中q_s为湿不稳定气块在绝热上升时变态饱和比湿，q_s0为其初值。

大气中还有因整(湿)层空气被抬升而出现对流不稳定能量^[12]。考虑到大气中水汽主要存在于对流层的中、低层，这里取p_s~900~800~700~600~500hPa共五层。当不稳定层被强迫抬升时，按定义是因其湿度层结的对流不稳定能量释放而造成了对流不稳定降水。则积分各层的对流不稳定降水：

$ {R_b} = \sum\limits_{i = 1}^5 {\int_{p_i^2}^{p_i^2} \Delta } {q_{si}}\frac{{{\rm{d}}p}}{g} $

(4)

式(4)中Δq_si=q_si²-q_si¹，q_si¹、q_si²分别是对应于p_i¹、p_i²各对流不稳定层经调整为对流稳定层后的下、上界面变态饱和比湿，R_b是前述五层之和。

因(3)、(4)式的气块不稳定能量和对流不稳定能量释放在理论上都是“瞬间”完成的，而实际应有一个重建不稳定湿层时间过程^[14]。可设，

$ {T_h} = \frac{{\int_{ps}^{{p_T}} {\frac{V}{S}qp\frac{{{\rm{d}}p}}{g}} }}{{\int_{{p_s}}^{{p_T}} {qp\frac{{{\rm{d}}p}}{g}} }} $

(5)

式(5)中S为气块以全风速V经过一个中-β尺度距离S(设S=100km)，而q、p分别为气块的比湿和气压。这里，认为T_h是重建中-β尺度湿度(权重)场平均时间。故上述的R_a、R_b可认为是在T_h时间内完成的理想对流性降水速率，则有：

$ {R_A} = \frac{{{R_a}}}{{{T_h}}} $

(3’)

$ {R_B} = \frac{{{R_b}}}{{{T_h}}} $

(4’)

R₃实际取上面R_A与R_B的线性组合，即：

$ {R_3} = {R_1} \cdot {R_A} + {R_2} \cdot {R_B} $

(6)

式(6)预报因子R₃认为，水汽通量散度降水(R₁)同时造成了气块不稳定降水(R_A)，而凝结函数降水(R₂)同时造成了对流不稳定降水(R_B)。又有：

$ \begin{array}{l} {R_4} = E\;\;\;\;\left( {E > 0} \right)\\ {R_4} = - 1\;\;\;\;\;\left( {E = 0} \right) \end{array} $

(7)

式(7)中的E为各个时次模式大气6小时预报降水量。

对于Kalman滤波，须使预报对象降水场R₅与R₄都相对地保持为连续量(场)。故这里的降水场定量预报历史样本选取条件是，只选取各个时次有一半以上预报站点发生降水、且平均降水量大于0.3mm/6h之降水过程为样本。从而有理由设定其它还没有发生降水站点，其降水量为负(-1mm/6h)。即有：

$ \begin{array}{l} {R_5} = {R_s}\;\;\;\;\left( {{R_s} \ge 0} \right)\\ {R_5} = - 1\;\;\;\;\;\left( {无降水} \right) \end{array} $

(8)

上式中的R_s为观测降水量。

在建立定量降水预报方程之前，还先设计用于定性降水预报方程的预报因子。设其4个预报因子分别为S₁、S₂、S₃、S₄和定性降水预报对象(有、无降水) S₅，则有：

$ {S_1}\left( {{R_1} \le 0} \right) = {S_2}\left( {{R_2} \le 0} \right) = {S_3}\left( {{R_3} \le 0} \right) = {S_4}\left( {{R_4} \le 0} \right) = \left( {{R_5} \le 0} \right) = - 1; $ ${S_1}\left( {{R_1} > 0} \right) = {S_2}\left( {{R_2} > 0} \right) = {S_3}\left( {{R_3} > 0} \right) = {S_4}\left( {{R_4} > 0} \right) = \left( {{R_5} > 0} \right) = + 1; $

(9)

其中：

$ {R_0} = \int_{{p_s}}^{{p_T}} {\mathop - \limits^\nabla } V{q_s}\frac{{{\rm{d}}p}}{g}\left( {RH < 80\% } \right) $

(10)

上面(9)式表明，定性降水预报方程的预报因子S₁、S₂、S₃和S₄与预报对象S₅，是将定量降水预报因子R₁、R₂、R₃和R₄与预报对象R₅予以“+1, -1”化，这样就可用Kalman滤波，以零为界，做降水定性“有”(S₅≥0.0)和“无”(S₅＜0)预报。

以上所有降水定量预报因子都为广义速率(mm/h)，而预报对象(降水(mm))为广义距离，故预报因子与预报对象之间是线性关系。使得在降水定量MOS预报方程物理意义上，求得预报因子的(数学期望)系数，便成为各自(降水过程)的统计“滤波”时间。因此，Kalman滤波(或回归)是在“寻找”模式大气中各种天气学意义降水过程的统计最优(或平均)时间。

2.2 地面气温预报因子

地面气温预报方程选取T213L31各个预报时次模式大气中4个预报因子，同样选用在数学上线性化、且其天气学意义明确的预报因子。确定各个预报时次的4个预报因子如下：

T₁：地面气温预报场(℃)；

T₂：预报850hPa层气温(℃，用作次日02时、08时预报)；和表征太阳白昼辐射的T₁₄－T₀₂与T₂₀－T₀₂(℃，用作次日14时、20时预报)；

T₃：预报700~1000hPa厚度场H，用静力平衡方程作垂直积分，计算出由H表示的1000hPa层预报气温(℃)：因

$ H = \frac{R}{{2g}} \cdot \ln \frac{{1000}}{{700}}\left( {{T_3} + {T_{700}}} \right) $

则

$ {T_3} = \frac{{2gH}}{{R\ln \frac{{1000}}{{700}}}} - {T_{700}} $

上面T₇₀₀是700hPa层的预报温度，R是湿空气比气体常数；

T₄：模式大气预报本站气温(℃)，模式大气并未直接预报本站气温，这里设在1000、925hPa层的气温与气压有“T-lnp”近似正比关系^[12]，则有：

$ {T_4} = A\ln {p_0} + B $

上面A、B均为常数，容易求得：

$ A = \left( {{T_{1000}} - {T_{925}}} \right)/\ln \left( {1000/925} \right) $

$ B = \left( {{T_{1000}} + {T_{925}}} \right)/2 - A\ln \left( {1000 \cdot 925} \right)/2 $

又上面p₀为本站气压，同理也用静力平衡方程作垂直积分：设

$ {H_0} = \frac{R}{{2g}} \cdot \ln \frac{{1000}}{{{p_0}}}\left( {{T_{1000}} + {T_{925}}} \right)$

则

$ \ln {p_0} = \ln 1000 - \frac{{2g{H_0}}}{{R\left( {{T_{1000}} + {T_{925}}} \right)}}$

这里，H₀为本站观测场的拔海高度与1000hPa层预报高度的高度差(可为负)，可见是在1000与925hPa层之间，用厚度H₀去内插(或外插)求出p₀及其对应的T₄。

T₅：MOS气温预报场(℃)。

显然，上面的T₁、T₂、T₃、T₄和预报对象T₅的量纲单位都是℃，故预报因子相对于预报对象是线性化的，且其天气学意义明确。这样求得的MOS预报方程，其预报因子系数，正是各自通过“滤波”或回归对于预报对象的权重贡献。

2.3 地面风场预报因子

地面风场预报方程选取T213L31各个预报时次模式大气中4个预报因子，同样应选用在数学意义上线性化、且其天气学意义明确的预报因子。确定各个预报时次的4个预报因子如下：

(u₁, v₁)：1000hPa层预报风场(m/s)；

(u₂, v₂)：预报1000~925hPa层的湍流粘性应力$\left| {\frac{{\partial u}}{{\partial z}}} \right|\frac{{\partial u}}{{\partial z}} $和$\left| {\frac{{\partial v}}{{\partial z}}} \right|\frac{{\partial v}}{{\partial z}} $，其定义是“单位面积空气作用应力”^[15]：由于在近地面层的风向与动量通量随高度变化很小，当气压梯度力、地转偏向力和湍流粘性应力三者平衡时，决定了地面风场的速率与方向，故(u₂, v₂)与地面风速相关；

(u₃, v₃)：预报700~1000hPa层厚度场H，且用天气学地转风(V_T)和热成风(V_g)公式^[15]，计算其成为1000hPa层的地转风形式(m/s)：

因

$ {V_{T700 - 1000}} = \frac{g}{f} \cdot k_\Lambda ^\nabla {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;\;H $

$ {V_{g700}} = \frac{g}{f} \cdot k_\Lambda ^\nabla {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;\;{H_{700}}$

则

$ \left( {{u_3},{v_3}} \right) = {V_{g1000}} = {V_{g700}} - {V_{T700 - 1000}} $

上面H₇₀₀是700hPa层预报高度，f是地转参数；

(u₄, v₄)：预报1000hPa层变压风^[15](V_D，m/s)：

则

$ \left( {{u_4},{v_4}} \right) = {V_{D1000}} = - \frac{g}{{{f^2}}} \cdot k_\Lambda ^\nabla \;\;\;\left( {\Delta H/\Delta t} \right) $

上面Δt=6h，ΔH取其前后两个时次模式大气预报1000hPa层高度差；

(u₅, v₅)：MOS地面风预报场(m/s)。

显然，上面的(u₁, v₁)、(u₂, v₂)、(u₃, v₃)、(u₄, v₄)和预报对象(u₅, v₅)的量纲单位都是m/s，故预报因子相对于预报对象也是线性化的。

3 样本与(相似)Kalman滤波初值参数

Kalman滤波循环递推方程的4个初值参数矩阵${{\tilde \beta }_0}$(数学期望回归系数)、C₀(${{\tilde \beta }_0}$的误差方差)、W(动态噪声方差)和V(量测噪声方差)计算，采用文献[1]中统计回归算法。因T213L31模式逐日提供4个时次(未来30、36、42、48小时)有效预报场，和提供分辨率为3小时预报降水场，“系统”运行时，仅采用各个预报时次历史样本，动态地为当天之前31个(天)，就可每天实时得到各自4个初值矩阵：即求各时次回归方程系数，是分别采用全31个样本和其中的后26个样本，以求出两套样本回归系数之差(Δβ)，同时求出各个Kalman滤波初值参数做预报，叫Kalman滤波；而运用统计相似，将全31个样本时序打乱，按相似(如相似距离、相似系数、或相关系数)程度重新排序，关键的最相似样本排在最后一个，与采用其中的后26个样本(淘汰排前的5个较不相似样本)，同样求出两套样本回归系数之差(Δβ)，并求出各个Kalman滤波初值参数做预报，叫相似-Kalman滤波。故Kalman滤波和相似-Kalman滤波都须用实况资料计算预报误差，用以每天订正系数${{\tilde \beta }_0}$，其重要差别在于后者用最相似样本做最后一次${{\tilde \beta }_0}$修正以作新的预报，而前者默认预报最相似样本即是当天时序最后样本。另外，当两套样本都全用31个相同样本时，则Δβ=0，使得${{\tilde \beta }_0}$≡${{\tilde \beta }_0}$，即滤波过程中${{\tilde \beta }_0}$不变，Kalman滤波就等同于“回归”，即始终是用${{\tilde \beta }_0}$作预报。但次日新样本进入之后，${{\tilde \beta }_0}$总是实时(每日)渐变的。显然，这里使用历史样本是每天更换一个，故上述Kalman滤波和相似-Kalman滤波循环递推方程的4个初值都是每日更新，并都只用于次日一次预报。以上若历史样本足够长，还可选用更多季节相似样本，以实时求得更好的初值参数。所以，初值参数${{\tilde \beta }_0}$、C₀、W和V每日实时确定、其回归性质与随机扰动特性在随季节变化。这里，“系统”仅用了一个月的样本，便可启动全年Kalman滤波和相似-Kalman滤波(或回归)统计MOS预报。

关于统计相似性^[16]：相似距离较好反映样本之间量值相似，相似系数较好反映样本之间形态相似，而相似系数与相关系数较为接近。因此，可以计算各历史样本多个预报因子(场)与当天预报因子(场)的相关系数、相似系数、和相似距离。更可按天气季节适当增加历史样本长度，从中找到真正最优相似样本。且在理论上，相似-Kalman滤波可以消除单纯Kalman滤波所谓“预报滞后”效应^[1-2]，但“相似”本身又将带来新的预报误差，则用相似-Kalman滤波做统计MOS预报，尚有进一步研究意义。

4 “系统”运行

“系统”作自动业务化运行。实际运作是用昨日(20时(北京时))T213L31数值模式预报产品，以定点、定时、定量地预报次日02、08、14、20时地面气温、风场、和前6小时降水。

“系统”运行过程中：

①实时计算和补充所需的各个预报因子与对应(地面气温、风场、和前6小时降雨量)观测资料新样本，形成逐日递补的实时历史样本数据文件。

②滤波初值参数每天更新一次，即31天样本中每天淘汰前一天又补充新一天样本，同时形成新的滤波初值${{\tilde \beta }}$_t做预报。

③实时计算各个预报时次的31个历史样本中的某个预报因子场与当天对应预报场的相关系数(或相似系数)，并可按相似程度将31个历史样本重新排序。

④可用Kalman滤波或相似-Kalman滤波做预报，或用回归做预报(试验表明回归较差)。

⑤分别对降水、地面气温和风场各预报时次31个历史样本的Kalman滤波或相似-Kalman滤波，实时计算它们的平均绝对误差(或均方差)，以做客观检验与效果分析。

⑥若缺当日T213L31数值预报产品某时次资料，则无当日相应时次MOS预报产品；若缺当日预报对象观测资料，则都不能进入历史样本；但都不影响“系统”明日继续运行。

5 “系统”预报试验

“系统”于2005、2006年6—8月为武汉中心气象台做了预报试验。“系统”做次日02、08、14、20时降水预报各为21、32、134、132个站点，而做地面气温和风场预报各为30、30、44、44个站点。现给出评分(参见表1、表2、表 3、表4，限于样本T213L31资料与对应各时次站点观测资料不可或缺，故表中各月都有少量的缺资料日)。

5.1 降水MOS预报

“系统”实际取定性S₂、R₂定量做相似计算，以做降水相似-Kalman滤波预报。

表1是逐日逐时次计算定性预报因子S₂历史样本与当天样本相关系数，做各个时次相似-Kalman滤波降水(定性)预报T_S评分。表1的月平均T_S评分表明，用相似-Kalman滤波预报降水已具一定评分水平。

表2是各个时次做相似-Kalman滤波降水定性预报晴雨月平均准确率。因“系统”预报有/无降水，都给出预报降水发生(雨/晴)概率，是用最新历史样本的预报晴雨准确率代替概率。表2表明，MOS预报可用有一定评分水平的预报晴雨准确率代替预报降水发生概率。

表3是各个时次中等以上降水定量预报T_S评分。表 3表明，“系统”对于6小时中雨、大雨具有一定预报水平，而对于定点(分县)、定时(6小时分辨率)、定量暴雨预报则少有水平。

5.2 地面气温MOS预报

“系统”实际计算样本T₁预报场的相似系数，做相似-Kalman滤波预报次日02时地面气温；和做单纯Kalman滤波预报次日08时、20时地面气温；而做有“条件”相似-Kalman滤波预报次日14时地面气温，即出现当日预报次日地面气温T₁全部站点平均值，相应地比昨日预报今日的平均值有大于2.5℃变化“条件”时，就计算与样本T₁预报场的相关系数，并做有“条件”相似-Kalman滤波，否则仍做单纯Kalman滤波。

这里仅对次日02时和14时(分别接近当天最低、最高气温)预报讨论：做相似-Kalman滤波预报次日02时各站地面气温月平均绝对误差都在2℃左右或以下，而做有“条件”相似-Kalman滤波预报次日14时各站地面气温月平均绝对误差大都在3℃以内(表略)。

5.3 地面风场MOS预报

“系统”实际采用有“条件”相似-Kalman滤波做次日各个时次地面风场预报，即分别出现当日预报次日预报因子u₁、v₁全部站点平均风场，相应地比昨日预报今日有大于1m·s^-1变化“条件”时，就计算各个时次预报风场相关系数，并做有“条件”相似-Kalman滤波，否则仍做单纯Kalman滤波。

用有“条件”相关系数相似-Kalman滤波，预报次日02、08、14和20时各站地面风场u、v分量的月平均绝对误差多在1m/s左右或以内(表略)。而预报各个时次全部站点“平均”地面风场的月平均绝对误差更小(表略)，相当于不定点“面”预报，增加了预报可用性。

6 小结

(1) 用Kalman滤波和相似-Kalman滤波，可作降水、地面气温和风场定点、定时、定量统计MOS预报。本文“系统”经预报试验表明，已经具备一定的定性与定量降水预报T_S评分水平，和一定的地面气温和风场预报评分水平。

(2) MOS预报用相似-Kalman滤波，或有“条件”相似-Kalman滤波，其效果可以好于用单纯Kalman滤波，是因为相似-Kalman滤波消去单纯Kalman滤波“预报滞后”效应，但“相似”预报本身又带来新误差，则各种统计相似之于Kalman滤波仍有进一步研究意义。

(3) MOS预报仅需样本为当天前31个(天)便可启动，其Kalman滤波初值参数可以高频率(每日)予以更新，以便适应自然天气季节变化，但历史样本短，对于局部暴雨、局地雷雨大风等小概率事件的MOS预报不会有水平。

(4) MOS预报可用历史样本雨/晴准确率，代替预报有/无降水发生概率，且有预报意义。

(5) MOS预报地面气温和风场的不定点“面”预报的平均绝对误差比定点预报误差小得多，增加了MOS预报可用性。

(6) MOS预报不出现奇异值，与选取对于预报对象是线性化因子有关，但线性化因子之间的独立性相对较差，将来或可考虑参用其它模式预报产品与因子。

(7) MOS预报水平随季节，实际随不同季节内的不同天气过程有所不同，天气越稳定，预报效果越好，表明有待增加/更换季节相似样本，或天气过程相似样本，并相应地对各种相似-Kalman滤波做深入研究。

(8) 本文“系统”容易移植到其他省、地气象台。