CN104484539A - 一种多参数时空气象场模拟方法及装置 - Google Patents

Abstract

本发明提供了一种多参数时空气象场模拟方法及装置，本方法中预先构建预设模型数据库，预设模型数据库是依据多站点多参数的历史气象数据进行主成分分析并对生成的每个主成分序列进行分频获得多个本征模态函数序列、然后对每个本征模态函数序列进行分段形成缀片，并对每个缀片进行模拟后得到的，所以预设模型数据库均能够融合历史气象数据的多个参数，并能够代表历史气象数据的规律。将所需仿真的目标时间段和太阳高度角集输入至预设模型数据库中，预设模型数据库根据自身已有的历史气象数据规律进行仿真，得到仿真的多站点多参数气象矩阵。本发明能够合理评估未来气象场，为环境风险分析学科提供技术支撑，使得评估某种合理的极端情景成为可能。

Description

一种多参数时空气象场模拟方法及装置

技术领域

本发明涉及环境风险分析领域，尤其涉及一种多参数时空气象场模拟方法及装置。

背景技术

环境是一个系统，对于污染物在环境迁移转化的过程而言，气象场是驱动该过程的驱动力。气象场的变化会加速(强化)或减缓(弱化)污染物的环境行为。气象场指的是空间各位置上多种气象参数的时间演变，具有鲜明的多参数系统动力学特征。例如，大气污染物扩散过程中，在不同时间段的气压、风速、风向、降水、太阳辐射等地面气象参数，以及不同时间段的风速、风向、温度、高度等若干高层气象参数都会污染源的扩散产生影响，将以上参数简称时间参数和空间参数，因此在评估污染源经由大气扩散导致的环境风险时，需要依赖气象场的空间特征和时间特征，即多参数时空气象场。

现有在环境风险分析中经常使用的方法为历史气象场截取法：从多参数时空气象场的历史观测记录中抽取若干片段作为环境风险评估所用的情景，假设区域时空气象场的历史记录完全代表了该区域时空气象场的变动性。因此，每个从中截取的片段都“适当”的代表了未来可能出现的气象场片段；即采用历史观测记录作为未来气象场。

有研究采取了一种改进的历史气象场截取法，在总结某特定气象参数，如降水强度历史变化规律基础上，合理估计特定气象参数在未来可能的变化——如均值变化或方差变化或自相关系数变化，在对该特定气象参数进行模拟预测基础上，将该特定气象参数的预测值代替历史气象参数，作为未来情景中的数值，这种方法可称之为单气象参数时空模拟法。

但单气象参数时空模拟法为降水、气温等单参数的时空模拟，而不是多参数时空气象场模拟，因为该方法未考虑多参数时空变化的协同效果，如果简单将单参数模拟结果相互叠加，则可能产生违反物理运动规律的“伪气象场”。

发明内容

本发明提供了一种多参数时空气象场模拟方法及装置，能够依据历史数据合理评估未来气象场，为环境风险分析学科提供技术支撑，使得评估某种合理的极端情景成为可能。

为了实现上述目的，本发明提供了以下技术手段：

一种多参数时空气象场模拟方法，包括：

获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列以及白噪声叠加得到一个主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；

按上述步骤在所述预设模型数据库中获得多个主成分序列的仿真序列；

对所述多个主成分序列的仿真序列进行主成分逆变换得到仿真后的标准气象矩阵；

对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。

优选的，所述在每个仿真缀片集内确认适合的多个仿真缀片包括：

针对每个本征模态函数序列的仿真缀片集中每个仿真缀片，判断是否满足以下三个条件：仿真缀片的主导频率与该本征模态函数序列的主导频率近似相等、仿真缀片的幅值与该本征模态函数序列的幅值近似相等且仿真缀片的相位与前一仿真缀片的尾端相位相容；

若上述三个条件均满足则确认该仿真缀片合适。

优选的，构建所述预设模型数据库的过程包括：

获取目标区域内若干个地面气象站和若干个高空气象站在若干年内不同时刻的多参数气象数据，并构成一个以时刻为行、所述气象数据为列的气象矩阵，对所述气象矩阵进行Z变换获得标准气象矩阵；

对所述标准气象矩阵进行主成分分析获得多个主成分向量和与每个主成分对应的主成分序列，多个主成分向量组成的矩阵为主成分矩阵；

对每个主成分序列执行经验模态分解，获得多个本征模态函数序列；

对每个本征模态函数序列分段获得多个缀片，利用状态子空间线性系统模型对每个缀片进行模拟，获得每个缀片的缀片模型；

将同一个本征模态函数序列的多个缀片模型作为一个缀片模型集；

与每个主成分序列的多个本征模态函数序列一一对应的多个缀片模型集组成预设模型数据库。

优选的，在对每个主成分序列执行经验模态分解之前还包括：

对每个主成分序列利用状态子空间线性系统模型进行模拟，得到每个主成分序列的子空间模型，若某一个或多个主成分序列的子空间模型的模拟结果小于允许误差则表明与该模型对应的主成分序列建模成功；

对未建模成功的每个主成分序列执行经验模态分解。

优选的，还包括：

将所述目标时间段和所述太阳高度角集输入至建模成功的一个或多个子空间模型；

将一个或多个子空间模型的输出结果作为对应的一个或多个主成分序列。

优选的，所述对每个本征模态函数序列分段获得多个缀片，利用状态子空间线性系统模型对每个缀片进行模拟，获得每个缀片的缀片模型包括：

通过频率分析在每个本征模态函数序列中获取一个主导频率，并计算得到主导频率所对应的周期；

以所述每个本征模态函数的周期的预设倍数为缀片长度，以所述周期的预设数量分之一为移动步长，按所述缀片长度和所述移动步长截取每个本征模态函数序列，获得每个本征模态函数序列的多个缀片；

利用状态子空间线性系统模型对每个本征模态函数序列的每个缀片进行模拟，获得每个本征模态函数序列中与每个缀片对应的缀片模型。

优选的，对所述标准气象矩阵进行主成分分析获得多个主成分序列包括：

对所述标准气象矩阵进行主成分分析得到若干个主成分序列；

在若干个主成分序列中选取荷载系数大于1，累计方差解释率达到预设值的前k个主成分序列作为有效主成分序列，其余无效主成分序列视为噪声；

将k个有效主成分序列作为最终的多个主成分序列。

优选的，所述对每个主成分序列执行经验模态分解，获得多个本征模态函数序列包括：

对每个主成分序列执行经验模态分解，获得若干个本征模态函数序列；

在若干个本征模态函数序列按主导频率由高至低选取多个本征模态函数序列，其余主导频率的本征模态函数序列视为噪声。

优选的，对所述多个主成分序列进行主成分逆变换得到仿真后的标准气象矩阵包括：

在多个主成分序列后加入若干维白噪声补足维数后，将补充白噪声后的多个主成分序列与主成分矩阵的逆相乘得到维数一致的标准气象矩阵。

一种多参数时空气象场模拟装置，包括：

获取数据单元，用于获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

获取主成分序列单元，用于将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列以及白噪声叠加得到一个主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；按上述步骤在所述预设模型数据库中获得多个主成分序列的仿真序列；

仿真单元，用于对所述多个主成分序列的仿真序列进行主成分逆变换得到仿真后的标准气象矩阵；对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。

本发明提供了一种多参数时空气象场模拟方法及装置，将所需仿真的目标时间段和太阳高度角集输入至预设模型数据库中，预设模型数据库根据自身已有的历史气象数据规律，对输入的仿真数据进行模拟，得到每个本征模态函数序列的仿真序列，多个本征模态函数序列的仿真序列形成一个主成分序列的仿真序列，多个主成分序列的仿真序列进行主成分逆变换、以及反Z变换即可得到仿真后的气象矩阵。

本方法中预先构建预设模型数据库，预设模型数据库包含多个缀片模型集，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个，所以预设模型数据库能够融合历史气象数据的多个参数，并能够代表历史气象数据的规律。本申请依据预设模型数据进行仿真模拟，所以得到仿真后的气象矩阵能够代表历史气象数据的规律，能够依据历史数据合理评估未来气象场，为环境风险分析学科提供技术支撑，使得评估某种合理的极端情景成为可能。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例公开的一种多参数时空气象场模拟方法的流程图；

图2为本发明实施例公开的一种多参数时空气象场模拟方法中建立预设模型数据库的流程图；

图3为本发明实施例公开的一种多参数时空气象场模拟方法中对每个本征模态函数序列分段模拟的流程图；

图4为本发明实施例公开的一种多参数时空气象场模拟装置结构图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明提供了一种多参数时空气象场模拟方法，包括：

步骤S101：获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

在本发明之前预先构建预设模型数据库，由于预设模型数据库的缀片模型是依据状态子空间线性系统模型构建的，该模型的输入是时间和太阳高度角，因此在需要仿真时，需要准备仿真所需的目标时间段和与目标时间段内每个时刻对应的太阳高度角，将目标时间段和太阳高度角集作为预设模型数据库的输入，经过预设模型数据库计算后可输出仿真的气象场。

步骤S102：将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列叠加得到一个主成分序列的仿真序列，如果模型集中还包括该主成分序列的噪声模型，则应再叠加上该噪声序列得到该主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；

预设模型数据库的构建过程为：依据多参数历史气象数据构成气象数据矩阵，然后对气象数据矩阵进行Z变换获得标准的气象数据矩阵，然后对气象数据矩阵进行主成分分析，获得代表气象数据矩阵的多个主成分序列，对每个主成分序列进行经验模态分解得到多个本征模态函数序列，每个本征模态函数序列只有一个主导频率和主导幅值，对每个本征模态函数序列进行时间分段模拟获得多个缀片模型。

所以预设模型数据库的架构为：多个主成分序列、每个主成分序列下具有多个本征模态函数序列、每个本征模态函数序列对应一个缀片模型集，一个缀片模型集内具有多个缀片模型。

以一个主成分序列为例进行说明，其余主成分序列的实现过程与此类似，将目标时间段和太阳高度角集分别输入至与每个缀片模型集的多个缀片模型，每个缀片模型均对目标时间段和太阳高度角集进行模拟，并输出自身模拟的仿真缀片，从而得到多个仿真缀片，多个仿真缀片组成仿真缀片集。

仿真缀片集中多个仿真缀片有的合适的、有的是不合适的，因此需要在仿真缀片集中进行筛选，筛选条件如下：仿真缀片的主导频率与该本征模态函数序列的主导频率是否近似相等、仿真缀片的幅值与该本征模态函数序列的幅值是否近似相等且仿真缀片的相位与前一仿真缀片的尾端相位是否相容；若上述三个条件均满足则确认该仿真缀片合适。

由于每个缀片模型为对一个本征模态函数序列进行时间分段获得的，对缀片模型仿真的结果需要尽量还原本征模态函数序列，每一个本征模态函数序列有一个主导频率和主导幅值，所以每个仿真缀片的主导频率要与本征模态函数序列的主导频率近似一致，否则说明该仿真缀片的仿真效果较差删除不用，仿真缀片的幅值与本征模态函数序列的幅值近似一致，否则说明该仿真缀片的仿真效果较差删除不用。每个仿真缀片均为一个时间片段，有起始时刻和最后时刻，每个仿真缀片在起始时刻的相位，需要与前一仿真缀片的最后时刻的相位近似相等，这样两者仿真缀片才能连缀起来。

因此只有频率、幅值、相位均达到要求的仿真缀片才是适合的仿真缀片，删除不满足要求的仿真缀片，这三条准则是本征模态函数序列信号特征的忠实反映——频率变化不大、振幅变化不大、连续无突变。

将每个仿真模型集中达到要求的仿真缀片连缀起来，形成一个本征模态函数序列的仿真序列，将多个仿真模型集的达到要求的仿真缀片连缀起来，便形成了多个本征模态函数序列的仿真序列，将一个主成分序列对应的所有本征模态函数序列的仿真序列组合在一起形成一个主成分序列的仿真序列。

步骤S103：按步骤S102在所述预设模型数据库中获得多个主成分序列；

按步骤S102依据预设模型数据库依次形成多个主成分序列，多个主成分序列即为按照目标时间段和太阳高度角集仿真得到的目标气象场的主成分序列。

步骤S104：对所述多个主成分序列进行主成分逆变换得到仿真后的标准气象矩阵；

由于预设模型数据库中的主成分序列是经过主成分分析后得到的，因此为了得到正常的意义上的气象场，将上述步骤S103中多个主成分序列进行主成分逆变换。具体的，首先在多个主成分序列后加入若干维白噪声补足维数后，主成分序列与主成分矩阵的逆相乘得到维数一致的标准气象矩阵。

步骤S105：对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。

再对步骤S104中标准气象矩阵进行反Z变换生成真正的仿真后的气象数据矩阵。

本发明提供了一种多参数时空气象场模拟方法，本方法中预先构建预设模型数据库，预设模型数据库包含多个缀片模型集，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个，所以预设模型数据库能够融合历史气象数据的多个参数，并能够代表历史气象数据的规律。本申请依据预设模型数据进行仿真模拟，所以得到仿真后的气象矩阵能够代表历史气象数据的规律，能够依据历史数据合理评估未来气象场，为环境风险分析学科提供技术支撑，使得评估某种合理的极端情景成为可能。

下面详细介绍构建预设模型数据库的过程，如图2所示，包括以下步骤：

步骤S201：获取目标区域内若干个地面气象站和若干个高空气象站在若干年内不同时刻的多参数气象数据，并构成一个以时刻为行、所述气象数据为列的气象矩阵，对所述气象矩阵进行Z变换获得标准气象矩阵；

以目标区域作为构建预设模型数据库的对象，后续输入至预设构建模型的目标时间段和太阳高度角集也应该为该目标区域内的，否则仿真结果不可信。按照时间序列，将目标区域内地面气象站(设共n个)记录的地面气象参数和区域内及区域周边高空气象站(设m个高空气象站)的探空气象参数并列，则可构成一个以时间为行、各地面及高空气象站参数为列的气象数据矩阵。

气象数据矩阵囊括了地面和高空的多个气象参数，因此，气象数据矩阵的各列不仅隐含了气象站的位置信息，而且隐含了各个气象参数之间、各个气象站不同气象参数之间的相关关系；气象数据矩阵的各行则进一步包含了这些不同位置不同气象参数之间关系的时间变化规律。换言之，气象数据矩阵包含了气象场的所有信息。

对气象数据矩阵每个参数进行Z转化，使之各列均为均值为0、方差为1的数列，该Z转化后的矩阵为标准的气象数据矩阵。

步骤S202：对所述标准气象矩阵进行主成分分析获得多个主成分向量和与每个主成分对应的主成分序列，多个主成分向量组成的矩阵为主成分矩阵；

由于标准气象矩阵列数过多，信息冗余，为了简化后续处理过程，在标准气象矩阵中主成分分析，提取有信息成分的主成分序列，具体执行过程如下：对所述标准气象矩阵进行主成分分析得到若干个主成分序列；在若干个主成分序列中选取荷载系数大于1，累计方差解释率达到预设值的前k个主成分序列作为有效主成分序列，其余主成分序列在本实施例中暂且视为噪声；将k个有效主成分序列作为最终的多个主成分序列。

荷载系数大于1说明该主成分序列为有效信息成分，荷载系数小于1说明该主成分序列的信息量不超过随机噪声，可舍弃，选取荷载系数大于1的多个主成分序列，并且多个主成分序列的累计方差解释率能够达到预设值，预设值为预先设定的值，如70％，即选取的主成分序列在所有主成分序列中所占的贡献率大于70％。

步骤S203：对每个主成分序列执行经验模态分解，获得多个本征模态函数序列；

由于天气受到的影响因素众多，所以主成分序列的频率特征较为复杂，通过简单的线性系统模型难以模拟，因此对每个主成分序列进行经验模态分解技术，将每个主成分序列分为若干个本征模态函数序列；对每个主成分序列执行经验模态分解，获得若干个本征模态函数序列。

本征模态函数序列IMF并非严格意义上的单频率成分，而是频率在相对缓慢变化之中的信号，这是气象场特征非平稳性的必然反映，不仅IMF的频率不固定，而且其振幅也是变化的。各IMF振幅的相对大小代表了其IMF成分的能量相对大小；某IMF信号随时间的振幅变化则反映了该IMF成分能量的时变性。尽管如此，某特定IMF的振幅和频率整体不会发生剧烈突变。

步骤S204：对每个本征模态函数序列分段获得多个缀片，利用状态子空间线性系统模型对每个缀片进行模拟，获得每个缀片的缀片模型；

尽管每个本征模态函数序列非平稳，但在某个较短的时间内，其特征是相对稳定的，是可以用较简单的系统模型进行模拟近似的。所以对本征模态函数序列分段得到多个缀片，对每个缀片利用线性子空间模型进行模拟，得到每个缀片的模型。

下面详细介绍步骤S204的具体执行过程，如图3所示，包括以下步骤：

步骤S301：通过频率分析在每个本征模态函数序列中获取一个主导频率，并计算得到主导频率所对应的周期；

步骤S302：以所述每个本征模态函数序列的周期的预设倍数为缀片长度，与所述周期的预设数量分之一为移动步长，按所述缀片长度和所述移动步长截取每个本征模态函数序列，获得每个本征模态函数序列的多个缀片；

步骤S303：利用状态子空间线性系统模型对每个本征模态函数序列的每个缀片进行模拟，获得每个本征模态函数序列中与每个缀片对应的缀片模型。

每个本征模态函数序列均有一个主导频率，通过频率计算得到该主导频率的主导周期；以主导周期的预设倍数作为该本征模态函数序列的长度，以主导周期的若干分之一作为移动步长，以长度和移动步长截取本征模态函数序列，从而获得多个长度一致的缀片，每个缀片的时间段较短，因此在一个缀片时间段内可以近似为线性系统，可以利用线性系统模型——状态子空间模型进行模拟。

步骤S205：将同一个本征模态函数序列的多个缀片模型作为一个缀片模型集，与每个主成分序列的多个本征模态函数序列一一对应的多个缀片模型集组成预设模型数据库。

同一个本征模态函数序列的多个缀片的缀片模型组成一个缀片模型集，多个本征模态函数序列形成一个主成分序列模型，多个主成分序列模型形成气象场模型，由此预设模型数据库构建完成，后续便可使用预设模型数据库进行后续的仿真过程。

图2所示的构建预设模型数据库的过程，在步骤S202与步骤S203之间还可以包括：

对每个主成分序列利用状态子空间线性系统模型进行模拟，得到每个主成分序列的子空间模型，若某一个或多个主成分序列的子空间模型的模拟结果小于允许误差则表明与该模型对应的主成分序列建模成功；对未建模成功的每个主成分序列执行步骤S203。

首先利用线性系统模型对所有主成分序列进行模拟，若某一个主成分序列的模拟结果与真实结果之间的误差小于允许误差，则说明该主成分序列仅以季节变化而变化，因此能够采用线性系统来表达，所以该主成分序列建模成功。

一般情况下，至少有一个主成分序列是仅与季节相关的，所以通过该步骤找出与季节相关的主成分序列，该主成分序列便可采用线性系统模型构建模型，由于线性系统模型非常简单，可以省去繁琐的计算过程，因此可以减少后续的大量计算过程。

在仿真时，将所述目标时间段和所述太阳高度角集输入至建模成功的一个或多个子空间模型；将一个或多个子空间模型的输出结果作为对应的一个或多个主成分序列。将该步骤中得到的主成分序列与步骤S203中得到的主成分序列一并作为仿真所需的主成分序列，再执行后续的主成分序列逆操作和反Z变换，从而得到气象数据矩阵。

下面介绍本发明的具体实施例：

本发明技术以区域内时空气象场的历史记录为数据基础。通常一个区域(如方圆数百公里)内都会有若干地面气象站，若干高空气象站，地面气象站会有各地面气象参数的时值记录或日值记录，高空气象站则按照规范在格林威治时间0点和12点(即北京时间8点和20点)一日两次记录探空气象参数记录。

(一)构建气象数据矩阵A

按照时间顺序，将区域内地面气象站(设共n个)记录的地面气象参数和区域内及区域周边高空气象站(设m个高空气象站)的探空气象参数并列，则可构成一个以时间为行、各地面及高空气象站参数为列的气象数据矩阵A。当然，由于地面气象站与高空气象站的记录频率不同，二者应分别通过内插函数插值至统一时间点上的气象参数值，如：如线性内插或状态子空间线性系统模拟内插等常用方法，插值至北京时间每个整点时刻，或每日的气象参数值。

气象数据矩阵A本身代表了气象场。因为各列代表的气象参数本身隐含了气象站的空间位置信息，各行则显在的表达了时间序列。这种思想是气象学中常用的经验正交函数EOF方法的理论基础。只是，常规EOF往往仅针对一个气象参数如降水，而此处A则囊括了地面和高空的多个气象参数，因此，A的各列不仅隐含了气象站的位置信息，而且隐含了各个气象参数之间、各个气象站不同气象参数之间的相关关系；A的各行则进一步包含了这些不同位置不同气象参数之间关系的时间变化规律。换言之，A包含了气象场的所有信息。

(二)标准化气象数据矩阵A

气象数据矩阵A内个参数首先案列进行Z转化，使之各列均为均值为0、方差为1的数列，该Z转化后的矩阵为标准化的气象数据矩阵B。

(三)主成分分析并析出有信息主成分序列

为得到气象场的时空分布规律，对标准化的气象数据矩阵B进行主成分分析(PCA)，以荷载系数为标准，选取荷载系数大于1且累计方差解释率达到70％的前k个主成分序列作为有信息主成分序列，而将剩余的主成分序列视作噪声。根据具体工作需要，累计方差解释率可以适当调整，如50％，90％等。根据主成分分析的基本理论，主成分向量代表了气象场的空间格局，多位置多参数的空间格局，而与主成分向量对应的主成分序列指代该空间格局的时间演变规律，称之BP_q,q＝1、2……k。

(四)尝试系统模拟每个有信息主成分序列BP_q

第q个主成分序列BP_q实际上是单变量时间序列，BP_q可能受到季节(日期)和日内时间(时刻)的影响，也可能与季节或时刻因素没有显在关系而是受制于地形、地表覆盖等地面因素与季节和时刻因素的综合作用。事实上，k个有信息主成分序列中总是会有至少一个是与季节性显著相关的；如果气象数据矩阵A或B的时间分辨率是12h或更短，则信息主成分序列中总是会有至少一个是与日内时刻显著相关的。

季节性可以用太阳高度角作为自变量，自变量定义域(0,90]来表示(太阳高度角的计算方法可参见气象数据规范)，时刻则可直接作为自变量，定义域[0,24)。采用状态子空间(State Subspace)线性系统模型对BP_q进行系统模拟，如果模拟结果令人满意，如决定系数r²>0.7，则该主成分序列BP_q建模成功，保留该模型M_q。

并将所需要仿真的输入信号(季节、时刻)输入M_q，产生BP_q的仿真结果后续步骤详见步骤(七)；若模拟结果不令人满意，则意味着BP_q并非由季节性或时间主导，则进入下一步分频模拟步骤。

(五)对有信息主成分序列BP_q进行分频并缀片化模拟

主要包括两个步骤，首先就是对BP_q频率分解为若干频率成分；其次是对每个频率成分进行模拟，而采取的模拟策略是缀片化，即分段模拟。

(1)对BP_q进行分频处理

第q个主成分序列BP_q不能被季节参数和时间参数模拟，直接原因是频率成分较为复杂而不是简单地以年为周期或以日为周期。理论上可以采用傅里叶变换实现分频，但傅里叶变换实际上只有在频率成分较简单时才能以较少的频率成分近似原信号。而气象信号本身的非平稳性特征使得常规傅里叶分频技术不再适合。

本发明中推荐使用希尔伯特-黄分频技术，即气候学和气象学中已经有应用的经验模态分解(EMD)技术。应用该技术，可以将BP_q分解为若干本征模态函数序列(IMF)之和；实际上EMD可以只关注有限个能量较大的以及具有用户感兴趣的频率特征的IMF(如式(1)中的w个IMF)，而将其余IMF之和整体视作一个噪音IMF，以r_w表示，如式(1)中的r_w：

${\mathrm{BP}}_q=\underset{1}{\overset{w+1}{\mathrm{\Σ}}}{\mathrm{IMF}}_v=\underset{1}{\overset{w}{\mathrm{\Σ}}}{\mathrm{IMF}}_v+r_w\·\·\·\·\·\·\left(1\right)$

所得到的有限个IMF中的任何一个IMF_v，都并非严格意义上的单频率成分，而是频率在相对缓慢变化之中的信号。这也正是气象场特征非平稳性的必然反映。不仅IMF的频率不固定，而且其振幅也是变化的。各IMF振幅的相对大小代表了其IMF成分的能量相对大小；某IMF信号随时间的振幅变化则反映了该IMF成分能量的时变性。尽管如此，某特定IMF的振幅和频率整体不会发生剧烈突变。

(2)对每个IMF_v进行缀片化模拟

尽管IMF_v非平稳，但在某个较短的时间内，其特征是相对稳定的，是可以用较简单的系统模型进行模拟近似的，第v个本征模态函数序列IMF_v的缀片化模拟方法，其核心思想是对IMF_v时间序列分段缀片，对每个缀片进行模拟。

(a)确定缀片长度

第v个本征模态函数序列IMF_v是时间序列信号，为方便在此将IMF_v记作原始信号S0。S0虽然不是单一频率信号，但一般存在一个主导信号，即能量最大的频率成分。通过频率分析得到该最大频率成分，设其周期长为P。建议以P的若干倍(如4倍)为缀片长度L；以P的若干分之一(如十分之一)为移动步长(缀片初始相位)t。按照L和t截取S0，生成一系列缀片信号集{S1_i}。

(b)以状态子空间模型模拟每个缀片信号

对于第i个缀片信号S1_i，以线性系统模型——状态子空间模型进行模拟，获得所需输入信号的太阳高度角、时刻以及采样顺序本身。首先优选确定子空间维度级数，从1级子空间模型到某级(一般可设为最高4级)子空间模型分别模拟S1_i，以模拟值与信号值的中值离差和平均离差为准则，优选出子空间模型级数。

然后以该级数的子空间模型对S1_i在进行优化模拟，记录模拟的误差指标FIT(平均相对误差)和子空间模型参数，作为缀片信号S1_i的模拟模型。上述过程对{S1_i}中每个缀片信号重复使用，获得一个本征模态函数序列的一系列备选模型{M_i}。

(六)仿真产生有信息主成分序列BP_q

这一步是步骤(五)相对应，也包括两个步骤，首先就是对BP_q的各频率成分进行仿真；其次复合每个频率成分产生BP_q的仿真结果。仿真意味着已经根据实际需要设定好了输入信号——太阳高度角集(季节信号)、时间段(日内时刻信号)、采样顺序(自然设定)。

(1)对每个IMF_v产生仿真信号

每个S0信号都对应一系列备选模型{M_i}。用每个M_i仿真输入信号序列，产生仿真信号缀片。不过，某一个仿真信号缀片是否采用，取决于下述评价：

仿真信号缀片的主导频率与原信号(S0)的主导频率近似相等；且

仿真信号缀片的振幅与原信号(S0)的振幅近似相等；且

仿真信号缀片的相位与前一仿真信号缀片的尾端相容。

只有以上三条准则(频率，振幅，相位)均达到可接受的程度，才采纳该仿真信号缀片；否则删除该备选模型。这三条准则也是IMF信号特征的忠实反映——频率变化不大、振幅变化不大、连续无突变。将符合条件的仿真信号缀片连缀起来组成一个本征模态函数序列。

值得指出，在M_i仿真时应允许M_i的状态初始值优化调整，以便适应仿真输入信号与原缀片信号S1_i的差异。在这个意义上，M_i只是决定了状态子空间的系统动力学特性，而不是直接固定初始状态。

将上述仿真过程按照设定的数目依次重复，缀片前后相连，则最终生成了设定输入长度上IMF_v的仿真信号

(2)对IMF_v仿真信号复合成BP_q的仿真结果

此处复合过程为：由式(1)可知，只需将IMF_v，v＝1、2……w的仿真信号相加即可，当然还应加上白噪声r_w，如式(2)所示。

(七)将所有仿真序列主成分逆变换

这一步是步骤(三)的反向操作，即将所有仿真序列q＝1、2…k与主成分矩阵的逆相乘得到仿真的标准化气象数据矩阵

值得提醒的是由于只有1..k个，因此应以噪声向量(白噪声)补足维数后再与主成分矩阵的逆相乘，否则矩阵维数不对应。噪声向量应按照步骤(三)中噪声方差大小设定，以补齐仿真数据方差。此处噪声方差系指累计方差解释率的方差之外的方差，当累计方差解释率取70％时，则噪声方差为总方差的30％。

(八)将仿真的标准化气象数据矩阵反Z变换

这一步是步骤(二)的反向操作。即将所有仿真序列乘以各列向量的标准偏差再以对应向量的平均进行平移，得到具有各气象参数单位和数值的仿真气象数据矩阵而中的各列即为各地面和高空气象站各气象参数的时间序列，当然，这些多维时间序列保持了气象场的时空格局。

到此，以所表示的仿真多参数时空气象场作为环境风险分析所需的气象场情景即可进行污染物迁移扩散，进而定量化分析环境风险的时空分布特征。

本技术方案中，基本上可以理解为“分析→分解→模拟→仿真→仿真结果反分解(复合)”的过程，前一半是获得模型，后一半是用模型仿真并对应的反向操作最终形成仿真的多参数时空气象场。

如图4所示，本发明还提供了一种多参数时空气象场模拟装置，包括：

获取数据单元31，用于获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

获取主成分序列单元32，用于将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列以及白噪声叠加得到一个主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；按上述步骤在所述预设模型数据库中获得多个主成分序列的仿真序列；

仿真单元33，用于对所述多个主成分序列的仿真序列进行主成分逆变换得到仿真后的标准气象矩阵；对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。

本发明提供了一种多参数时空气象场模拟装置，本装置中预先构建预设模型数据库，预设模型数据库包含多个缀片模型集，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个，所以预设模型数据库能够融合历史气象数据的多个参数，并能够代表历史气象数据的规律。本申请依据预设模型数据进行仿真模拟，所以得到仿真后的气象矩阵能够代表历史气象数据的规律，能够依据历史数据合理评估未来气象场，为环境风险分析学科提供技术支撑，使得评估某种合理的极端情景成为可能。

下面以华北地区的详细实验过程为例对本发明进行验证：

应用华北地区UTM平面投影600km×600km区域上的多参数时空气象场模拟对本技术的应用做一说明。

该区域根据中国气象数据共享平台网提供的历史记录数据，选定了在2001-2010年共10年间的9个地面气象站(标准站)和范围覆盖本区域的共14个探空气象站的历史记录数据用作气象场原始数据。选择这些站以及数据年限的原因在于这些站在此期间气象数据相对丰富，记录缺失较少。地面气象参数为3h时值(0，3,6，……21)数据，包括9个参数：气压、风向、风速、气温、相对湿度、总云量、低云量、云底高度、降水强度；探空数据为格林尼治时间0点和12点数据，本应用案例中采用了925hPa，850hPa，700hPa，500hPa，300hPa共5个等压面，每层包括4个参数：高程、温度、风向、风速。则本区域气象场包括9*9+14*5*4＝81+280＝361参数，其中地面参数81个，高空参数280个。

应用本技术进行了建模。在主成分序列分解后，按照70％的方差解释率选择10个主成分序列(10PC)进行了模拟和仿真。每个主成分序列在EMD分解过程中按照高频到低频最多分解21个IMF。作为可能的建模简化以及本技术应用潜力，在本案例中也用60％方差解释率和50％方差解释率分别选择前5个主成分序列(5PC)和前3个主成分序列(3PC)进行了建模，并进行了比较(见下文)。

本技术允许两种环境风险分析用气象场情景生成方法：①对某种感兴趣的气象特征进行整体缩放；②对某种感兴趣的气象特征的某频段成分进行缩放。无论哪种情景生成方法，都会对原气象场参数间的交叉相关关系进行一定的扭曲。但是，如果这种扭曲没有在统计上达到显著的水平，则这种生成的情景仍然是可信的。在本应用案例中，通过对PC的识别，以第2主成分序列(PC2)为例，分别整体上扩大5％，10％，25％，或以其高频成分(日内变化)第1-3IMF扩大5％，10％，25％分别产生情景，与不予调整的默认仿真情景以其共产生了7中情景。PC2表达了扣除主导季风性影响后的参与气象场的西风带特征。

整体上，上述7种情景与采用5PC和3PC(无缩放调整)情景，共产生了9种情景，每个情景模拟出10年数据，形成9个数据集，如表1所示。

表1 气象场模拟情景

本技术的模拟能力主要以以下三个方面的能力进行评价：

(1)361个参数的每个参数的数值分布是否与原历史数据的数值分布相符；

(2)361个参数相互之间的相关关系(交叉相关关系)与原历史数据集参数之间的交叉相关系相符；

(3)361个参数每个参数的自相关性与原历史数据集中对应参数的自相关性相符。

以上三个方面只有全部满足，才能证实本技术的确实实现了多参数时空气象场的模拟能力，即保留了数值分布、空间结构、时间动态的特征。

表2和表3分别显示了9个地面气象站和14个探空气象站中通过前述检验的数量。

表2 模拟的9个地面气象站各参数通过检验的个数及模拟值均值与记录值均值之差

表3 模拟的14个探空气象站各参数通过检验的个数及模拟值均值与记录值均值之差

可见，地面气象参数有48％，高空气象参数有33％通过统计检验。为进一步探究模拟数据集的可用性，特别是针对那些未通过检验的参数，各参数模拟值与记录值的平均差别也列出，见表2和表3的“均值差”一列。从均值差可见，这些差别都是微不足道的。比如气压均值差为74Pa,风向均值差为-16°，湿度均值差为2.5个百分点等。云底高度均值差为-185m，考虑到历史记录数据中该参数的巨大误差，这个均值差也在可容忍范围内。高空气象参数中高程相差1m左右，气温相差1℃左右，风向相差1°做有关，风速相差1°左右，这些差别与实际记录误差相比，微不足道。

综上各参数百分位点对比图、均值差、统计检验结果，可以得出结论，就数值分布而言，无论地面还是高空参数，模拟结果与历史记录结果相符。个别统计上显著不同的参数，相较于历史记录数据误差(该误差在记录数据中不能反映出来)都是微不足道的。

本实施例方法所述的功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算设备可读取存储介质中。基于这样的理解，本发明实施例对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该软件产品存储在一个存储介质中，包括若干指令用以使得一台计算设备(可以是个人计算机，服务器，移动计算设备或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (10)

1.一种多参数时空气象场模拟方法，其特征在于，包括：

获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列以及白噪声叠加得到一个主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；

按上述步骤在所述预设模型数据库中获得多个主成分序列的仿真序列；

对所述多个主成分序列的仿真序列进行主成分逆变换得到仿真后的标准气象矩阵；

对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。

2.如权利要求1所述的方法，其特征在于，所述在每个仿真缀片集内确认适合的多个仿真缀片包括：

针对每个本征模态函数序列的仿真缀片集中每个仿真缀片，判断是否满足以下三个条件：仿真缀片的主导频率与该本征模态函数序列的主导频率近似相等、仿真缀片的幅值与该本征模态函数序列的幅值近似相等且仿真缀片的相位与前一仿真缀片的尾端相位相容；

若上述三个条件均满足则确认该仿真缀片合适。

3.如权利要求1所述的方法，其特征在于，构建所述预设模型数据库的过程包括：

获取目标区域内若干个地面气象站和若干个高空气象站在若干年内不同时刻的多参数气象数据，并构成一个以时刻为行、所述气象数据为列的气象矩阵，对所述气象矩阵进行Z变换获得标准气象矩阵；

对所述标准气象矩阵进行主成分分析获得多个主成分向量和与每个主成分对应的主成分序列，多个主成分向量组成的矩阵为主成分矩阵；

对每个主成分序列执行经验模态分解，获得多个本征模态函数序列；

对每个本征模态函数序列分段获得多个缀片，利用状态子空间线性系统模型对每个缀片进行模拟，获得每个缀片的缀片模型；

将同一个本征模态函数序列的多个缀片模型作为一个缀片模型集；

与每个主成分序列的多个本征模态函数序列一一对应的多个缀片模型集组成预设模型数据库。

4.如权利要求3所述的方法，其特征在于，在对每个主成分序列执行经验模态分解之前还包括：

对每个主成分序列利用状态子空间线性系统模型进行模拟，得到每个主成分序列的子空间模型，若某一个或多个主成分序列的子空间模型的模拟结果小于允许误差则表明与该模型对应的主成分序列建模成功；

对未建模成功的每个主成分序列执行经验模态分解。

5.如权利要求4所述的方法，其特征在于，还包括：

将所述目标时间段和所述太阳高度角集输入至建模成功的一个或多个子空间模型；

将一个或多个子空间模型的输出结果作为对应的一个或多个主成分序列。

6.如权利要求3所述的方法，其特征在于，所述对每个本征模态函数序列分段获得多个缀片，利用状态子空间线性系统模型对每个缀片进行模拟，获得每个缀片的缀片模型包括：

通过频率分析在每个本征模态函数序列中获取一个主导频率，并计算得到主导频率所对应的周期；

以所述每个本征模态函数序列的周期的预设倍数为缀片长度，以所述周期的预设数量分之一为移动步长，按所述缀片长度和所述移动步长截取每个本征模态函数序列，获得每个本征模态函数序列的多个缀片；

利用状态子空间线性系统模型对每个本征模态函数序列的每个缀片进行模拟，获得每个本征模态函数序列中与每个缀片对应的缀片模型。

7.如权利要求3所述的方法，其特征在于，对所述标准气象矩阵进行主成分分析获得多个主成分序列包括：

对所述标准气象矩阵进行主成分分析得到若干个主成分序列；

在若干个主成分序列中选取荷载系数大于1，累计方差解释率达到预设值的前k个主成分序列作为有效主成分序列，其余无效主成分序列视为噪声；

将k个有效主成分序列作为最终的多个主成分序列。

8.如权利要求3所述的方法，其特征在于，所述对每个主成分序列执行经验模态分解，获得多个本征模态函数序列包括：

对每个主成分序列执行经验模态分解，获得若干个本征模态函数序列；

在若干个本征模态函数序列按主导频率由高至低选取多个本征模态函数序列，其余主导频率的本征模态函数序列视为噪声。

9.如权利要求3所述的方法，其特征在于，对所述多个主成分序列进行主成分逆变换得到仿真后的标准气象矩阵包括：

在多个主成分序列后加入若干维白噪声补足维数后，将补充白噪声后的多个主成分序列与主成分矩阵的逆相乘得到维数一致的标准气象矩阵。

10.一种多参数时空气象场模拟装置，其特征在于，包括：

获取数据单元，用于获取目标时间段和太阳高度角集，所述太阳高度角集中每个太阳高度角与所述目标时间段内每个时刻一一对应；

获取主成分序列单元，用于将所述目标时间段和所述太阳高度角集、输入至预设模型数据库的每个本征模态函数序列的缀片模型集内，获得每个本征模态函数序列的仿真缀片集，在每个仿真缀片集内确认适合的多个仿真缀片，将适合的多个仿真缀片连缀得到一个本征模态函数序列的仿真序列，将一个主成分序列所含的所有本征模态函数序列的仿真序列以及白噪声叠加得到一个主成分序列的仿真序列；所述预设模型数据库包含多个缀片模型集，每个缀片模型集内具有多个缀片模型，且每个缀片模型集与一个本征模态函数序列对应，每个本征模态函数序列为对多参数历史气象数据进行主成分分析获得每个主成分序列后、对每个主成分序列进行经验模态分解得到的多个本征模态函数序列中的一个；按上述步骤在所述预设模型数据库中获得多个主成分序列的仿真序列；

仿真单元，用于对所述多个主成分序列的仿真序列进行主成分逆变换得到仿真后的标准气象矩阵；对所述标准气象矩阵进行反Z变换获得仿真的气象数据矩阵。