一种基于经验模态分解的自动气象站数据质量控制方法

技术领域

本发明涉及自动气象站采集数据的质量控制领域，特别发明了一种对自动 气象站采集的气温数据的质量控制方法。

背景技术

目前国内的自动气象站具有站点分布密集、地形差异大、测站环境恶劣、 数据采集和传输自动化程度高、资料实时性强、中小尺度天气现象明显等特点, 其观测数据的质量问题比人工常规地面站观测资料显得更复杂和严重，因而必须 对自动气象站的观测数据进行质量控制。我国现采集数据的质量控制在两个层面 进行,一是在观测台站进行基本质量控制；第二是在数据处理中心进行,数据处 理中心除采用综合时间一致性、内部一致性检查外,还应对仪器传感器失效、长 期漂移进行评估。根据我国气象业务体制,除基层台站进行基本质量控制外,省 级应对自动站实时数据进行全面系统质量控制。但是，自动气象站在进行数据采 集过程中容易受到大量随机噪声、系统噪声、微气象噪声和粗大噪声的干扰；这 些噪声是目前质量控制方法不能消除的，从而影响自动气象站观测数据的质量， 不能很好的服务于中尺度、短时效的精细天气预报。

传统基于经验模态分解的去噪方法是将IMF分量的高频部分直接去掉再重 构去噪，这种方法忽视了去掉高频部分存在的真实信息，在去掉噪声的同时也去 掉了有用的信息；同时，IMF分量间的高频是相对的局部高频，在低频的分量中 仍然存在相对的高频和噪声，这是传统基于模态分解去噪方法不足的地方。

发明内容

本发明所要解决的技术问题是：

提供一种基于经验模态分解的自动气象站数据质量控制方法，解决了目前质 量控制方法不能消除随机噪声、系统噪声、微气象噪声和粗大噪声的问题，提高 了自动气象站观测数据的质量。

本发明为解决上述技术问题，采用如下技术方案：

一种基于经验模态分解的自动气象站数据质量控制方法，包括以下步骤：

步骤1.采集自动气象站温度观测数据x(t)，其中t为采样时间，且t为大于0 的自然数；

步骤2.对x(t)进行基本质量控制，得到基本控制后的数据

步骤3.根据公式对基本控制后的数据进行经验模态分 解，得到本征模分量IMF_i与的趋势项r；其中k为大于1的自然数，i∈k；

步骤4.计算本征模分量IMF_i的能量E_i，根据E_i绘制能量柱状图，并观察能 量柱状图，得到本征模分量能量的阶跃点M，根据该阶跃点确定存在噪声的本 征模分量IMF_j；其中M为自然数，且1<M≤k；j为自然数，且1≤j≤M-1；

步骤5.对基本控制后的数据与本征模分量IMF_j进行时间序列相关性分 析得到相关系数v_j，根据相关系数v_j计算赋权本征模分量NEWIMF_j和赋权后的 趋势项NEWr；

步骤6.根据赋权本征模分量NEWIMF_j和赋权后的趋势项NEWr进行数据重 构得到去噪声后的数据，完成去噪质量控制。

步骤2中所述基本质量控制包括界限值检查、气候极值检查、内部一致性检 查、时变一致性检查、空间一致性检查、特殊天气事件检查。

步骤4中所述计算本征模分量的能量E_i运用如下公式计算：

$E_i=\underset{i=1}{\overset{k}{\Sigma }}{\mathrm{IMF}}_i{\left(t\right)}^2,$

其中IMF_i(t)为第i个本征模分量的时间序列，t为采样时间点。

步骤5中所述计算赋权本征模分量运用如下公式计算：

NEWIMF_j=v_j·IMF_j；j=1,2…M-1

所述赋权后的趋势项运用如下公式计算：

$\mathrm{NEWr}=\underset{p=M}{\overset{k}{\Sigma }}{\mathrm{IMF}}_p+r.$

v_j运用如下公式计算：

$v_j=\frac{\Sigma {\mathrm{IMF}}_j·\tilde{x}\left(t\right)-\frac{\Sigma {\mathrm{IMF}}_j\Sigma \tilde{x}\left(t\right)}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)}}{\sqrt{(\Sigma {\tilde{x}}^2\left(t\right)-\frac{{\left(\Sigma \tilde{x}\left(t\right)\right)}^2}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)})(\Sigma {\mathrm{IMF}}_j^2-\frac{{\left(\Sigma {\mathrm{IMF}}_j\right)}^2}{\mathrm{size}\left({\mathrm{IMF}}_j\right)})}},$

其中size(IMF_j)分别表示与IMF_j的时间序列长度。

所述步骤6中数据重构运用如下公式计算：其中为去噪声后的数据。

与现有技术相比，本发明具有如下有益效果：

通过消除观测数据的随机噪声和高斯噪声来实现对自动气象站观测数据的 去噪质量控制，弥补了传统质量控制方法的不足，能够提高自动气象站观测数据 在长期气候研究和短时预报的可靠性。

附图说明

图1为本发明方法的流程图。

图2为本发明基本质量控制后的数据本征模分量及趋势项波形图。

图3为本发明基本质量控制后的数据本征模分量能量变化图。

图4为本发明赋权本征模分量及赋权后的趋势项波形图。

图5为本发明去噪前后数据的希尔伯特谱图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明：

为本发明方法的流程，如图1所示，对单站数据进行采集，然后对采集的数 据做基本质量控制后，再进行EEMD（经验模态）分解得到本征模分量和趋势项， 进一步对本征模分量和趋势项做加权处理，得到赋权本征模分量和赋权后的趋势 项，最后对赋权后的本征模分量和赋权后的趋势项做数据重构，得到去噪后的数 据，去噪过程完成。

以下将对南京编号为58238站点2007年1月份每隔10分钟观测一次的温 度数据进行实施例分析，进一步说明本发明：

1、采集自动气象站温度观测数据x(t)，t为采样时间，且t为大于0的自然 数，本实施例选择t=1，2…4464。

2、对x(t)经过界限值检查、气候极值检查、内部一致性检查、内部一致性 检查、时变检查、空间一致性检查等基本质量控制，得到基本控制后的数据

3、根据公式对基本控制后的数据进行经验模态分解， 得到本征模分量IMF_i与的趋势项r；其中k为大于1的自然数；

4、利用能量计算公式求取本征模分量的能量，绘制能量柱 状图，并观察能量柱状图，得出E_i的能量阶跃点M，根据该分界点确定存在噪 声的本征模分量IMF_j；其中IMF_i(t)为第i个本征模分量的时间序列，t为采样时 间点，M为自然数，且1<M≤k；j为自然数，且1≤j≤M-1；本征模分量的 能量阶跃点M，即本征模分量的能量变化突变的点。

5、对基本控制后的数据与本征模分量IMF_j进行时间序列相关性分析得 到相关系数v_j，v_j运用如下公式计算：

$v_j=\frac{\Sigma {\mathrm{IMF}}_j·\tilde{x}\left(t\right)-\frac{\Sigma {\mathrm{IMF}}_j\Sigma \tilde{x}\left(t\right)}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)}}{\sqrt{(\Sigma {\tilde{x}}^2\left(t\right)-\frac{{\left(\Sigma \tilde{x}\left(t\right)\right)}^2}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)})(\Sigma {\mathrm{IMF}}_j^2-\frac{{\left(\Sigma {\mathrm{IMF}}_j\right)}^2}{\mathrm{size}\left({\mathrm{IMF}}_j\right)})}},$其中size(IMF_j)分别表示与IMF_j的时间序列长度；根据相关系数v_j计算赋权本征 模分量NEWIMF_j和赋权后的趋势项NEWr，计算赋权本征模分量运用如下公式 计算：NEWIMF_j=v_j·IMF_j(j=1,2…M-1)；赋权后的趋势项运用如下公式计算：

$\mathrm{NEWr}=\underset{p=M}{\overset{k}{\Sigma }}{\mathrm{IMF}}_p+r.$

6、根据赋权本征模分量NEWIMF_j和赋权后的趋势项NEWr进行数据重构得 到去噪声后的数据，数据重构运用如下公式计算： $\tilde{x}{\left(t\right)}_{\mathrm{new}}=\underset{j=1}{\overset{M-1}{\Sigma }}{\mathrm{NEWIMF}}_j+\mathrm{NEWr},$其中为去噪声后的数据，完成去噪质量 控制。

本实施例将分解为5个本征模分量IMF₁、IMF₂、IMF₃、IMF₄、IMF₅和 1个趋势项r，其相对频率关系为：IMF₁＞IMF₂＞IMF₃＞IMF₄＞IMF₅＞r；如 图2所示，横坐标是采样数据点间隔，纵坐标是本征模分量对应的数值，为方便 阅读波形图，每个量依次累计减10处理。

利用公式计算五个本征模分量IMF₁、IMF₂、IMF₃、IMF₄、 IMF₅和1个趋势项r的能量，其中E₁～E₅分别为五个本征模分量IMF₁、IMF₂、 IMF₃、IMF₄、IMF₅的能量；E_r为1个趋势项r的能量；根据E₁～E₅、E_r绘制 能量变化柱状图，如图3所示，根据柱状图显示，在E₅处发生能量跳跃，因此 得到能量阶跃点M=5，判定IMF₁、IMF₂、IMF₃、IMF₄为存在噪声的分量。

根据公式$v_j=\frac{\Sigma {\mathrm{IMF}}_j·\tilde{x}\left(t\right)-\frac{\Sigma {\mathrm{IMF}}_j\Sigma \tilde{x}\left(t\right)}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)}}{\sqrt{(\Sigma {\tilde{x}}^2\left(t\right)-\frac{{\left(\Sigma \tilde{x}\left(t\right)\right)}^2}{\mathrm{size}\left(\tilde{x}\left(t\right)\right)})(\Sigma {\mathrm{IMF}}_j^2-\frac{{\left(\Sigma {\mathrm{IMF}}_j\right)}^2}{\mathrm{size}\left({\mathrm{IMF}}_j\right)})}}$分别计算与IMF₁、IMF₂、IMF₃、IMF₄的相关系数v_j,(j=1,2,3,4)，v₁表示与IMF₁的 相关系数，v₁=0.8057；v₂表示与IMF₂的相关系数，v₂=0.1068；v₃表示与IMF₃的相关系数，v₃=0.2049；v₄表示与IMF₄的相关系数，v₄=0.2024。

利用相关系数v₁～v₄和本征模分量IMF₁～IMF₄根据公式 NEWIMF_j=v_j·IMF_j,j=1,2,3,4计算赋权本征模分量NEWIMF₁～NEWIMF₄；利 用IMF₅、趋势项r根据公式NEWr=IMF₅+r计算赋权后的趋势项NEWr，如图 4所示。

利用公式对NEWIMF_j和NEWr进行数据重构， 为南京编号为58238站点2007年1月份每隔10分钟观测的温度数据经 过本系统去噪后的数据，整个去噪过程完成。

为分析本发明的有益效果，此处用美国NOAA国家气候数据中心的同时段再 分析气温数据作为分析数据，本发明方法数据分别与小波去噪、传统EMD去噪进 行信噪比和均方根误差进行对比，对比结果表格如下表所示：

本方法的信噪比明显高于小波去噪法、传统EMD去噪法；均方根远小于小波 去噪法、传统EMD去噪法。

气温数据去噪前后的希尔伯特谱图如图5所示，去噪前图上噪点密集，去噪 后波形平滑、干净，主要信息保持良好。

通过上表及图5对比分析，本发明方法去噪效果明显。