基于时空插值的高分辨温盐数据生成方法及装置

1.本发明涉及海洋环境建模技术领域，更具体的说是涉及基于时空插值的高分辨温盐数据生成方法及装置。


背景技术：

2.我国的海洋装备对海洋环境的耐久性和适应性问题十分突出，在装备试验环节，真实、复杂的海洋环境模型是提高装备试验效果、有效评估装备性能的重要支撑。因此，如何构建真实、复杂的海洋环境模型，为海洋装备虚拟试验提供丰富、真实的海洋环境数据，对海洋装备试验具有重要意义。
3.目前，在海洋环境模型的构建上，只依据了空间维度的观测数据，而没有考虑时间维度上的观测数据；同时，观测数据在时空上分布离散、稀疏的特点，成为基于时空数据构建海洋环境模型的又一瓶颈。
4.因此，如何突破瓶颈，同时考虑时间维度和空间维度的观测数据，构建海洋环境模型是本领域技术人员亟需解决的问题。


技术实现要素：

5.有鉴于此，本发明提供了一种基于时空插值的高分辨温盐数据生成方法及装置，针对部分观测数据时空分布离散、稀疏的情况，研究基于时空插值的海洋环境建模方法，实现对观测数据的重构，生成符合所需时空分布规律的、高分辨的网格化海洋环境数据。
6.为了实现上述目的，本发明采用的技术方案包括：
7.第一、一种基于时空插值的高分辨温盐数据生成方法，包括，
8.获取温盐的时空监测数据，并对所述时空监测数据的总体分布形式进行平稳性分析，提取满足平稳性的温盐时空数据；
9.根据所述温盐时空数据中的时间数据得到时间变异函数，根据所述时空数据中的空间数据得到空间变异函数；
10.根据所述时间变异函数和所述空间变异函数，构建时空变异函数；
11.利用时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据，所述时空kriging插值法的表达式为：
[0012][0013]
其中，(ti,si)表示样本点i的时空域坐标，λi(i＝1,2,...,n)表示待求的权重系数。
[0014]
第二、一种基于时空插值的高分辨温盐数据生成装置，包括，
[0015]
温盐时空数据获取单元，用于获取温盐时空监测数据，并提取满足平稳性要求的温盐时空数据；
[0016]
时间和空间变异函数拟合单元，用于根据温盐时空数据中的时间数据拟合得到时
间变异函数，根据温盐时空数据中的空间数据拟合得到空间变异函数；
[0017]
时空变异函数构建单元，用于根据所述时间变异函数和空间变异函数，构建时空变异函数；
[0018]
高分辨率温盐数据生成单元，用于通过时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据。
[0019]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于时空插值的高分辨温盐数据生成方法及装置，通过获取温盐时空监测数据，并进行平稳性分析，从而获得满足平稳性要求的时空数据，采用适合的模型与方法拟合得到时间变异函数以及空间变异函数后，进一步构建时空变异函数，有效地利用时空邻近样本进行插值。
[0020]
本发明中利用时空插值方法对时空变异函数进行插值，可有效克服现有技术中缺陷，即同时考虑时间维度和空间维度的观测数据，并基于时空分布离散、稀疏的观测数据生成标准化的、具有高分辨率的温度和盐度网格数据；以便于实现对海洋环境模型中的海水温度、盐度参数的生成，以及依据生成的温盐数据采用工程模型实现海水密度、声速参数的生成。
[0021]
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。
[0022]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0023]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0024]
图1为本发明基于时空插值的高分辨温盐数据生成方法流程图；
[0025]
图2为本发明时间变异函数拟合过程示意图；
[0026]
图3为本发明时空kriging插值流程图；
[0027]
图4为本发明argo数据筛选前后的空间分布，其中(a)为筛选前的数据空间分布，(b)为筛选后的数据空间分布；
[0028]
图5为本发明插值前后数据对比图，(a)为插值前后的盐度数据对比图，(b)为插值前后的温度数据对比图；
[0029]
图6为本发明趋势面拟合前的温度数据直方图分布情况；
[0030]
图7为本发明趋势面拟合后的温度数据直方图分布情况；
[0031]
图8为本发明所选时空范围20m深度处的温度数据趋势面拟合结果图；
[0032]
图9为本发明中2016年5月预留观测数据点的空间分布及插值误差图；
[0033]
图10为本发明各插值点的误差图。
具体实施方式
[0034]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完
整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0035]
本发明实施例公开了一种基于时空插值的高分辨温盐数据生成方法，其中，时空插值法主要通过构建时空变异函数，从而有效利用时空邻近样本进行插值，以得到标准化的、具有高分辨率的温度和盐度网格数据；进而，实现对海洋环境模型中的海水温度、盐度参数的生成，以及依据生成的温盐数据采用工程模型实现海水密度、声速参数的生成。
[0036]
实施例一
[0037]
本实施例中，基于时空插值的高分辨温盐数据生成方法包括，
[0038]
获取温盐的时空监测数据，并对所述时空监测数据的总体分布形式进行平稳性分析，提取满足平稳性的温盐时空数据；
[0039]
根据所述温盐时空数据中的时间数据得到时间变异函数，根据所述时空数据中的空间数据得到空间变异函数，
[0040]
根据所述时间变异函数和所述空间变异函数，构建时空变异函数；
[0041]
利用改进的时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据，所述改进的时空kriging插值法的表达式为：
[0042][0043]
其中，(ti,si)表示样本点i的时空域坐标，λi(i＝1,2,...,n)表示待求的权重系数。
[0044]
kriging插值法是一种在地理统计学中常用到的插值方法。主要基于区域化变量的理论知识，通过构建空间变异函数的模型以及对采样数据空间上的自相关分析，按照自然现象的空间变化规律得到无偏最优估计量，从而进行数据插值。
[0045]
本实施例中，时空kriging插值法是在kriging插值法的基础上，通过对空间域到时空域的扩展实现了空间变量向时空变量的扩展，将原本的空间变异函数拓展为了时空变异函数，并综合利用邻近数据进行插值。以此既改善了原有的kriging法指定时刻下空间采样点不足所导致的插值精度差的情况，同时也充分考虑了采样数据的时空范围内的变化规律。
[0046]
因此，本发明中时空插值法以采样点的时空自相关分析为基础，依据观测数据的时空分布和变化规律，综合考虑每个采样点间的时空关系来进行插值，适用于具有离散的时空分布特征的数据集，即本技术中时空kriging插值法，可有效解决海洋观测数据在时空分布上稀疏散布的问题，并能实现海洋观测数据的高分辨的标准化的网格构建。
[0047]
实施例二
[0048]
下面结合具体实施例对本发明基于时空插值的高分辨温盐数据生成方法中各步骤进行说明。
[0049]
第一步、获取温盐的时空监测数据，
[0050]
本发明中时空监测数据类似于空间变量的分析，即用dr表示空间域度，r表示空间维度(如二维空间)，t表示时间域度，那么当存在点s＝(x,y)∈dr,t∈t，z(s,t)所表示的就是时空随机变量。
[0051]
本发明中，时空kriging插值法的插值精度主要取决于时空变异函数的拟合效果，而拟合效果则受纯时间变异函数和纯空间变异函数的拟合结果的直接影响。当用于插值的时空数据具有稀疏散布的特征时，就会导致两类变异函数难以拟合。
[0052]
具体的，对于时间监测数据，一般情况下，稀疏散布的时空数据通常由移动设备进行观测采样得到的，如argo浮标数据，使得获得的时空数据通常不具有相对固定的空间坐标，即空间分布具有随机性，也就无法获取到固定空间位置下的时间观测序列；相比于移动采样设备，固定观测站采样的数据则具有一个稳定的空间位置，同时数据的采样频率也保持一致(如tao锚式浮标数据)，能够容易地获取到稳定的时间观测序列，计算拟合时间变异函数，但观测站数量往往较少。
[0053]
因此，本发明提出借助时空范围内的多个观测站获取的观测数据，来计算估计时间变异函数的实验值，并用于时间变异函数的拟合，对于空间坐标处的观测站，其时间变异函数的实验值计算公式如下：
[0054][0055]
由于固定观测站的数量较少，同时插值点处的时间变异情况可能同k个观测站的时间变异情况存在差异，对于插值点处的时间变异函数的拟合，可以先分别拟合出k个观测站所处的空间位置下的时间变异函数γ
t
(sj,h
t
)(j＝1,2,
…
,k)，之后对各时间变异函数的块金系数c0、拱高c和变程a通过反距离加权估计的方式(各观测站相对插值点处的距离进行反距离加权)得到插值点处的实际的时间变异函数γ
t
(h
t
)，具有以下表达式：
[0056][0057]
λj为各观测站的距离权重系数，系数大小由采用的反距离加权法和到插值点的距离共同决定，相应的拟合过程如图2所示。
[0058]
对于空间监测数据，实际中，采样点数是有限的，本技术为了使拟合的空间变异函数能较好地反映空间变化特征，要求在变程内的各个空间距离下的采样点对的数量不少于20。
[0059]
但由于移动采样设备的观测时刻和观测频率不具有一致性，无法做到同步采样，使同一时刻下的采样点数过少，此时拟合得到的空间变异函数不具备代表性，无法准确反映时空变量的空间变化特征。而将所有时段的时空数据叠置进行空间变异函数的拟合求解，又会无法排除时间变异过大所带来的时间变异函数的影响。
[0060]
对此，本技术提出通过将多个时间切片的采样点数进行叠置，计算空间变异函数的实验值，并进行空间变异函数的拟合，从而有效避免采样点数过少的同时，将时间变异所带来的影响减少至最小，其中，多时段叠置的空间变异函数的实验值计算公式为：
[0061][0062]
由式可以看出，当m＝1时该式即为单个时段下的空间变异函数实验值计算公式，m的取值为时间切片数量时则是全时段下的空间变异函数实验值计算公式。因此，m的引入实现了原有的单一时间切片下的空间变异函数的实验值在时间维上向两侧的适度延展，同时
延展的大小取决于m的取值(通常不大于3)。对于延展后得到的各个时段的时间切片分别计算空间变异函数的实验值后，将所有子时段的计算结果叠置并进行空间变异函数的拟合，得到空间变异函数γs(hs)。这样在较短的时段内，通过多个时段的时空数据的叠置保证了空间变异函数拟合所需的数据量的同时，极大地减小了时间变异函数对空间变异函数拟合效果的干扰。
[0063]
第二步、对所述时空监测数据的总体分布形式进行平稳性分析，提取满足平稳性的温盐时空数据；
[0064]
时空kriging插值法能够进行使用的前提是时空数据具有平稳性。这是因为时空kriging插值法在计算过程中(求解变异函数的过程)，要求区域化变量必须存在若干对z(x)和z(x+h)样本的实现，而实际课题中的区域化变量的实际取值是唯一且不重复的，无法满足这一要求。此时，就要求区域化变量具有平稳性，即满足平稳假设。
[0065]
平稳性分析是时空数据分析前必须进行的一步，本发明中采用二阶平稳假设进行检验，
[0066]
另一种实施例中，也可以采用平稳假设、二阶平稳假设或本征假设进行时空数据的平稳性分析。
[0067]
其中，二阶平稳假设对于研究区域内的区域化变量z(x)具有如下约束条件：
[0068]
(a)z(x)在所处的研究区域内具有恒定的数学期望，即：
[0069]
e[z(x)]＝m
[0070]
(b)z(x)在所处的研究区域内存在协方差函数，且协方差函数平稳，数值只与位移量h有关，与当前所处的位置x无关，即：
[0071]
cov[z(x),z(x+h)]＝c(h)
[0072]
依据二阶平稳假设的相关约束条件，就可以进行对时空变量的平稳假设的检验，即对样本的总体分布形式进行检验。
[0073]
本实施例中，采用统计直方图法和趋势面分析法对样本的总体分布形式进行分析。
[0074]
其中，统计直方图法是先将变量xi所对应的观测值按大小进行排序的同时，按照一定的间距对其进行分组。之后统计各个分组下的变量xi的观测值出现的次数，以此绘制出直方图并观察样本的频率分布。若样本呈正态分布的特性，变量xi就具有区域化变量的随机性特征，满足二阶平稳假设。
[0075]
趋势面分析法则是通过构造一个几何曲面，模拟数据在空间或时间上的变化趋势。趋势面分析法采用回归分析法，以实际观测值来推算物理趋势面，使得残差平方和达到最小，即满足如下形式：
[0076][0077]
其中，zi(xi，yi)表示观测值，表示趋势值。
[0078]
通常计算趋势面的方法有多项式函数和傅里叶级数，其中多项式函数计算过程简单，易于操作，常用于趋势面分析法中。趋势面分析法通过对直方图不满足正态分布的数据构造曲面拟合后，求解变量观测值的残差值，再采用统计直方图法进行平稳检验，此时若残
差值满足平稳假设要求，即可将其视为区域化变量，进行时空kriging插值。这种通过趋势面拟合后进行时空kriging插值的方法也可称为残差kriging插值法。
[0079]
第三步、根据所述温盐时空数据中的时间数据得到时间变异函数，根据所述时空数据中的空间数据得到空间变异函数，
[0080]
本发明中，由于时空kriging插值法由于向时空域进行了扩展，无法直接采用曲线拟合进行得到，因此根据观测数据计算得到的实验变异函数值，选择合适的理论模型(即相应的曲线模型)以及拟合方法，分别对时间变异函数和空间变异函数进行曲线拟合。具体包括：
[0081]
(a)变异函数的理论模型选择
[0082]
变异函数的理论模型反映的是变异函数的曲线的整体变化规律，根据曲线最终是否趋于恒定则可以分为有基台值模型(曲线最终趋于一个恒定值)和无基台值模型。有基台值模型通常包含纯块金效应模型、高斯模型、球状模型等，无基台值模型则包含线性无基台值模型。实际中，需要根据时空数据的特点，选择相应的理论模型进行拟合。
[0083]
(b)拟合方法选择
[0084]
在地统计学中，变异函数的实质就是进行曲线拟合，拟合的方法可以采用最小二乘法或者加权回归法。最小二乘法在实现上较为简单，但最终的拟合效果可能会相对较差，而加权回归法则考虑了各个用于的拟合点上所包含的点对数量，以此作为权重进行加权拟合，对于点对数量多的拟合点，其结果更准确，统计特征更好，具有的权重也更大，此时的拟合结果也更优，缺点则是实现上较为复杂。
[0085]
第四步、根据所述时间变异函数和所述空间变异函数，构建时空变异函数；
[0086]
基于第三步获得的空间变异函数和时间变异函数，采用分离型或非分离型的模型进行时空变异函数的构建。其中，其中积合式模型是最常用到的非分离型模型，构建方法如下：
[0087]
γ(hs,h
t
)＝(k1c
t
(0)+k2)γs(hs)+(k1cs(0)+k3)γ
t
(h
t
)-k1γs(hs)γ
t
(h
t
)
[0088][0089]
式中，γs(hs,h
t
)是所构建的时空变异函数，c
st
(0,0)、cs(0)、c
t
(0)分别是时空变异函数、空间变异函数和时间变异函数的基台值。
[0090]
第五步、利用改进的时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据，
[0091]
时空kriging插值在kriging插值的基础上向时空维进行了扩展，相应的时空kriging插值公式也对空间kriging插值公式进行了时空域的扩展；
[0092]
具体的，改进的时空kriging插值法的表达式为：
[0093][0094]
其中，(ti,si)表示样本点i的时空域坐标，λi(i＝1,2,...,n)表示待求的权重系数。
[0095]
对于权重系数λi的求解，需要满足z
*
(s,t)为z(s,t)的无偏估计量并具有最小的估计方差。在区域化变量满足二阶平稳的情况下，λi可以通过时空kriging方程组计算求解，该方程组表达式如下：
[0096][0097]
式中，μ为估计方差，γ((si,ti),(sj,tj))为时空变异函数，其中，(s,t)为时空域中样本点坐标的展开形式。
[0098]
上式构成了一个n元方程组，其中是需要满足的条件函数。该方程组本质是求解下式所表示的条件函数：
[0099][0100]
(s0,t0)表示某一样本点在时空域中的坐标。
[0101]
对该式引入拉格朗日系数进行推导，可以得到如下的矩阵形式；
[0102][0103]
将第四步构建的时空变异函数，代入到上述矩阵中，便可把时空kriging方程组求解转化为矩阵运算问题，通过矩阵求逆和矩阵乘积运算便可求解得到相关的权重系数，这样极大地方便了时空kriging插值法的实现。
[0104]
实施例三
[0105]
本发明进一步公开了一种基于时空插值的高分辨温盐数据生成装置，包括，
[0106]
温盐时空数据获取单元，用于获取温盐时空监测数据，并提取满足平稳性要求的温盐时空数据；
[0107]
时间和空间变异函数拟合单元，用于根据温盐时空数据中的时间数据拟合得到时间变异函数，根据温盐时空数据中的空间数据拟合得到空间变异函数；
[0108]
时空变异函数构建单元，用于根据所述时间变异函数和空间变异函数，构建时空变异函数；
[0109]
高分辨率温盐数据生成单元，用于通过时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据。
[0110]
可选的，还包括数据叠置拟合单元，用于对多时段、多空间的监测数据进行拟合。
[0111]
本技术中上述装置，能够在装备试验过程中，依据试验的需求生成指定时空位置下的海水温度、盐度数据，依据工程模型又可以计算得到对应的海水密度、声速信息，为装备试验，尤其是水下航行设备运行过程的试验，提供需要的参数信息。
[0112]
进一步，还可开发成基于时空插值的海洋环境建模软件，软件能根据使用者的需求生成特定时空范围、指定分辨率的网格化海洋环境数据。
[0113]
实施例四
[0114]
为使本领域技术人员清楚本发明方案的实施过程，本技术进一步提供一种利用实施例二中方法进行温盐数据差值的实施例。
[0115]
本实施例，以argo浮标温盐数据作为主数据，用于实现空间变异函数的拟合并参与到最终的数据插值过程，以tao浮标温盐数据则作为辅助数据，用于实现时间变异函数的拟合。
[0116]
对于时空范围的选择，这里选定了140
°
e～160
°
e,10
°
n～30
°
n的经纬度范围的海区作为数据插值的空间范围，时间范围则是2016年5月。所使用的观测数据以及最后插值生成的温盐数据都处于这个时空范围内，同时最终生成的温盐数据的时空分辨率为1天
×1°×1°
。
[0117]
本实施例中插值过程如图3所示，包括数据预处理、时空变异函数建模求解以及时空kriging插值三个部分；
[0118]
1)数据预处理
[0119]
argo数据为netcdf格式的海洋数据，本实施例采用通过计算机仿真软件先进行读取和预处理后，按日期转存为相应的txt文件，再通过c++编程读取txt文件并进行后续的插值建模。
[0120]
argo数据读取方式可以通过计算机仿真软件配备的netcdf工具箱实现所需数据的读取，包含数据的经度、纬度、温度、盐度、压力(或深度)。读取完成后，由于读取的argo数据包含了整个太平洋区域内的温盐数据，空间范围较选定海区范围大，因此需要对读取数据做相应的空间范围筛选。如图4所示，分别展示的是2016年5月5日的argo数据空间分布情况(a)和筛选后的数据分布情况(b)。
[0121]
由于argo数据在垂直方向上，即深度方向上是离散分布的，不具有统一的深度层。不同空间位置、不同时间下采样得到的温盐数据，采样的深度层也不同。因此，需要将剖面数据的观测值插值到统一的标准层，即确立一个统一的深度层。对于浅层的海水，由于会受到光照和风浪的影响，数据的变化幅度较大；而对于深层的海水，由于所处的海洋环境相对稳定，数据的变化也趋于恒定。此时，为了减少argo数据在深度层上的数据量的同时，不丢失数据深度层方向的变化趋势的细节，可以采用不等间距深度层设置的方式，并通过对原始的观测数据在深度层方向进行数据插值，得到统一标准深度层的温盐数据。
[0122]
深度层方向的数据插值可以采用akima插值法，该方法可以将离散的数据点拟合成光滑曲线，且具有相比三次样条插值拟合曲线更光滑、自然的优势。
[0123]
对于标准深度层的设置，本技术，共设置了57个深度层，深度层间隔随着深度的增加逐渐递增。标准层设置如表1所示，
[0124]
表1
[0125]
序号深度/m序号深度/m序号深度/m序号深度/m序号深度/m15131202530037650491250210141302632038700501300320151402734039750511400
4301615028360408005215005401716029380418505316006501817030400429005417007601918031420439505518008702020032440441000561900980212203346045105057197510902224034500461100111002326035550471150121102428036600481200
[0126]
进一步，图5分别展示了2016年5月5日单个观测点插值前后的盐度(a)、温度数据对比图(b)。
[0127]
可以看出，插值后生成的标准层数据仍能很好地反映原始温度、盐度数据的变化情况，确保了在减小数据量的同时仍能为装备试验提供指定空间位置、指定时刻下真实的温度、盐度的深度变化情况。
[0128]
akima插值到标准层后，需要做二阶平稳假设检验。区域化变量的二阶平稳假设检验是kriging插值前的必要工作。对于海水温度特性，其在空间分布上会存在明显的分布趋势，即海水温度由低纬度向高纬度呈降低趋势。因此，在kriging插值前，需要对数据进行曲面拟合，将温盐数据分解成趋势量和剩余量，再用剩余量进行插值，将插值结果与趋势量相加得到最终结果。
[0129]
对于二阶平稳假设检验，本课题采用统计直方图法进行检验，即通过观察数据的直方图分布情况确定。二阶平稳数据的直方图一般呈正态分布。温盐数据由于其空间分布的趋势影响，其分布可能不呈正态分布。而经过趋势面拟合后，剩余量的分布满足正态分布规律，符合二阶平稳假设要求。如图6、图7所示则是进行趋势面拟合前后的所设定的时空范围内的20m深度层下的温度数据的直方图分布情况。可以看出，在趋势面拟合前的温度数据的直方图难以观察到正态分布的规律，而进行趋势面拟合后的温度数据的直方图则具有明显的正态分布规律，也证明了趋势面拟合后的温度数据满足二阶平稳假设，可以采用时空kriging插值法进行数据插值。
[0130]
趋势面拟合过程中，采用三次多项式曲面拟合的效果要比二次多项式好。因此，本课题采用的是三次多项式拟合，并基于最小二乘法获得最佳拟合曲面。其表达式如下：
[0131]
z＝a0+a1x+a2y+a3x2+a4xy+a5y2+a6x3+a7x2y+a8xy2+a9y3[0132]
其中，分别是拟合的趋势量，拟合数据的经度，拟合数据的纬度。
[0133]
如图8所示则给出了2016年5月下的所选海区范围内的20m水深处的温度数据的趋势面拟合结果。观察拟合曲面可以看出，拟合的趋势面的温度值由低纬度向高纬度呈下降趋势，这与海水温度数据的空间变化特征吻合。因此，可以认为趋势面拟合有效剔除了空间分布对温度数据的影响，而通过趋势面拟合，采用剩余量进行时空插值是具有合理性的。
[0134]
以上为针对argo温盐数据的预处理过程，而对于tao温盐数据，由于其时空分布规律较为稳定，具有固定的空间位置、稳定的时间观测频率和相对统一的深度层。同时tao温盐数据只是作为辅助数据，用于辅助计算拟合时间变异函数。因此在预处理过程，只需要筛选出符合时空范围的tao温盐数据，同时剔除无效数据，即可满足课题的要求。
[0135]
2)时空变异函数建模求解
[0136]
该建模求解过程可分为空间变异函数拟合、时间变异函数拟合、积合式模型构建三部分。
[0137]
在进行空间变异函数拟合和时间变异函数拟合前，首先需要依据时空数据的特点，使用恰当的变异函数理论模型及拟合方法。对于变异函数理论模型，常用的有：纯块金效应模型、高斯模型、球状模型和指数模型等，分别反映了不同的变异函数曲线的变化特征。拟合方法则有最小二乘法、加权回归法。
[0138]
实际上，根据对使用的时空数据的分析，大部分的变异函数的拟合都可以采用球状模型进行拟合。因此，本实施例采用球状模型进行空间变异函数、时间变异函数的拟合，其表达式如下：
[0139][0140]
其中，c0为块金值，c为拱高，a为变程。
[0141]
下面以经度148.968
°
、纬度26.362
°
、时间2016年5月5日、深度层为20m的点为待插值的数据点，进行温度的时空插值。
[0142]
(1)时间变异函数拟合
[0143]
由于tao温盐数据具有固定空间位置，且观测周期稳定为1天，非常便于时间变异函数的构建和拟合，故采用tao温盐数据进行时间变异函数的拟合。首先是读取选择海区附近的tao数据，即上述的预处理过程，时间范围为5月份的数据，之后按天为单位利用式计算各个固定点的变异函数值，并对相同变程(即曲线的自变量值，数据间的时间间隔值)的数值进行加权后，基于最小二乘法用球状模型拟合得到各固定点的时间变异函数γ
t
(si,h
t
)。
[0144]
对于拟合得到的各个固定点的时间变异函数，利用式按与待插值点的空间距离进行反距离加权参数求和，得到最终的时间变异函数表达式，所采用的反距离权重系数λi的计算公式具有如下表达式：
[0145][0146]
其中，hj为各个固定点到插值点的空间距离。
[0147]
对于待插值的数据点处的温度数据，拟合得到的时间变异函数如下：(块金指数c0＝0.00700637，变程a＝16，拱高c＝0.01252363)
[0148][0149]
(2)空间变异函数拟合
[0150]
数据采用预处理后的argo温盐数据(空间范围即为提及的经纬度范围，时间范围
则是2016年5月)，由于单个时间段，即单天内的温盐数据的数据量较少(平均只有5个数据点)，这里为提高可用于拟合的数据量，采用多时段叠置的空间变异寒素拟合方法，将多个时间段内的数据叠加，计算叠加后数据之间的空间变异函数数值以及相对应的地理距离，最后综合计算的函数值，对函数值按对应的地理距离所处区间进行划分(这里以每50km进行一个区间划分)，对区间内的变异函数值进行加权平均后，得到单个地理距离下的空间变异函数值，再基于最小二乘法采用球状模型拟合得到空间变异函数曲线。其中，两个数据点间的地理距离计算采用haversine公式，表达式如下所示：
[0151][0152]
其中，r为地球半径，一般为6371km，为两数据点的纬度，δλ为两点经度差值。
[0153]
同样对于所选时空范围下的深度层20m的argo温度数据，通过上述方法拟合得到的空间变异函数表达式如下：(块金指数c0＝0.430662，变程a＝247，拱高c＝0.05976)
[0154][0155]
(3)时空变异函数拟合
[0156]
基于上面拟合得到的时间变异函数和空间变异函数结果，可知时间变异函数基台值c
t
(0)＝0.01953，空间变异函数基台值cs(0)＝0.430662；时空变异函数的基台值则可以取时空样本的最大值c
s,t
(0,0)＝0.666918，之后代入积合式模型的求解公式，可以得到参数k1＝-25.7675,k2＝1.5032,k3＝12.0971，最后化简得到的时空变异函数模型表达式如下：
[0157]
γ(hs,h
t
)＝γs(hs)+γ
t
(h
t
)+25.7675γs(hs)γ
t
(h
t
)
[0158]
3)时空kriging插值
[0159]
其中各参与插值的数据点的选择，使用了待插值点的空间范围方圆10
°
以内的、时间间隔
±
3天内的数据点，si,ti分别为各参与插值的数据点的空间位置和所处时间。
[0160]
选取好参与插值的数据点后，则是对于权重系数λi的求解，依据式构建权重系数λi关于时空变异函数的方程组，并转化为如式所示的矩阵形式进行求解，通过计算矩阵γ
ij
和列向量γ
0i
的值，并进行矩阵逆运算和矩阵乘积，即可求解得到列向量λi的值，也就是权重系数λi的值。
[0161]
针对该插值点的时空插值，共使用了插值点时空范围附近的28个数据点，并计算得到了各个数据点的权重系数，最后插值得到了该点的数据值，并与该点的趋势量相加后，得到了该插值点的插值温度值为25.2989℃。
[0162]
对于插值结果的准确性，还需要对插值的精度评价进行衡量，由于一般温度变化具有较为明显的特征，本发明精度评价主要采用温度数据进行。以选定的经纬度范围的海区作为数据插值精度评价的空间范围，时间范围则是2016年5月1日到2016年5月30日内。
[0163]
插值精度评价的过程按照时空kriging插值的流程进行(如图3所示)，并采用“留一法”进行验证，即按日期进行逐次查询处于时空范围内的数据点，每次都预留一个观测数据点作为待插值的点进行时空kriging插值。之后将插值结果同实际观测结果进行比对计算实际误差，依据误差的大小进行评估。各个预留数据点的误差大小及空间分布情况如图9、图10所示。
[0164]
进一步，本技术基于不同评价指标计算差值结果和实际结果间的误差大小，包括平均绝对误差、均方根误差和平均绝对百分比误差。其评价结果如表2所示，表2
[0165]
参数类型maermsemape温度(℃)0.91561.43433.62％盐度(psu)0.16840.24510.48％密度(kg/m3)0.34900.46690.034％声速(m/s)2.24293.62241.5％
[0166]
综合表内的各类指标下的计算结果可以看出，该插值法具有较高的插值精度，可以满足课题要求的相对误差小于15％的指标。
[0167]
同时采用该时空插值方法生成温盐数据，能够在装备试验过程中，依据试验的需求生成指定时空位置下的海水温度、盐度数据，依据工程模型又可以计算得到对应的海水密度、声速信息，为装备试验，尤其是水下航行设备运行过程的试验，提供需要的参数信息。
[0168]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0169]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

技术特征：
1.一种基于时空插值的高分辨温盐数据生成方法，其特征在于，包括，获取温盐的时空监测数据，并对所述时空监测数据的总体分布形式进行平稳性检验，提取满足平稳性的温盐时空数据；根据所述温盐时空数据中的时间数据获得时间变异函数，根据所述温盐时空数据中的空间数据获得空间变异函数；根据所述时间变异函数和所述空间变异函数，构建时空变异函数；利用时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据，所述时空kriging插值法的表达式为：其中，z
*
(s
i
,t
i
)表示插值后的点，z(s
i
,t
i
)表示样本点i，(s
i
,t
i
)表示样本点i的时空域坐标，λ
i
(i＝1,2,...,n)表示样本点i的权重系数。2.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，采用二阶平稳假设对时空监测数据进行平稳性分析，其中，约束条件如下：所述时空监测数据在所处的研究区域内具有恒定的数学期望；所述时空监测数据在所处的研究区域内存在协方差函数，且所述协方差函数平稳，数值只与位移量有关。3.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，采用平稳假设或本征假设对所述时空监测数据的总体分布形式进行平稳性分析。4.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，所述时空监测数据的总体分布形式通过统计直方图法和/或趋势面分析法进行分析。5.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，所述时间变异函数和所述空间变异函数，采用如下理论模型和拟合方法获得；所述理论模型包括有基台值模型和/或无基台值模型，其中，有基台值模型包含纯块金效应模型、高斯模型、和/或球状模型；无基台值模型包含线性无基台值模型；所述拟合方法包括最小二乘法和/或加权回归法。6.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，采用积合式模型构建时空变异函数，构建方法如下：γ(h
s
,h
t
)＝(k1c
t
(0)+k2)γ
s
(h
s
)+(k1c
s
(0)+k3)γ
t
(h
t
)-k1γ
s
(h
s
)γ
t
(h
t
)式中，γ
s
(h
s
,h
t
)是所构建的时空变异函数，h
s
和h
t
为时空域中样本点坐标；c
st
(0,0)、c
s
(0)、c
t
(0)分别是时空变异函数、空间变异函数和时间变异函数的基台值，γ
s
(h
s
)为空间变异函数，γ
t
(h
t
)为时间变异函数。7.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，使z
*
(s,t)为z(s,t)的无偏估计量并具有最小的估计方差，并根据如下公式对所述λ
i
求解，
式中，μ为估计方差，γ((s
i
,t
i
),(s
j
,t
j
))为时空变异函数，其中，(s,t)为时空域中样本点坐标的展开形式。8.根据权利要求1所述的一种基于时空插值的高分辨温盐数据生成方法，其特征在于，所述时空监测数据包括时间监测数据和空间监测数据；所述时间监测数据经多个观测点，通过如下公式获得：式中，h
t
为对应的时间间隔距离，z(s,t
i
)为样本点(s,t
i
)的观测值，z(s,t
i
+h
t
)是与z(s,t
i
)时间距离相隔h
t
的样本观测值；n为同一观测点不同时刻的观测点的总数；所述空间监测数据通过对多个时间切片的采样点数进行叠置获得，叠置公式如下：式中，m表示时间切片的个数，k为空间观测点的数量，h
s
为时间间隔距离，j代表第j个时间切片。9.一种基于时空插值的高分辨温盐数据生成装置，其特征在于，包括，温盐时空数据获取单元，用于获取温盐时空监测数据，并提取满足平稳性要求的温盐时空数据；时间和空间变异函数拟合单元，用于根据温盐时空数据中的时间数据拟合得到时间变异函数，根据温盐时空数据中的空间数据拟合得到空间变异函数；时空变异函数构建单元，用于根据所述时间变异函数和空间变异函数，构建时空变异函数；高分辨率温盐数据生成单元，用于通过时空kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据。10.根据权利要求9所述的一种基于时空插值的高分辨温盐数据生成装置，其特征在于，还包括数据叠置拟合单元，用于对多时段、多空间的监测数据进行拟合。

技术总结
本发明公开了一种基于时空插值的高分辨温盐数据生成方法及装置，方法包括获取温盐的时空监测数据，对所述时空监测数据的总体分布形式进行检验，提取满足平稳性的温盐时空数据；根据所述温盐时空数据中的时间数据得到时间变异函数，根据所述时空数据中的空间数据得到空间变异函数；根据所述时间变异函数和所述空间变异函数，构建时空变异函数；利用时空Kriging插值法对所述时空变异函数进行插值，得到高分辨率的温盐网格数据。本发明中公开的方法或装置可同时考虑时间维度和空间维度的观测数据，并基于时空分布离散、稀疏的观测数据生成标准化的、具有高分辨率的温度和盐度网格数据。格数据。格数据。


技术研发人员：林连雷 于航懿 杨昌达
受保护的技术使用者：哈尔滨工业大学
技术研发日：2023.06.19
技术公布日：2023/10/11