基于温度场的气象时空数据动态可视化

徐晓文，聂芸

（华北计算技术研究所北京100083）

随着时代的发展，计算机的技术越发的先进，同时数据观测和采集仪器的进步，可以生成大量的气象数据。气象数据具有很强的时空性，表现在一个气象信息处理系统中的数据源既有同时间不同空间的数据，也有同空间不同时间的数据[1]。同时，通过不同尺度的观测也会有不同的比例和精度。

气象可视化科学研究已经比较成熟，能够对气象观测的要素数据、分析预报的产品进行按空间高度层次、按时次的分析和显示，如等值线分析显示、填色分析与显示等，在气象水文预报业务和气象水文信息应用系统中广泛使用[2]。随着气象水文信息对日常生活、行业工作的影响日益被重视，在应用系统中提供气象水文环境信息的感知能力成为迫切的需求，同时也为了更好地进行气象环境应用分析，更加充分、直观地展现气象水文信息，对气象信息可视化提出了更高的要求。以往单层次、单时次的静止数据展现不能充分体现气象数据时空连续的特征；天气图等面向专业人员的方式也难以被一般用户所接受。

本论文研究基于已有的气象可视化技术基础，在不改变可视化的具体实现的过程情况下，利用现有的可视化数据，通过多种方法实现中间插值，从而获得气象时空上连续的可视化效果。

1 系统总体设计

1.1 系统实现整体步骤

在实现气象时空数据动态可视化过程中，需要每一个中间图像生成的数据，通过这些图像的连续播放，从而实现气象时空数据的动态可视化。在很多情况下，气象数据都是单时次，不连续的，所以需要将这些不连续的单时次的数据组合起来，如果直接生成可视化图像进行连续播放，那么大部分的气象数据在感官上都存在明显的撕裂感。所以需要利用这些原本的单时次，不连续的数据生成一连串的数据来补充这个画面减少撕裂感的产生。这时候就需要应用到插值技术。

由于是做气象时空数据动态可视化，所以需要进行插值做出中间动画来实现动态展示。考虑到分段低次插值函数都有一致收敛性，但光滑性较差，对于气象时空数据动态可视化要求有一定的光滑度，即最好保证存在二阶连续导数。所以本文采用最常用的三次样条函数[3]。三次样条插值（Cubic Spline Interpolation）简称Spline插值，我们选择的是三次样条插值[3]。这种插值技术比较容易生成平滑的曲线，并且在多个关键点上数据相关，在保证关键的准确的获取到的数据的精度上，可以合理的生成出相关的数据[4]。同时由于数据是矩阵的形式，所以利用近似最邻近特征匹配方法进行采样，获取出一连串的点，从而构建出合理的三次样条函数[5]。

之后通过线性条插值补全数据。通过这种方法形成的新的矩阵，不仅在数据上连贯性良好，还在初始气象数据矩阵十分精确。

1.2 系统实现流程设计

图1表示了实验的具体流程。根据流程图示意，首先数据格式要求为矩阵形式。那么初始数据为给定的一连串的时间上的矩阵数据。

假定初始存在5个数据场，数据矩阵分别为M1，…，M5。这是一个固定时间间隔的序列矩阵，分别依照时间顺序编号1，2，…，5。如图2所示。

图1 系统流程图

图2 初始矩阵序列图

根据图2所示，存在一个矩阵序列，包含依照时间排序的5个数据矩阵（其中2和4没有画出来），矩阵大小全部都为m×n，它们依照时间顺序编号为1，2，3，4，5，本文主要实现的就是在这一个连续的数据矩阵序列中插入中间数据并利用这些数据和原始的数据生成一个动态的可视化效果。根据生成动画的原理，为了获得流畅的动画效果，需要在这些数据矩阵插入中间数据。

2 系统实现步骤

2.1 系统初始数据的获取

文中研究中假定初始数据为矩阵型网格数据。这种类型的数据是实际业务中最常见的格式。实际使用中，通过气象水文资料的数据处理，就能够直接得到这类数据，通常按照行列顺序存储。对于气象实况等基于观测站点的离散数据，可以通过网格化处理形成矩阵型数据。

气象数据的网格节点数（行列数）一般是按照在一定地理范围内固定间隔分布来计算的，比如全球预报数据，按照4度、1度划分，可形成90*45，360*180的数据矩阵[6]。

本章以4度划分的矩阵为例进行算法说明。

同时，为了获取到比较合理的数据，本文在实验中所采用的数据为NOAA（美国国家海洋和大气管理局）的数据（https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html），下 载 Air Temperature数据并提取文件里的数据，缩小矩阵规模进行实验。

2.2 提取矩阵的关键点

文中存在5个初始数据矩阵M1，…，M5，这是一个依照时间顺序排序的序列矩阵，分别赋予编号1，2，3，4，5。现在本文需要在这些矩阵中建立一定的映射关系用来方便插值。取矩阵M1和M2作为样本，都是同等大小的m×n矩阵，要求矩阵M1的M1（x，y）{（0＜x＜90，0＜y＜45，且x，y为整数，其中x，y对应的是矩阵中的坐标位置）}求对应到矩阵M2的最近临近相似点[7]。

考虑到数据矩阵空间上的对应与时间上的连续，所以需要考虑空间对应关系以及时间上的变化关系[8]。

本文所采用的方法就是求取

这个由M2矩阵的一部分元素组成的新的矩阵与M1(x,y)最相近的，此处最相近不仅仅代表空间上的相近，也有时间上数值的相近。由于采用的是近似最临近特征匹配，所以本文需要考虑两个矩阵中的距离因素，即为空间要素，假设这个距离为D12。为了确定距离的具体数值，本文做了以下的定义：

在M2中求取M1(x,y)中的近似最临近特征，依照本文设定距离D12的数值为M1(x,y)与M2(x,y)对应，那么距离D12的数值为0，然后与M2(x-1,y)，M2(x,y-1)，M2(x+1,y)与M2(x,y+1)对应，那么距离D12的数值为1，其他距离D12的数值为2。

之后应该计算所有条件下对每个对应匹配占有的影响数值I，所利用的公式I=（M2（x，y）-M1（x，y））+D12。其中M1(x,y)可以分别变换成

矩阵中任何的元素。之后选取I最小的点作为匹配。实现过程如下所示[9]：

1）首先将原始数据矩阵M1与M2进行分割，通过多区域采样来提高数据的精确度。如图3所示，两个矩阵按照红线所示分成等分的9个部分，分别是取横轴和纵轴的三分之一作为一个分划位置。这些分给出来的矩阵依照先从左到右再从上到下的规则 分 别 命 名 为M1p1，M1p2，…，M1p9与M2p1，M2p2，…，M2p9。之后各个相应的区域找对应关系，比如M1p1与M2p1对应；

图3 原始数据矩阵区域划分示意图

2）以M1p1与M2p1对应为例，本文采用的对应关系如图4所示。为圆点到星点的匹配过程。分割后M1p1与M2p1均为m/3×n/3大小的数据矩阵，其中圆点的取值为这个区域的中心点M1p1（m/6，n/6），星点的位置对应矩阵如下：

之后的对应关系重新依照这个方法对应。同理M2与M3之间的关系重新匹配。每个关键点取由前一个中心点确定的作为样本数据。即第一个样本点的数值是Nm1p1，根据这个关系找到的是Nm2p1，但是Nm2p1可能已经不是矩阵M2p1(m/6,n/6)的数据了。但是寻找Nm3p1的时候仍然以矩阵M2p1的M2p1(m/6,n/6)元素开始参考运算。Nm2p1并不会参与到M2p1与M3p1的运算中；

图4 数据关键点的查找示意图

3）后运算影响数值数据I，当存在多个相同的I时，并且刚好I的取值为最小，那么根据气象时空数据时间和空间都具有变换的关系，采用匹配点的顺时针方向顺序选取第一个I值最小的点作为匹配点。

通过上述方法，可以在数据矩阵M1和M2之间建立出一个对应的关系。本文开始分割矩阵M1与M2，得到了多个数据矩阵，分别为M1p1，M1p2，…，M1p9与M2p1，M2p2，…，M2p9。之后以M1p1与M2p1为例，通过M1p1的中心点找到对应的矩阵M2p1的所需求的数据和位置，之后再利用矩阵M2p1的中心点找到矩阵M3p1对应的数据和位置，之后就可以得到一连串的5个数值，分别为Nm1p1，Nm2p1，… ，Nm5p1。N就是全部 5个序列矩阵演变关系的一个基准向量，建立了5个矩阵的对应关系。同理可以获得其他的分块矩阵的基准向量。本文以之前的矩阵获取到的基准向量为例。

2.3 获取三次样条函数

根据2.2可知，通过近似最临近特征匹配本文在初始数据上获取到了一连串的5个数值，分别为Nm1p1，Nm2p1，… ，Nm5p1。那么接下来为了构建出这些关键点的中间插值数据，本文考虑到气象数据的流体性质，采用三次样条插值[10]。首先根据初始数据是一个存在时间顺序关系的数据矩阵，那么根据这个信息，赋予一个变量T。那么假设T=[1，2，3，4，5]。同理存在随着这个序数改变的变量，即为上一小节获取的 5 个数值Nm=[Nm1p1，Nm2p1，…，Nm5p1]，那么就可以根据这个条件得到三次样条函数。主要步骤如下所示[11]：

1）三次样条插值多项式Sn(x)在每个小区间[xi-1,xi]上是3次多项式，其在此区间上的表达式如下：

其中，hi=xi+1-xi。因此，只要确定了Mi的值，就确定了整个表达式；

2）在此处x的取值为T，S（x）的取值为Nm。令：

则Mi满足如下n-1个方程：

本文采用第一种边界方法，所以得到以下的方程：

3）根据上文信息，存在n+1个方程，有着n+1个Mi未知数，那么根据得到的方程解出Mi，就可以得到方程的解。

2.4 得到插值数据

根据2.3得到的三次样条插值的方程，假设本文需要求取M1p1和M2p1之间的插值数据。假设两个矩阵之间插值数据为20个。这些数据和M1p1和M2p1是等时间间隔的，所以假设序号也是等间隔的。那么假设这些数据为N1l2p1n1，N1l2p1n2，…，N1l2p1n20。那根据时间次序的序列的划分T=[1，1.05，1.10，…，2]。根据上一小节的方程，可以得到M1p1和M2p1之间的三次样条函数，其中自变量为T，通过T的数值来确定对应的插值点的数值。

如图5所示，在两个矩阵之间插入数据，不过现在插入的只是一些关键点的数据，所以插入矩阵的数据不完全，它插入的是星号的数据，所以需要对剩余部分进行填充。

图5 中间数据的插入过程

以N1l2p1n1为例子，假设已经取到了N1l2p1n1，现在利用这个数据进行补全成插值矩阵M1l2p1n1，这个插值矩阵对应的序号为1.05。

图6 补全插入数据的过程

2.5 进行动态可视化

如上所示，M1l2n1已经获取到了，同理可以得到更多的插值数据。本文采用的是现有的可视化手段，插值矩阵和原始数据构造相同，所以可以直接进行可视化展示[14]。

3 系统成果与分析

3.1 系统成果展现

本文采用最常用的等值线方法进行可视化。可视化效果如图7所示。

为了方便观察，图8为图7通过提取的左上角部分区域的可视化显示图片，可以观测到等值线有变化过程，但是大致的区域的变动不大。说明这个数据插值进行可视化是有效果的。同时在运行过程中并没有感受到明显的卡顿，而且动画的过度也比较流畅，证明了优化的效果。

3.2 系统分析

本文研究需要大量的成熟数值方法、图形图像处理方法支撑，为了能够快速实验多种方法，使用了科学计算工具matlab作为实验平台。

实验中采取5个连续45×90的初始数据矩阵，每个数据之间插入300个数据，去除可视化步骤。在实验中内存占用如下：

图7 可视化显示序列图

图8 可视化显示部分序列图

物理内存：4135 MB。

交换页面：5036 MB。

虚拟内存：6341 MB。

在使用e3处理器，16 g内存的环境下，所需要的运行时间为15.051630 s；即为每个数据插入平均耗时不到0.02 s，满足可视化要求。

4 结论

本论文研究基于已有的气象可视化技术基础，在不改变可视化的具体实现的过程情况下，利用现有的可视化数据，通过多种方法实现中间插值，利用插值的数据进行可视化展现，然后连续播放这些可视化效果，从而获得气象时空上连续的可视化效果。论文设计的方法，分离了可视化方法和数据处理方法，使得动态展现的效果能够针对多种可视化方法实施，便于扩展应用；数据处理方法简单，可并行执行，便于应用实现。论文有些工作还有待进一步展开[15]，包括：

1）可视化效果可扩展。考虑到与现有系统的结合使用，本文在说明气象时空数据可视化方法时，仅选取了等值线场和色彩场方法。对于时空数据可视化的呈现方式本身，也是值得研究的问题。

2）方法验证和改进。论文采用的实验数据类型有限，还需要采用不同类型不同规模的真实数据，对方法进行测试，更好地确定方法的适应性、有效性。