一种海洋气候研究中的滤波方法及其应用

1.本发明涉及海洋及大气数据处理技术领域，尤其涉及一种海洋气候研究中的滤波方法及其应用。


背景技术：

2.海洋和大气中普遍存在多种尺度的运动过程。通过空间尺度的区分可以分为行星尺度运动、大尺度运动、中尺度运动、小尺度运动及微尺度运动等。其也能够通过时间尺度进行区分。通常来说，空间尺度越大的运动过程，其对应的时间尺度也越大。目前在此领域获得广泛应用的滤波方法为滑动平均法、粗粒化法和高斯函数或窗函数法，通过对傅里叶变换得到的频域或波数域观察发现，所得到的低通滤波结果在高频或高波数部分仍具有能量分布，通过原始数据减去低通滤波的结果即为高通滤波结果。同样地，通过对傅里叶变换得到的频域或波数域观察发现，所得到的高通滤波结果在低频或低波数部分也具有能量分布。尤其是高通滤波所得结果的总能量小，采用这些方法得到的高通滤波结果在小于区分尺度对应的部分能量大小与大于区分尺度对应的部分能量大小量级相当。这就为后续的计算研究引入了很大的误差，进而影响了分析该物理量对地球海洋气候系统的作用结果的准确性。因此，现亟需一种用于海洋气候研究不同时间、空间尺度，引入在频域或波数域的截断处理方式，以及对截断结果进行判断提升尺度区分准确度的滤波方法及应用。


技术实现要素：

3.本发明的目的在于克服上述现有技术中的不足，旨在提供一种用于海洋气候研究不同时间、空间尺度，引入在频域或波数域的截断处理方式，以及对截断结果进行判断提升尺度区分准确度的滤波方法及应用。
4.为达到上述目的，本发明是通过下述技术方案予以实现的：
5.一种海洋气候研究中的滤波方法，包括如下步骤：
6.(1)对时间或空间原始数据x求傅里叶变换得到频域或波数域上的结果
7.(2)利用误差函数erf、激活函数tanh或sigmoid对频域或波数域上的进行截断处理，分别得到低通滤波频域结果和高通滤波频域结果
8.其中，误差函数erf的表达式为
9.激活函数tanh的表达式为
10.激活函数sigmoid的表达式为
11.误差函数erf对时空数据进行尺度区分，即在三维频率-波数域上截断，得到的低通滤波部分和高通滤波部分表达式分别表示为：
[0012][0013][0014]
激活函数tanh在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：
[0015][0016][0017]
激活函数sigmoid在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：
[0018][0019][0020]
其中，ω为频率，k和l分别为纬向和经向的波数分量。a为时间或空间区分尺度的倒数，b为截断边缘的倾斜度。p、q和r为函数参数，取值均为0或1，当p＝q＝r＝1时，此函数在三个方向上的截断大小相等；当p、q和r中任意一个等于0时，上式表示在二维情况下截断；若p＝0，上式表示在二维波数域上截断；若q＝0或r＝0，上式表示在二维频率-波数域上截断；
[0021]
(3)对低通和高通滤波结果的能量进行计算分析，为量化滤波质量，以a作为目标区分频率或波数，则理想低通滤波结果的能量大小为理想高通滤波结果的能量大小为定义有效部分表示滤波后在频域期望区分尺度范围之内的能量部分，误差部分表示滤波后在频域期望区分尺度范围之外的能量部分，可以得到以下四个比例：
[0022]
①
低通滤波有效比例为：
[0023]
②
低通滤波误差比例为：
[0024]
③
高通滤波有效比例为：
[0025]
④
高通滤波误差比例为：
[0026]
其中，滤波有效比例越接近1，表示滤波效果越好，滤波误差比例越接近0，表示滤波效果越好，为得到一个表征整体滤波质量的参数，定义滤波质量指数为：
[0027]
q＝|1-valid
l
|+|1-validh|+|error
l
|+|errorh|
ꢀꢀ
(1-7)
[0028]
q为滤波质量指数；
[0029]
(4)对低通滤波频域或波数域结果进行傅里叶逆变换得到时间或空间域上大尺度的特征分布x
l
，x
l
即为低通滤波结果，通过原始数据x减去低通滤波后的结果x
l
得到小尺
度的特征分布xh，xh即为高通滤波结果；
[0030]
(5)当滤波结果满足步骤(3)中得到的频率或波数域上的滤波质量指数q小于20％，即0.2，且满足步骤(4)中经过逆变换后时间或空间域上的低通滤波结果特征分布不存在超过原始数据大小的异常值，则步骤(2)中的函数参数属于适合尺度区分的函数参数所在范围，如不能同时满足以上两个判定标准，则步骤(2)中的函数参数不属于适合尺度区分的函数参数所在范围，需要返回步骤(2)改变函数参数大小重新计算；
[0031]
(6)通过步骤(2)至步骤(5)，q越小得到适合尺度区分的函数参数越适合尺度区分，将得到的函数参数带入到步骤(2)的函数中，即可得到滤波结果；
[0032]
(7)将滤波结果带入地球海洋气候系统中，提高对地球海洋气候系统的作用结果的准确性。
[0033]
还包括，一种滤波方法的应用，基于一种海洋气候研究中的滤波方法，应用于尺度分析、能量串级、不同尺度运动过程的相互作用、涡旋特性研究、天气预报、滤除内波影响、揭示海洋条带结构。
[0034]
相对于现有技术，本发明具有以下有益效果：
[0035]
本发明中的滤波方法在频率域或波数域的可利用截断函数包括误差函数erf、激活函数tanh或sigmoid，截断的表达式简单直观。另外此方法能够克服现存尺度区分方法中滤波不干净的问题，通过频率或波数域和时间、空间域两个方面的评判标准进行判定，使得尺度区分结果更精确。
附图说明
[0036]
图1为本发明的滤波方法流程图；
[0037]
图2为实施例中得到的频域低通滤波部分和高通滤波部分示意图；
[0038]
图3为实施例中得到的低通有效比例、高通有效比例、低通误差比例、高通误差比例随滤波参数的变化情况示意图；
[0039]
图4为实施例中得到的滤波质量指数q随滤波参数的变化情况及其与背景技术中提及的现有方法的结果比对示意图；
[0040]
图5为实施例中得到的最优滤波质量指数q情况下低通滤波误差比例和高通滤波误差比例及其与背景技术中提及的现存方法的结果比对示意图；
[0041]
图6为实施例中时域低通滤波时间序列及其与背景技术中提及的现存方法的结果比对示意图。
具体实施方式
[0042]
下面将参考附图并结合实施例来详细说明本发明。
[0043]
如图1所示，一种海洋气候研究中的滤波方法，包括如下步骤：
[0044]
(1)对时间或空间原始数据x求傅里叶变换得到频域或波数域上的结果
[0045]
(2)利用误差函数erf、激活函数tanh或sigmoid对频域或波数域上的进行截断处理，分别得到低通滤波频域结果和高通滤波频域结果
[0046]
其中，误差函数erf的表达式为
[0047]
激活函数tanh的表达式为
[0048]
激活函数sigmoid的表达式为
[0049]
误差函数erf对时空数据进行尺度区分，即在三维频率-波数域上截断，得到的低通滤波部分和高通滤波部分表达式分别表示为：
[0050][0051][0052]
激活函数tanh在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：
[0053][0054][0055]
激活函数sigmoid在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：
[0056][0057][0058]
其中，ω为频率，k和l分别为纬向和经向的波数分量。a为时间或空间区分尺度的倒数，b为截断边缘的倾斜度。p、q和r为函数参数，取值均为0或1，当p＝q＝r＝1时，此函数在三个方向上的截断大小相等；当p、q和r中任意一个等于0时，上式表示在二维情况下截断；若p＝0，上式表示在二维波数域上截断；若q＝0或r＝0，上式表示在二维频率-波数域上截断；
[0059]
(3)对低通和高通滤波结果的能量进行计算分析，为量化滤波质量，以a作为目标区分频率或波数，则理想低通滤波结果的能量大小为理想高通滤波结果的能量大小为定义有效部分表示滤波后在频域期望区分尺度范围之内的能量部分，误差部分表示滤波后在频域期望区分尺度范围之外的能量部分，可以得到以下四个比例：
[0060]
①
低通滤波有效比例为：
[0061]
②
低通滤波误差比例为：
[0062]
③
高通滤波有效比例为：
[0063]
④
高通滤波误差比例为：
[0064]
其中，滤波有效比例越接近1，表示滤波效果越好，滤波误差比例越接近0，表示滤波效果越好，为得到一个表征整体滤波质量的参数，定义滤波质量指数为：
[0065]
q＝|1-valid
l
|+|1-validh|+|error
l
|+|errorh|
ꢀꢀ
(1-7)
[0066]
q为滤波质量指数；
[0067]
(4)对低通滤波频域或波数域结果进行傅里叶逆变换得到时间或空间域上大尺度的特征分布x
l
，x
l
即为低通滤波结果，通过原始数据x减去低通滤波后的结果x
l
得到小尺度的特征分布xh，xh即为高通滤波结果；
[0068]
(5)当滤波结果满足步骤(3)中得到的频率或波数域上的滤波质量指数q小于20％，即0.2，且满足步骤(4)中经过逆变换后时间或空间域上的低通滤波结果特征分布不存在超过原始数据大小的异常值，则步骤(2)中的函数参数属于适合尺度区分的函数参数所在范围，如不能同时满足以上两个判定标准，则步骤(2)中的函数参数不属于适合尺度区分的函数参数所在范围，需要返回步骤(2)改变函数参数大小重新计算；
[0069]
(6)通过步骤(2)至步骤(5)，q越小得到适合尺度区分的函数参数越适合尺度区分，将得到的函数参数带入到步骤(2)的函数中，即可得到滤波结果；
[0070]
(7)将滤波结果带入地球海洋气候系统中，提高对地球海洋气候系统的作用结果的准确性。
[0071]
本实施例针对26年内的aviso卫星数据进行处理，该数据的时间分辨率为1天，空间分辨率为1/4
°
，分析选中的位置处海洋表面高度变化情况，以将一年作为时间区分尺度为例，进行低通和高通滤波。
[0072]
(1)对该卫星数据一个位置点处约26年内每天的海洋表面高度数据时间序列η求傅里叶变换得到频率域上的结果
[0073]
(2)期望将一年，即365天，作为时间区分尺度，对原始数据进行低通和高通滤波。利用误差函数erf的表达式：
[0074][0075][0076]
对公式(1-1)和(1-2)中的参数q和r取0，对频域进行截断处理，其中公式(1-1)和(1-2)中参数a取1/365，并将实际数值5
×
104代入控制截断边缘倾斜度的参数b进行计算，分别得到低通滤波频率域结果和高通滤波频率域结果低通滤波部分和高通滤波部分表达式分别表示为：
[0077][0078][0079]
其中即认为是滤除了时间尺度小于365天的年内波动影响的低通滤波部分，是仅保留了时间尺度小于365天的年内波动影响的高通滤波部分，由此得到的频率域原始
数据结果低通滤波部分和高通滤波部分见附图2。
[0080]
(3)对低通和高通滤波结果的能量进行计算分析，为量化滤波质量，以1/365(/天)作为理想区分频率，则理想低通滤波结果的能量大小为理想高通滤波结果的能量大小为又有步骤2中得到的低通滤波频域结果和高通滤波频域结果定义有效部分表示滤波后在频域期望区分尺度范围之内的能量部分，误差部分表示滤波后在频域期望区分尺度范围之外的能量部分，可以得到以下四个比例：
[0081]
①
低通滤波有效比例为：
[0082]
②
低通滤波误差比例为：
[0083]
③
高通滤波有效比例为：
[0084]
④
高通滤波误差比例为：
[0085]
其中，滤波有效比例越接近1，表示滤波效果越好，滤波误差比例越接近0，表示滤波效果越好。在此种比例分析方法情况下，定义滤波参数为l’，其对应的频率域参数为a’＝1/l’。将变化的a’带入到公式(1-8)、(1-9)中，可以得到当期望区分尺度为365天时以上四种比例随滤波参数l’变化的结果，本方法得到的四种比例随滤波参数l’变化的结果如附图3所示。可以由附图3看出低通滤波有效比例和低通滤波误差比例随滤波参数改变而产生的变化较小，而高通滤波有效比例和高通滤波误差比例随滤波参数改变而产生的变化则较大。同理可以得到背景技术中提及现存常用的滑动平均法和汉宁窗方法的四种比例随滤波参数l’变化的结果，此处不单独将其结果附于附图中。为得到一个表征低通滤波和高通滤波整体滤波质量的参数，定义滤波质量指数q为：
[0086]
q＝|1-valid
l
|+|1-validh|+|error
l
|+|errorh|
ꢀꢀ
(1-7)
[0087]
q越小，滤波效果越好。
[0088]
(4)对低通滤波频域结果进行傅里叶逆变换可以得到时间域上低通滤波结果η
l
的特征分布。通过原始数据η减去低通滤波后的结果η
l
得到时间域上高通滤波结果ηh的特征分布。如图6所示，为时域低通滤波得到的时间序列η
l
及其与背景技术中提及现存常用的滑动平均法和汉宁窗方法的低通滤波结果比对，由于高通滤波时间序列ηh在图中不能很好地用肉眼评判结果，故不单独将ηh的结果列在附图中。
[0089]
(5)判定方法需同时满足步骤(3)中频率域上的滤波质量指数q小于20％，即0.2，且步骤(4)中经过逆变换后时间域上的低通滤波结果η
l
的特征分布不存在超过原始数据大小的异常值，则步骤(2)中的函数参数属于适合尺度区分的函数参数所在范围。
[0090]
本实施例中滤波质量指数q随滤波参数l’改变而产生的变化情况及其与背景技术中提及的现存常用的滑动平均法和汉宁窗方法的结果比对见附图4所示。滑动平均法和汉
宁窗方法的滤波质量指数q在所有滤波参数l’范围内均大于0.2。而本实施例的滤波质量指数q在滤波参数l’为340.9至444.8区间内均小于0.1，即此区间内的滤波质量指数q均在可接受范围20％内，并在滤波参数l’为381.3时滤波质量指数q达到最小值，即为本方法频率域的最优情况。最优滤波质量指数q情况下低通滤波误差比例和高通滤波误差比例及其与背景技术中提及的现存常用的滑动平均法和汉宁窗方法的滤波结果比对见附图5。本实施例中得到的低通和高通滤波误差比例均比另外两种方法小，尤其是高通滤波误差比例比另外两种方法小一个数量级，且在可接受范围内，能够满足频率域评判标准。观察时域低通滤波结果η
l
的特征分布图即附图6中的时间序列，其特征分布不存在超过原始数据大小的异常值，即本实施例得到的η
l
变化平滑且和其他两种方法变化趋势一致，能够满足时间域评判标准。
[0091]
同时满足以上两个判定标准，则步骤2中截断函数的函数参数的实际数值属于适合尺度区分的函数参数所在范围。如不能同时满足以上两个判定标准，则步骤2中截断函数的函数参数的实际数值不属于适合尺度区分的函数参数所在范围，需要返回步骤2改变函数参数的大小重新计算，直至通过附图5和图6发现选取的函数参数的实际数值同时满足以上两个判定标准，即可确定是适合进行尺度区分的函数参数。
[0092]
(6)通过步骤(2)至步骤(5)，可以得到适合以时间尺度为一年进行低通和高通滤波，即以365天为区分尺度的函数参数的大小，将函数参数带入到步骤(2)的公式(1-8)和(1-9)中，即可得到较为优良的尺度区分结果，即滤波结果，本实施例得到的较为优良的低通滤波结果见附图6中的时间序列η
l
。
[0093]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种海洋气候研究中的滤波方法，其特征在于：包括如下步骤：(1)对时间或空间原始数据x求傅里叶变换得到频域或波数域上的结果(2)利用误差函数erf、激活函数tanh或sigmoid对频域或波数域上的进行截断处理，分别得到低通滤波频域结果和高通滤波频域结果其中，误差函数erf的表达式为激活函数tanh的表达式为激活函数sigmoid的表达式为误差函数erf对时空数据进行尺度区分，即在三维频率-波数域上截断，得到的低通滤波部分和高通滤波部分表达式分别表示为：波部分和高通滤波部分表达式分别表示为：激活函数tanh在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：示为：激活函数sigmoid在进行截断时，三维情况下低通滤波部分和高通滤波部分表达式分别表示为：别表示为：其中，ω为频率，k和l分别为纬向和经向的波数分量。a为时间或空间区分尺度的倒数，b为截断边缘的倾斜度。p、q和r为函数参数，取值均为0或1，当p＝q＝r＝1时，此函数在三个方向上的截断大小相等；当p、q和r中任意一个等于0时，上式表示在二维情况下截断；若p＝0，上式表示在二维波数域上截断；若q＝0或r＝0，上式表示在二维频率-波数域上截断；(3)对低通和高通滤波结果的能量进行计算分析，为量化滤波质量，以a作为目标区分频率或波数，则理想低通滤波结果的能量大小为理想高通滤波结果的能量大小为定义有效部分表示滤波后在频域期望区分尺度范围之内的能量部分，误差部分表示滤波后在频域期望区分尺度范围之外的能量部分，可以得到以下四个比例：
①
低通滤波有效比例为：
②
低通滤波误差比例为：
③
高通滤波有效比例为：
④
高通滤波误差比例为：其中，滤波有效比例越接近1，表示滤波效果越好，滤波误差比例越接近0，表示滤波效果越好，为得到一个表征整体滤波质量的参数，定义滤波质量指数为：q＝1-valid
l
+1-valid
h
+error
l
+error
h
(1-7)q为滤波质量指数；(4)对低通滤波频域或波数域结果进行傅里叶逆变换得到时间或空间域上大尺度的特征分布x
l
，x
l
即为低通滤波结果，通过原始数据x减去低通滤波后的结果x
l
得到小尺度的特征分布x
h
，x
h
即为高通滤波结果；(5)当滤波结果满足步骤(3)中得到的频率或波数域上的滤波质量指数q小于20％，即0.2，且满足步骤(4)中经过逆变换后时间或空间域上的低通滤波结果特征分布不存在超过原始数据大小的异常值，则步骤(2)中的函数参数属于适合尺度区分的函数参数所在范围，如不能同时满足以上两个判定标准，则步骤(2)中的函数参数不属于适合尺度区分的函数参数所在范围，需要返回步骤(2)改变函数参数大小重新计算；(6)通过步骤(2)至步骤(5)，q越小得到适合尺度区分的函数参数越适合尺度区分，将得到的函数参数带入到步骤(2)的函数中，即可得到滤波结果；(7)将滤波结果带入地球海洋气候系统中，提高对地球海洋气候系统的作用结果的准确性。2.一种滤波方法的应用，基于一种海洋气候研究中的滤波方法，其特征在于：应用于尺度分析、能量串级、不同尺度运动过程的相互作用、涡旋特性研究、天气预报、滤除内波影响、揭示海洋条带结构。

技术总结
本发明公开了一种海洋气候研究中的滤波方法及其应用，包括如下步骤：(1)对时间或空间数据x求傅里叶变换得到频域或波数域上的结果(2)利用函数对频域或波数域上的进行截断处理，分别得到低通和高通滤波频域结果；(3)对低通和高通滤波结果的能量进行计算分析；(4)对低通滤波频域或波数域结果进行傅里叶逆变换可以得到时间或空间域上大尺度的特征分布；(5)当滤波结果满足步骤(3)中得到的频率或波数域上的滤波质量指数20％，则步骤(2)中的函数参数属于适合尺度区分的函数参数所在范围；(6)将滤波结果带入地球海洋气候系统中提高准确性。本发明旨在提供一种用于海洋气候研究中进行判断提升尺度区分准确度的滤波方法及应用。用。用。


技术研发人员：陈阳 陈儒
受保护的技术使用者：天津大学
技术研发日：2023.07.20
技术公布日：2023/10/19