台风多变量集合敏感性分析方法、台风预报方法及其系统

摘要

本发明公开了台风多变量集合敏感性分析方法、台风预报方法及其系统。所述台风多变量集合敏感性分析方法，基于台风的集合预报，以集合的预报响应函数和分析场的影响因子建立多元线性回归模型，计算多变量回归系数。当给定分析场的模式某格点某变量的扰动时，通过该格点该变量与影响因子中除去该格点该变量的其他变量的相关性计算分析场的增量，并以局地化函数进行约束，再结合多变量回归系数，估计得到预报响应函数的变化。因此，本发明可以简单、高效地指出影响台风预报的初始关键区域和变量，有效提高台风的敏感性估计和预报水平。

权利要求

1.一种台风多变量集合敏感性的分析方法，其特征在于，包括以下步骤：

（1）确定预报响应函数J及其影响因子x

给定一个台风的集合预报，于引导时间(t1，t2)内，比较台风的集合预报结果和台风的实际观测数值，得到台风的预报误差；其中，t1时刻指集合预报的给定分析场的分析时刻；t2时刻指集合预报给定分析场的预报时刻；

根据预报误差，构建预报响应函数J(x)：预报响应函数J(x)能够表征集合预报结果，且在引导时间(t1，t2)内，预报响应函数J(x)与预报响应函数的影响因子x之间的关系符合多元线性回归模型；影响因子x为集合预报初始场中对集合预报结果有影响的模式各格点上的各物理量，影响因子x有P个物理量，分别为x1,x2,…,xP,P为正整数；集合预报初始场指分析时刻t1对应的集合预报给定的分析场；

影响因子x选定为模式区域内各模式高度层上各模式格点的位温T、水汽混合比Q、切向风TW、径向风RW和气柱扰动干空气质量MU；

（2）获取预报响应函数J对影响因子x的多变量集合敏感性

在集合预报给定的分析场中，基于预报响应函数J(x)与影响因子x之间的多元线性回归模型，结合t1时刻对应的集合预报初始场的影响因子x和影响因子x在t2时刻对应的集合预报结果J(x)，计算出该多元线性回归模型的多变量回归系数，从而表征预报响应函数J(x)对影响因子x的多变量集合敏感性；

（3）获取预报响应函数变化的估计

在t1时刻，针对集合预报初始场的影响因子x中位于模式某格点A位置处的某变量xp，给定扰动σp或增量δxp时，通过格点A处的该变量xp与影响因子x中除该变量xp外的其余变量之间的相关性，并以局地化函数进行约束，获得t1时刻的集合预报初始场增量δx；再结合步骤(2)中计算出的多元线性回归模型的多变量回归系数，得到对预报响应函数变化的估计。

2.根据权利要求1所述的台风多变量集合敏感性分析方法，其特征在于，预报响应函数J与影响因子x之间的多元线性回归模型为：J＝b1x1+b2x2+…+bPxP+ε           (1)

其中：x1,x2,…,xP为P个分析时刻t1的影响因子；b1,b2,…,bP为多元回归系数，表征预报响应函数J对各影响因子x的敏感性；ε为多元线性回归模型的残差。

3.根据权利要求2所述的台风多变量集合敏感性分析方法，其特征在于，当集合预报在分析时刻t1的某影响因子xp发生扰动σp时，多变量集合敏感性估计的预报响应函数的变化δJ根据下式计算：式中：ρ表示局地化函数，表示Hadamard积，δx是由δx1,δx2,...,δxP构成的P×1矩阵；

公式(2)中增量δx的表达式为：

其中，上标e表示由K个集合成员各自距平的xi构成的1×K矩阵，cov为协方差，var为方差，i表示影响因子序号；

公式(2)中多变量回归系数β的表达式为：

其中，β为多元回归系数b1,b2,...,bP组成的P×1矩阵；x为影响因子x1,x2,…,xP组成的P×1矩阵；X为集合扰动 构成的P×K矩阵；J为K个集合成员的预报响应函数构成的K×1矩阵；上标T表示矩阵转置，上标-1表示矩阵求逆。

4.一种基于台风多变量集合敏感性分析的台风预报方法，其特征在于，包括以下步骤：（1）确定预报响应函数J及其影响因子x

在给定分析场中，于预报时长(t1，t2)内，比较台风的集合预报和台风的实际观测，得到台风的预报误差；其中，t1时刻指集合预报的给定分析场的分析时刻；t2时刻指集合预报给定分析场的预报时刻；

根据预报误差，构建预报响应函数J(x)：预报响应函数J(x)能够表征集合预报结果，且在预报时长(t1，t2)内，预报响应函数J(x)与预报响应函数的影响因子x之间的关系符合多元线性回归模型；影响因子x为集合预报初始场中对集合预报结果有影响的模式各格点上的各物理量，影响因子x有P个物理量，分别为x1,x2,…,xP,P为正整数；集合预报初始场指分析时刻t1对应的集合预报给定的分析场；

影响因子x选定为模式区域内各模式高度层上各模式格点的位温T、水汽混合比Q、切向风TW、径向风RW和气柱扰动干空气质量MU；

（2）获取预报响应函数J对影响因子x的多变量集合敏感性

在集合预报给定的分析场中，基于预报响应函数J(x)与影响因子x之间的多元线性回归模型，结合t1时刻对应的集合预报初始场的影响因子x和影响因子x在t2时刻对应的集合预报结果J(x)，计算出该多元线性回归模型的多变量回归系数，从而表征预报响应函数J(x)对影响因子x的多变量集合敏感性；

（3）获取预报响应函数变化的估计

在t1时刻，针对集合预报初始场的影响因子x中位于模式某格点A位置处的某变量xp，给定扰动σp或增量δxp时，通过格点A处的该变量xp与影响因子x中除该变量xp外的其余变量之间的相关性，并以局地化函数进行约束，获得t1时刻的集合预报初始场增量δx；再结合步骤(2)中计算出的多元线性回归模型的多变量回归系数，得到对预报响应函数变化的估计；

（4）台风预报

根据步骤(3)所得到的预报响应函数变化的估计，找到预报时长(t1，t2)内，对预报结果影响最大的关键区域和变量；将给定此关键变量扰动σp或者增量δxp带来的相应的集合预报初始场增量δx与t1时刻的影响因子x结合，得到t1时刻改进的集合预报初始场；

从该改进的集合预报初始场出发进行台风的数值预报，即可获得改进的台风数值预报。

5.根据权利要求4所述的基于台风多变量集合敏感性分析的台风预报方法，其特征在于，预报响应函数J与影响因子x之间的多元线性回归模型为：J＝b1x1+b2x2+…+bPxP+ε     (5)

其中：x1,x2,…,xP为P个分析时刻t1的影响因子；b1,b2,...,bP为多元回归系数，表征预报响应函数J对各影响因子x的敏感性；ε为多元线性回归模型的残差。

6.根据权利要求5所述的基于台风多变量集合敏感性分析的台风预报方法，其特征在于，步骤(3)中，当集合预报在分析时刻t1的某影响因子xp发生扰动σp时，多变量集合敏感性估计的预报响应函数的变化δJ根据下式计算：式中：ρ表示局地化函数，表示Hadamard积，δx是由δx1,δx2,...,δxP构成的P×1矩阵；

公式(6)中增量δx的表达式为：

其中，上标e表示由K个集合成员各自距平的xi构成的1×K矩阵，cov为协方差，var为方差，i表示影响因子序号；

公式(6)中多变量回归系数β的表达式为：

其中，β为多元回归系数b1,b2,…,bP组成的P×1矩阵；x为影响因子x1,x2,…,xP组成的P×1矩阵；X为集合扰动 构成的P×K矩阵；J为K个集合成员的预报响应函数构成的K×1矩阵；上标T表示矩阵转置，上标-1表示矩阵求逆。

7.一种基于台风多变量集合敏感性分析系统，其特征在于，包括中央处理单元，该中央处理单元中运行有第一计算机程序，该第一计算机程序可被执行以实现如权利要求1-3中任一项所述的方法。

8.一种基于台风多变量集合敏感性分析的台风预报系统，其特征在于，包括中央处理单元，该中央处理单元中运行有第二计算机程序，该第二计算机程序可被执行以实现如权利要求4-6中任一项所述的方法。

9.一种计算机可读介质，其特征在于，存储有第一计算机程序，该第一计算机程序可被执行以实现如权利要求1-3中任一项所述的方法。

10.一种计算机可读介质，其特征在于，存储有第二计算机程序，该第二计算机程序可被执行以实现如权利要求4-6中任一项所述的方法。

说明书

台风多变量集合敏感性分析方法、台风预报方法及其系统

技术领域

[0001] 本发明涉及一种台风多变量集合敏感性分析方法及其系统，属于数值天气预报领域。

[0002] 本发明还涉及一种基于台风多变量集合敏感性分析的台风预报方法及其系统。

背景技术

[0003] 在数值天气预报中，敏感性分析可以定量估计初始场扰动对预报结果的影响。通过敏感性分析，可以找到对预报结果影响最大的关键变量及区域，进而提高预报能力，并增强对天气系统误差传播的理解。

[0004] 一种常用的敏感性分析方法，是在主观选择的相对重要的初始场变量及区域添加扰动，通过积分数值模式，直接获得预报结果的改变。这种方式的局限在于，人工判断关键变量和区域比较困难，并且需要花费大量的计算资源用于数值模拟。

[0005] 相对于主观的敏感性分析方法，客观的敏感性方法不需要反复积分数值模式，因此可以满足敏感性分析的时效，并节约大量的计算资源。客观的敏感性分析方法主要包括伴随敏感性和集合敏感性。其中：

[0006] 伴随敏感性，采用非线性数值模式的切线性模式及其伴随模式，通过反向积分获得给定预报改变所需要的初始场扰动。目前国内外已有不少研究将伴随敏感性应用于提高台风预报的准确性。但不足的是，构建非线性数值模式的切线性模式及其伴随模式非常困难，特别是在边界层和微物理参数化的“开关”变量等问题上。

[0007] 集合敏感性，利用集合预报建立回归模型，计算回归系数，即可简单、高效地获得预报结果对初始场的敏感性。常用的集合敏感性分析为单变量集合敏感性，其计算模式某一格点上某一变量的改变带来的预报结果的变化，没有考虑同一格点上不同变量间的相互作用以及相邻格点间不同变量间的相互作用。Hacker andLei(2015)指出，基于一元线性回归的单变量集合敏感性忽略影响因子之间的相互作用，往往会导致敏感性高估，特别是在中小尺度强非线性天气过程中。由此，Hacker and Lei(2015)提出基于多元线性回归的多变量集合敏感性，这种方法在Lorenz(2005)理想模式中被证明可以提供更准确的敏感性估计，但缺乏真实天气过程的应用。

发明内容

[0008] 本发明将多变量集合敏感性技术首次应用于台风预报的初始敏感性分析中，提出了一种台风多变量集合敏感性分析的台风预报方法。该方法基于台风的集合预报，以集合的初始场与集合预报响应函数建立多元线性回归模型，计算多变量回归系数。当给定初始场模式某格点某变量的增量时，通过该格点该变量与该格点的其他变量和相邻格点的各变量的相关性，并受到局地化函数的约束，再结合多变量回归系数，估计得到预报响应函数的变化。因此，本发明可以简单、高效地指出影响台风预报的关键区域和变量，有效提高台风的敏感性估计和预报水平。

[0009] 为实现上述的技术目的，本发明将采取如下的技术方案：

[0010] 一种台风多变量集合敏感性的分析方法，其特征在于，包括以下步骤：

[0011] (1)确定预报响应函数J及其影响因子x

[0012] 给定一个台风的集合预报，于引导时间(t1，t2)内，比较台风的集合预报结果和台风的实际观测数值，得到台风的预报误差；其中，t1时刻指集合预报的给定分析场的分析时刻；t2时刻指集合预报给定分析场的预报时刻；

[0013] 根据预报误差，构建预报响应函数J(x)：预报响应函数J(x)能够表征集合预报结果，且在引导时间(t1，t2)内，预报响应函数J(x)与预报响应函数的影响因子x之间的关系符合多元线性回归模型；影响因子x为集合预报初始场中对集合预报结果有影响的模式各格点上的各物理量，影响因子x有P个物理量，分别为x1,x2,…,xP,P为正整数；集合预报初始场指分析时刻t1对应的集合预报给定的分析场；

[0014] (2)获取预报响应函数J对影响因子x的多变量集合敏感性

[0015] 在集合预报给定的分析场中，基于预报响应函数J(x)与影响因子x之间的多元线性回归模型，结合t1时刻对应的集合预报初始场的影响因子x和影响因子x在t2时刻对应的集合预报结果J(x)，计算出该多元线性回归模型的多变量回归系数，从而表征预报响应函数J(x)对影响因子x的多变量集合敏感性；

[0016] (3)获取预报响应函数变化的估计

[0017] 在t1时刻，针对集合预报初始场的影响因子x中位于模式某格点A位置处的某变量xp，给定扰动σp或增量δxp时，通过格点A处的该变量xp与影响因子x中除该变量xp外的其余变量之间的相关性，并以局地化函数进行约束，获得t1时刻的集合预报初始场增量δx；再结合步骤(2)中计算出的多元线性回归模型的多变量回归系数，得到对预报响应函数变化的估计。

[0018] 本发明的再一个技术目的是提供一种基于台风多变量集合敏感性分析的台风预报方法，包括以下步骤：

[0019] (1)确定预报响应函数J及其影响因子x

[0020] 在给定分析场中，于预报时长(t1，t2)内，比较台风的集合预报和台风的实际观测，得到台风的预报误差；其中，t1时刻指集合预报的给定分析场的分析时刻；t2时刻指集合预报给定分析场的预报时刻；

[0021] 根据预报误差，构建预报响应函数J(x)：预报响应函数J(x)能够表征集合预报结果，且在预报时长(t1，t2)内，预报响应函数J(x)与预报响应函数的影响因子x之间的关系符合多元线性回归模型；影响因子x为集合预报初始场中对集合预报结果有影响的模式各格点上的各物理量，影响因子x有P个物理量，分别为x1,x2,…,xP,P为正整数；集合预报初始场指分析时刻t1对应的集合预报给定的分析场；

[0022] (2)获取预报响应函数J对影响因子x的多变量集合敏感性

[0023] 在集合预报给定的分析场中，基于预报响应函数J(x)与影响因子x之间的多元线性回归模型，结合t1时刻对应的集合预报初始场的影响因子x和影响因子x在t2时刻对应的集合预报结果J(x)，计算出该多元线性回归模型的多变量回归系数，从而表征预报响应函数J(x)对影响因子x的多变量集合敏感性；

[0024] (3)获取预报响应函数变化的估计

[0025] 在t1时刻，针对集合预报初始场的影响因子x中位于模式某格点A位置处的某变量xp，给定扰动σp或增量δxp时，通过格点A处的该变量xp与影响因子x中除该变量xp外的其余变量之间的相关性，并以局地化函数进行约束，获得t1时刻的集合预报初始场增量δx；再结合步骤(2)中计算出的多元线性回归模型的多变量回归系数，得到对预报响应函数变化的估计。

[0026] (4)台风预报

[0027] 根据步骤(3)所得到的预报响应函数变化的估计，找到预报时长(t1，t2)内，对预报结果影响最大的关键区域和变量；将给定此关键变量扰动σp或者增量δxp带来的相应的集合预报初始场增量δx与t1时刻的影响因子x结合，得到t1时刻改进的集合预报初始场；

[0028] 从该改进的集合预报初始场出发进行台风的数值预报，即可获得改进的台风数值预报。

[0029] 一种基于台风多变量集合敏感性分析系统，包括中央处理单元，该中央处理单元中运行有第一计算机程序，该第一计算机程序可被执行以实现上述的台风多变量集合敏感性的分析方法。

[0030] 一种基于台风多变量集合敏感性分析的台风预报系统，包括中央处理单元，该中央处理单元中运行有第二计算机程序，该第二计算机程序可被执行以实现如上述的基于台风多变量集合敏感性分析的台风预报方法。

[0031] 一种计算机可读介质，，存储有第一计算机程序，该第一计算机程序可被执行以实现上述的台风多变量集合敏感性的分析方法。

[0032] 一种计算机可读介质，存储有第二计算机程序，该第二计算机程序可被执行以实现上述的基于台风多变量集合敏感性分析的台风预报方法。

[0033] 根据上述的技术方案，相对于现有技术，本发明具有如下的有益效果：

[0034] 本发明将多变量集合敏感性首次运用于台风初始敏感性分析中，相比主观选择分析而言节约了大量计算资源，相比伴随敏感性更便于理解和操作、并且不需要建立数值模式的线性模式和伴随模式，相比单变量集合敏感性则有效解决了其无法精确估计敏感性区域和变量的关键问题，因此可以高效地、准确地获得敏感性的估计。

[0035] 通过验证，经本发明所述的预报方法获得的台风数值预报结果，与实际的观测值相比，更为准确。

附图说明

[0036] 图1为本发明所述台风多变量集合敏感性分析方法的流程示意图。

[0037] 图2为超强台风“海燕”集合预报在2013年11月4日00时的三层嵌套网格示意图。

[0038] 图3为超强台风“海燕”从2013年11月4日00时起，历时126小时的路径观测值(虚线)和集合预报结果(细实线：各集合成员，粗实线：集合平均)。

[0039] 图4为超强台风“海燕”从2013年11月4日00时起，历时126小时的台风中心海平面气压(Sea level pressure,单位hPa)观测值(粗虚线)和集合预报结果(细实线：各集合成员，粗实线：集合平均)；细虚线标出了选择分析的预报响应函数时刻。

[0040] 图5为扰动在预报响应函数前24小时的位温所得的基于集合敏感性分析的预报响应函数的变化(单位hPa)；其中，图5a-图5c展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过单变量集合敏感性分析得到的结果；图5d-图5f展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过多变量集合敏感性分析得到的结果。

[0041] 图6为扰动在预报响应函数前24小时的水汽混合比所得的基于集合敏感性分析的预报响应函数的变化(单位hPa)；其中：图6a、图6b展示的分别为添加扰动到最靠近850hPa、500hPa模式层时，单变量集合敏感性分析得到的结果；图6c、图6d展示的分别为添加扰动到最靠近850hPa、500hPa模式层时，通过多变量集合敏感性分析得到的结果。

[0042] 图7为扰动在预报响应函数前24小时的切向风所得的基于集合敏感性分析的预报响应函数的变化(单位hPa)；其中，图7a-图7c展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过单变量集合敏感性分析得到的结果；图7d-图7f展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过多变量集合敏感性分析得到的结果。

[0043] 图8为扰动在预报响应函数前24小时的径向风所得的基于集合敏感性分析的预报响应函数的变化(单位hPa)；其中：图8a-图8c展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过单变量集合敏感性分析得到的结果；图8d-图8f展示的分别为添加扰动到500-600km环状区域、气压平均最靠近850hPa、500hPa、200hPa模式层时，通过多变量集合敏感性分析得到的结果。

[0044] 图9a展示的为扰动在预报响应函数前48小时的850hPa位温(T)所得的基于单变量集合敏感性分析的预报响应函数的变化(单位hPa)；而图9e展示的为扰动在预报响应函数前48小时的850hPa位温(T)所得的基于多变量集合敏感性分析的预报响应函数的变化(单位hPa)；

[0045] 图9b展示的为扰动在预报响应函数前48小时的850hPa水汽混合比(Q)所得的基于单变量集合敏感性分析的预报响应函数的变化(单位hPa)；而图9f展示的为扰动在预报响应函数前48小时前的850hPa水汽混合比(Q)所得的基于多变量集合敏感性分析的预报响应函数的变化(单位hPa)；

[0046] 图9c展示的为扰动在预报响应函数前248小时的850hPa切向风(TW)所得的基于单变量集合敏感性分析的预报响应函数的变化(单位hPa)；图9g展示的为扰动在预报响应函数前48小时的850hPa切向风(TW)所得的基于多变量集合敏感性分析的预报响应函数的变化(单位hPa)；

[0047] 图9d展示的为扰动在预报响应函数前48小时的850hPa径向风(RW)所得的基于单变量集合敏感性分析的预报响应函数的变化(单位hPa)；图9h展示的为扰动在预报响应函数前48小时的850hPa径向风(RW)所得的基于多变量集合敏感性分析的预报响应函数的变化(单位hPa)。

[0048] 图10a展示的为扰动模式初始场获得的真实预报响应函数的变化与单变量集合敏感性估计值的对比图(引导时间24小时)；

[0049] 图10b展示的为扰动模式初始场中获得的真实预报响应函数的变化与单变量集合敏感性估计值的对比图(引导时间48小时)；

[0050] 图10c展示的为扰动模式初始场中获得的真实预报响应函数的变化与多变量集合敏感性估计值的对比图(引导时间24小时)；

[0051] 图10d展示的为扰动模式初始场中获得的真实预报响应函数的变化与多变量集合敏感性估计值的对比图(引导时间48小时)；

[0052] 图10a-图10d中，虚线为最小二乘法回归线，RMSE是真值与估计值的均方根误差。

具体实施方式

[0053] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、表达式和数值不限制本发明的范围。同时，应当明白，为了便于描述，附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。

[0054] 实施例1

[0055] 图1是一种计算机程序的流程框图，其中，不包含虚线框，即仅包括实线框的部分，为本发明所述台风多变量集合敏感性分析方法的流程框图；以下将结合实线框对应的流程图部分，详细地说明本发明所述的台风多变量集合敏感性分析方法，包括以下步骤：

[0056] (1)确定预报响应函数J及其影响因子x

[0057] 1.1、台风集合预报及观测对比

[0058] 对目标台风进行集合预报，获得目标台风的集合预报数据，该集合预报数据包括K个集合成员；

[0059] 获取目标台风的观测数据，包括台风位置和强度等，将目标台风的集合预报数据与观测数据进行对比分析，获得目标台风集合预报的误差；

[0060] 1.2确定预报响应函数和影响因子

[0061] 根据步骤1.1得到的预报误差，首先选定能够表征该预报结果的预报响应函数；例如在目标时刻代表台风强度预报误差的台风中心海平面气压的误差、或者最大风速的误差，或者代表台风路径预报误差的台风中心偏移距离等。

[0062] 然后选择集合预报初始场中对集合预报结果有影响的各模式格点上的各物理量作为影响因子，例如每一模式格点的位温、水汽混合比、经、纬向风、海洋表面温度等。

[0063] (2)获取预报响应函数J对影响因子x的多变量集合敏感性

[0064] 2.1、构建预报响应函数J与影响因子x之间的多元线性回归模型

[0065] 在48小时预报时间内，可以认为集合预报初始场扰动与集合预报的结果满足线性关系，P个影响因子x1,x2,…,xP和预报响应函数J之间的关系符合多元线性回归模型：

[0066] J＝b1x1+b2x2+…+bPxP+ε

[0067] 其中：x1,x2,…,xP为P个影响因子；b1,b2,…,bP为回归系数，表征预报响应函数对各影响因子的敏感性；ε为多元线性回归模型的残差；

[0068] 2.2、根据预报响应函数J与影响因子x之间的多元线性回归模型，计算回归系数b1,b2,…,bP组成的矩阵β

[0069] 当拥有K个集的初始场和其对应的集合预报时，相当于拥有了计算多元线性回归模型回归系数的K个样本，回归系数可解。将x1,x2,…,xP写作P×1矩阵x，将集合距平的K个x收集进入P×K矩阵X中，将K个集合距平的J收集进入K×1矩阵J中，回归系数b1,b2,…,bP组成的矩阵β计算公式如下：

[0070]

[0071] 其中：上标T表示矩阵转置，上标-1表示矩阵求逆；

[0072] (3)估计初始场扰动带来的预报结果变化

[0073] 当集合预报初始场某影响因子xp发生扰动σp(或者增量δxp，若为增量δxp，所涉及公式中与扰动σp对应的部分均修改为δxp)时，多变量集合敏感性估计得到的预报响应函数变化δJ根据下式计算：

[0074]

[0075] 式中：ρ表示局地化函数，表示Hadamard积，δx是由δx1,δx2,…,δxP构成的P×1矩阵；

[0076] 其中：

[0077]

[0078] 其中，上标e表示由K个集合成员各自距平的xi构成的1×K矩阵，cov为协方差，var为方差，i表示影响因子序号。

[0079] 还需要特别说明的是，预报响应函数的变化δJ的表达式不唯一，比如也可以不包含上述的局地化函数，可使用近似的矩阵转置表达式来替代。增量δx的表达式不唯一，比如也可以使用矩阵的表达式或矩阵求逆的表达式等。

[0080] 本发明还提供一种基于台风多变量集合敏感性分析系统，包括中央处理单元，该中央处理单元中运行有第一计算机程序，该第一计算机程序可被执行以实现本实施例1所述的台风多变量集合敏感性的分析方法。

[0081] 本发明还提供一种计算机可读介质，，存储有第一计算机程序，该第一计算机程序可被执行以实现实施例1所述的台风多变量集合敏感性的分析方法。

[0082] 实施例2

[0083] 附图1中，实线框与虚线框一起构成的整个流程图，公开了本发明所述的一种基于台风多变量集合敏感性分析的台风预报方法，其与实施例1相比，具有相同的步骤(1)至(3)，不同仅在于，在实施1中步骤(3)施行后，还有一个步骤(4)，具体地，本发明所述的台风预报方法，包括以下步骤：

[0084] (1)确定预报响应函数J及其影响因子x

[0085] 在给定分析场中，于预报时长(t1，t2)内，比较台风的集合预报和台风的实际观测，得到台风的预报误差；其中，t1时刻指集合预报的给定分析场的分析时刻；t2时刻指集合预报给定分析场的预报时刻；

[0086] 根据预报误差，构建预报响应函数J(x)：预报响应函数J(x)能够表征集合预报结果，且在预报时长(t1，t2)内，预报响应函数J(x)与预报响应函数的影响因子x之间的关系符合多元线性回归模型；影响因子x为集合预报初始场中对集合预报结果有影响的模式各格点上的各物理量，影响因子x有P个物理量，分别为x1，x2，…，xP，P为正整数；集合预报初始场指分析时刻t1对应的集合预报给定的分析场；其中：预报响应函数J与影响因子x之间的多元线性回归模型为：

[0087] J＝b1x1+b2x2+…+bPxP+ε    (5)

[0088] 其中：x1，x2，…，xP为P个分析时刻t1的影响因子；b1，b2，…，bP为多元回归系数，表征预报响应函数J对各影响因子x的敏感性；ε为多元线性回归模型的残差。

[0089] (2)获取预报响应函数J对影响因子x的多变量集合敏感性

[0090] 在集合预报给定的分析场中，基于预报响应函数J(x)与影响因子x之间的多元线性回归模型，结合t1时刻对应的集合预报初始场的影响因子x和影响因子x在t2时刻对应的集合预报结果J(x)，计算出该多元线性回归模型的多变量回归系数，从而表征预报响应函数J(x)对影响因子x的多变量集合敏感性。

[0091] (3)获取预报响应函数变化的估计

[0092] 在t1时刻，针对集合预报初始场的影响因子x中位于模式某格点A位置处的某变量xp，给定扰动或增量δxp时，通过格点A处的该变量xp与影响因子x中除该变量xp外的其余变量之间的相关性，并以局地化函数进行约束，获得t1时刻的集合预报初始场增量δx；再结合步骤(2)中计算出的多元线性回归模型的多变量回归系数，得到对预报响应函数变化的估计。

[0093] 具体地，当集合预报在分析时刻t1的某影响因子xp发生扰动σp(或者给定增量δxp)时，多变量集合敏感性估计的预报响应函数的变化δJ根据下式计算：

[0094]

[0095] 式中：p表示局地化函数，表示Hadamard积，δx是由δx1，δx2，...，δxP构成的P×1矩阵；

[0096] 公式(6)中增量δx的表达式为：

[0097]

[0098] 其中，上标e表示由K个集合成员各自距平的xi构成的1×K矩阵，cov为协方差，var为方差，i表示影响因子序号；

[0099] 公式(6)中多变量回归系数β的表达式为：

[0100]

[0101] 其中，β为多元回归系数b1，b2，…，bP组成的P×1矩阵；x为影响因子x1，x2，…，xP组成的P×1矩阵；x为集合扰动 构成的P×K矩阵；J为K个集合成员的预报响应函数构成的K×1矩阵；上标T表示矩阵转置，上标-1表示矩阵求逆。

[0102] (4)台风预报

[0103] 根据步骤(3)所得到的预报响应函数变化的估计，找到预报时长(t1，t2)内，对预报结果影响最大的关键区域和变量；将给定此关键变量扰动σp或者增量δxp带来的相应的集合预报初始场增量δx与t1时刻的影响因子x结合，得到t1时刻改进的集合预报初始场；

[0104] 从该改进的集合预报初始场出发进行台风的数值预报，即可获得改进的台风数值预报。

[0105] 将根据步骤(4)得到的改进的预报结果，与实际的台风观测比较，可以进一步验证本发明所述的多变量集合敏感性的可靠性。

[0106] 本发明还提供一种基于台风多变量集合敏感性分析的台风预报系统，包括中央处理单元，该中央处理单元中运行有第二计算机程序，该第二计算机程序可被执行以实现如实施例2所述的基于台风多变量集合敏感性分析的台风预报方法。

[0107] 本发明还提供一种计算机可读介质，存储有第二计算机程序，该第二计算机程序可被执行以实现实施例2所述的基于台风多变量集合敏感性分析的台风预报方法。

[0108] 应用例

[0109] 超强台风“海燕”是2013年全球最强的台风，截至目前仍保持着西北太平洋台风的最大风速记录。根据美国国家飓风中心热带气旋关键数据集(Tropical Cyclone Vitals，TCVitals)，“海燕”于2013年11月4日06时(协调世界时，下同)命名，台风中心位于6.0°N，150.2°E。随后，其快速增强并向西北移动，于11月7日18时到达峰值强度，中心气压仅
895hPa，最大风速达87m/s。数小时后，“海燕”登陆菲律宾，造成当地约36000人员伤亡失踪和巨额经济损失。

[0110] 在中央气象台业务预报中，“海燕”的路径预报与实际观测比较符合，但最大强度预报存在明显低估。通过敏感性分析提高“海燕”强度预报准确性，对理解台风动力结构、提高数值预报能力、增强人民生命财产安全保障都具有重大的意义。

[0111] 因此，将本发明所述的多变量集合敏感性分析方法应用于台风“海燕”强度预报的敏感性分析中，并与传统的单变量集合敏感性进行对比，具体如下：

[0112] 1.台风集合预报及观测对比

[0113] 台风集合预报采用美国天气研究与预报模式高级研究版本(The Advanced Research Weather Research and Forecasting Model，WRF-ARW)3.4版本。

[0114] 图2为超强台风“海燕”集合预报的三层嵌套网格设置示意图。最外层网格(Domain 01,D01)固定，内侧网格(D02,D03)随台风中心移动。D01、D02、D03网格数分别为320×270、
198×198、360×360，水平分辨率分别为27km，9km和3km。垂直设置56个模式层，模式层顶为
10hPa。

[0115] 物理参数化方面，本示例采用快速辐射传输模式(Rapid Radiative Transfer Model，RRTM)长波和短波辐射方案、统一Noah陆面模式(The Unified Noah land-surface model)、延世大学行星边界层方案(The Yonsei University planetary boundary layer scheme)、WRF单参数(WRF single-moment,WSM)六种粒子微物理方案(6-class microphysics scheme)。此外，积云参数化方案仅在D01中使用，采用修正的Tiedtke积云方案(The modified Tiedtke cumulus scheme)。

[0116] 对于初始条件和边界条件，美国国家环境预测中心(National Center for Environmental Prediction，NCEP)的全球预报系统(Global Forecast System,GFS)发布的0.25°×0.25°的分析场资料经过插值和扰动提供每6小时一次的集合边界条件和11月1日00时的初始条件。扰动方式基于固定协方差扰动技术(Fixed-covariance perturbation technique)，利用WRF资料同化系统三维变分同化(WRFDA-3DVAR)对气候态的背景误差协方差进行随机采样，产生满足气候态背景误差协方差的随机扰动。后续时刻的初始场来自于循环的集合卡尔曼滤波系统(Ensemble Kalman filter system)。

[0117] 在11月4日00时进行长时间的集合预报，获得从11月4日00时至11月9日06时共126小时的集合预报，其中涵盖了台风“海燕”的强度峰值时刻。集合成员共80个。

[0118] 观测资料使用的是TCVitals数据集。

[0119] 图3和图4分别为“海燕”2013年11月4日00时起126小时的路径和最低海平面气压的观测值(粗虚线)及集合预报(细实线：集合成员，粗实线：集合平均)。可以看出，集合预报对路径的模拟略有偏差，但整体还是比较准确的。但强度预报明显偏弱，集合成员之间的差异也较大，表明“海燕”的强度预报存在较大的困难。

[0120] 2.确定预报响应函数和影响因子

[0121] 为分析台风强度预报于初始场的敏感性，选择在靠近台风峰值强度的11月8日00时(图4中细虚线指示)，代表强度预报误差的台风中心最低海平面气压误差(预报值减观测值)作为预报响应函数。选择分别在11月7日00时和11月6日00时，模式区域D03内各模式高度层上各模式格点的位温(T)、水汽混合比(Q)、切向风(TW，逆时针环流为正)、径向风(RW，出流为正)和气柱扰动干空气质量(MU)作为影响因子，共计360×360×56×4+360×360＝29160000个。其中，TW和RW由D03模式输出的经纬向风插值到T所在格点获得。另外，影响因子的引导时间分别为24h和48h。

[0122] 3.计算敏感性

[0123] 按公式计算敏感性：

[0124] 4.估计初始场扰动带来的预报结果改变

[0125] 给初始场影响因子一定扰动，根据 可以算得其扰动带来的预报响应函数变化。

[0126] 在本示例中，4.1中扰动的影响因子是在目标时刻11月8日00时前24、48小时，代表台风低层(850hPa)、中层(500hPa)、高层(200hPa)的模式垂直层上纬向每五个格点的T、TW、RW等变量，以及850hPa、500hPa模式垂直层上每五个格点的Q变量，共计2×360/5×360×(3×3+2)＝570240次扰动试验。4.2中扰动影响因子是在各引导时间时，在950hPa、850hPa、700hPa、500hPa、200hPa、100hPa模式垂直层上的T、TW、RW等变量，以及950hPa、850hPa、
700hPa、500hPa模式垂直层上的Q变量中随机选择的100个影响因子；共计2×100＝200次扰动试验。

[0127] 扰动的大小为该影响因子的集合离散度(其集合预报的标准差)。采用的ρ是常用的Gaspariand Cohn(GC)局地化函数，局地化尺度水平为2000km，垂直为1.5ln(hPa)。

[0128] 4.1与单变量集合敏感性对比

[0129] 为了阐明多变量集合敏感性的优势，首先将其与单变量集合敏感性进行对比。

[0130] 图5为扰动在目标时刻的24小时前的位温时，通过(a)-(c)单变量集合敏感性和(d)-(f)多变量集合敏感性估计得到的预报响应函数的变化(单位hPa)。第1到3行分别为添加扰动到500-600km环状区域气压平均最靠近850hPa、500hPa、200hPa的模式层的结果。注意，单变量与多变量集合敏感性的填色范围是不同的。

[0131] 可以看出，多变量集合敏感性的形态分布与单变量是相似的。负大值区是以台风中心为圆心、半径几百到近千公里随高度上升而变宽的近似圆形的区域，这表示该区域在目标时刻24小时前，位温正扰动会带来最低海平面气压预报误差的下降，台风强度预报的增强。这与台风暖心增强将会带来后续强度增强的特征一致，说明该结果在物理上是合理的。

[0132] 但是单变量估计的敏感性明显高于多变量，特别是在对流层高层。例如，在单变量集合敏感性估计中，850hPa模式高度层上台风中心位置添加该点一倍集合离散度扰动最高可以带来-12.6hPa预报响应函数的变化，但多变量集合敏感性对相同扰动的相应估计值仅为-10.8hPa。到500hPa模式高度层上，单变量和多变量集合敏感性估计的极值则分别变为了-12.1hPa和-2.9hPa。到200hPa模式高度层上，单变量和多变量集合敏感性的差异进一步拉大，两种方法的敏感性估计值分别变为了-13.3hPa和-0.4hPa。

[0133] 图6为扰动在目标时刻24小时前的水汽混合比，通过(a)-(b)单变量集合敏感性和(c)-(d)多变量集合敏感性分别得到的预报响应函数变化(单位hPa)。第1、2行分别展示添加扰动到最靠近850hPa、500hPa的模式高度层的结果。

[0134] 相似于位温的扰动，单变量与多变量集合敏感性拥有类似的形态分布，但单变量集合敏感性相对于多变量集合敏感性存在对敏感性明显高估的现象。这两种方法的负大值区主要集中在台风中心到几百公里随高度而变宽的区域内，即这片区域的增湿会带来预报误差下降，即台风预报强度的增强。这与台风眼墙区域增湿会带来后续强度增强的特征是吻合的。但在扰动850hPa高度上的水汽时，单变量集合敏感性估计的预报响应函数变化的极值为-10.6hPa，多变量集合敏感性的极值为-8.8hPa。而在500hPa高度上，单变量、多变量的极值则分别为-10.3hPa和-2.5hPa，充分展现了单变量集合敏感性高估敏感性的缺点。

[0135] 图7和图8同图5，仅扰动变量替换为在目标时刻24小时前的切向风和径向风。

[0136] 单变量、多变量集合敏感性一致证明，在目标时刻24小时前，台风气旋式环流增强会带来后续台风增强。而850hPa台风东侧、500hPa台风西侧和东南侧的入流增强，高层出流增强会带来台风增强。同样，从数值上看，多变量集合敏感性有效的更正了单变量集合敏感性对敏感性高估的问题，特别是在对流层的高层。

[0137] 图9为扰动在目标时刻48小时前的850hPa高度上的(a)(e)位温(T)、(b)(f)水汽混合比(Q)、(c)(g)切向风(TW)、(d)(h)径向风(RW)，利用(a)-(d)单变量集合敏感性和(e)-(h)多变量集合敏感性得到的预报响应函数的变化(单位hPa)。

[0138] 当引导时间延长至48小时，单变量、多变量集合敏感性分析的主要物理结论与引导时间为24小时相似，但是敏感性的估计值略有减小。多变量集合敏感性的预报响应估计仍然小于单变量集合敏感性。

[0139] 总的来说，单变量和多变量集合敏感性对敏感性的形态分布的估计是相似的，但多变量集合敏感性可以较单变量集合敏感性更加精确地估计敏感性，特别是在台风高层。这可能是由于影响因子的相关尺度(Correlation length scale)随高度上升增大，因此高层存在相关性的影响因子增多，则多变量集合敏感性考虑因子之间相互作用的优势随高度上升而愈发显著。

[0140] 4.2与真实模式相应对比

[0141] 图10为添加扰动到初始场，并积分模式所得的真实的预报响应函数的变化(横坐标)与集合敏感性所估计的预报响应函数的变化(纵坐标)的对比，(a)-(b)单变量集合敏感性、(c)-(d)多变量集合敏感性。第1、2列分别是引导时间为目标时刻前24、48小时的结果。图中虚线为最小二乘法回归线，RMSE为估计值相对于真值的均方根误差。

[0142] 当引导时间为目标时刻前24小时时，多变量集合敏感性最小二乘法回归线更接近对角线，即其与模式真值的比更接近1：1，RMSE也只有不到单变量的一半，说明多变量集合敏感性对预报响应函数的估计要明显由于单变量集合敏感性。当引导时间为目标时刻前48小时时，多变量集合敏感性的回归线虽然出现了少量偏差，但仍然比单变量集合敏感性准确，RMSE也仅有单变量的40％。

[0143] 总而言之，由模式模拟所得的真实的预报响应函数的变化验证表明，多变量集合敏感性对敏感性的估计较传统的单变量集合敏感性更精确和可靠。