本发明属于海洋灾害监测技术领域，尤其涉及一种台风中心电离层tec异常的滑动三次曲线探测方法。
背景技术：
台风作为对流层内的一种强天气系统，能够在大气中激发重力波并通过能量、动量输送至电离层高度，从而造成电离层扰动(黄朝松等，1992)。早在1960年，hines(hines，1960)就提出了高层大气存在重力波的理论，同时证明了电离层与中低层大气之间存在耦合机制。一些学者利用卫星和电离层多普勒数据与台风资料进行分析和对比，发现了台风对电离层有一定的扰动作用(fabriceetal.，2002；garneretal.，2008；vanina-dartetal.，2011；余涛等，2010；高泽等，2017；邓忠新等，2012)。肖赛冠等(2006)对1988年和1990年两次强台风影响期间的电离层多普勒记录及相应的台风资料进行了细致的相关分析，研究表明，在这两次台风影响期间，电离层形态中除有明显的波状扰动(中尺度声重波)出现外，还显示了声重波在激发电离层不规则结构方面的种子作用。xiao等(2007)收集并分析了1987到1992的24个强台风完整的记录和与之相应的电离层高频多普勒数据。通过对连续性高频多普勒观测数据的分析发现这些反应的时空演变及共同特点。研究结果表明，台风确实对对流层的波状扰动发挥了重要作用尤其是当台风在大陆沿海登陆时更为明显。isaev等(2007)利用cosmos-1809卫星数据分析了电离层与台风的相关性。结果表明，低纬电离层的电场增强与热带风暴和台风的形成和演化有关。tao等(2010)研究了2007年在厦门登陆的3次台风事件对电离层参数f0f2的影响，结果显示，3次台风在登陆前后的时间段内，电离层参数f0f2均出现扰动。随着gps技术的发展，越来越多的学者开始利用gps卫星数据计算电离层tec，从而研究台风对电离层的影响(kunitsynetal.，2011；perevalovaetal.，2011，张东和等，2000；张小红等，2013)。毛田等(2010)利用50余个gps台站的观测资料，研究了台风“麦莎”对电离层tec的影响，分析表明，台风在登陆前已经影响到电离层，tec出现增大趋势。王心智等(2016)利用国际igspimo和bjfs站的gps数据计算了utor台风期间的电离层vtec值，通过分析发现，pimo站的vtec值在台风utor发生期间，出现了明显的增加。
现有技术中关于台风影响下的电离层tec异常探测方法具有如下缺点：
(1)观测资料的局限性：现有技术基本都是利用设置在陆地的gps接收机接收gps卫星信号来反演电离层tec变化，或者利用设置在陆地的电离层高频多普勒台站来获取电离层资料电离层tec异常情况。然而，无论是gps接收机还是电离层高频多普勒台站都位于陆地，而台风一般都发源于热带海面，那里温度高，大量的海水被蒸发到了空中，形成一个低气压中心。随着气压的变化和地球自身的运动，流入的空气也旋转起来，形成一个逆时针旋转的空气漩涡，这就是热带气旋。只要气温不下降，这个热带气旋就会越来越强大，最后形成了台风。台风在海面的生成过程中对陆地的影响较小，仅仅是台风接近陆地或登陆时，陆地上的观测仪器才会观测到相应的数据变化。因此，对于台风生成及发展阶段，无论是陆地上的gps接收机还是电离层高频多普勒台站都难以准确地探测电离层tec异常。
(2)探测方法的局限性：由于台风活动期间，会出现电离层tec异常，因此对电离层tec异常探测具有十分重要的意义。目前，现有技术中对于电离层tec异常值探测的主要方法有平均值法、中位数法、滑动均值法、滑动四分位法等。
其中，平均值法和中位数法属于静态探测方法，其对照及检测的资料在时间上已经覆盖了时间点之前和之后的效应，未能剔除样本数据中的异常数据，使得计算的背景值不合理，从而导致探测的异常值不可靠。滑动均值法是在平均值法的基础上加以改进，它通过选取适当的滑动窗口，计算该窗口下的均值和标准差，以均值为背景值，k倍标准差作为异常值检测依据。滑动四分位法与滑动均值法类似，首先选取适当的滑动时间窗口，将该窗口下相同时刻的tec值从小到大排列并分为四等分，其等分点依次表示为q1、q2、q3，则四分位距值iqr＝q3-q1，q2作为该窗口下的tec背景值，取相应的k倍四分位距值作为异常检验边界。滑动均值法和滑动四分位法相比平均值法和中位数法，能够较为准确地计算出背景值，进而依据背景值计算出电离层tec的上下限值，从而检测异常值的存在。但由于台风影响下的电离层tec数据通常具有非线性和非平稳的复杂信号特征，导致滑动均值法和滑动四分位法计算的背景值的精度有限，因此，探测电离层tec异常结果的可靠性也不高。
技术实现要素：
本发明的目的在于提出一种台风中心电离层tec异常的滑动三次曲线探测方法，以解决现有观测资料局限性和探测方法局限性的问题。
本发明为了实现上述目的，采用如下技术方案：
台风中心电离层tec异常的滑动三次曲线探测方法，包括如下步骤：
s1.搜索台风中心电离层tec网格点
根据日本气象厅发布的台风观测数据，以及欧洲定轨中心提供的电离层tec网格点，搜索台风风速最大时，距离台风中心最近的tec网格点；
s2.构建电离层tec时间序列
以距离该台风中心最近的tec网格点为观测点，当探测该tec网格点第y天的tec背景值时，选择第y天之前d天和当天的tec数据构成一个序列；
其中，d被称为滑动时间窗口长，滑动时间窗口d为27天；
第y天和前27天电离层tec序列，一共构成28天的电离层tec序列{xz}，z＝1,2,…,28；
由于欧洲定轨中心提供的全球电离层地图时间分辨率是2小时，因此，构成的电离层tec时间序列总长度数目n＝12*28＝336；
s3.电离层tec时间序列的分离
s3.1.搜索电离层tec时间序列极值点
设电离层tec时间序列的第a个值为teca，在一个邻域t内的任意点，都有tecb<teca，则teca是电离层tec时间序列的一个极大值；同理，若在一个邻域t内的任意点，都有tecb>teca，则teca是电离层tec时间序列的一个极小值；
其中，t为tec时间序列的平均周期，tecb为teca一个邻域t内的任意值；
s3.2.极值点线段的分段
按照时间顺序分别将电离层tec极大值点和极小值点进行编号，同时分别将电离层tec的极大值和极小值连接起来，构成由多个极值点组成的直线段；
对数据进行分段时，按照编号由前向后的顺序每五个数据分为一段，直至将所有数据都分成不同的线段，若最后数据不足五个时，则与前一条线段组成线段；
s3.3.最小二乘三次曲线拟合
某极值点线段数据的三次拟合曲线函数为：wt＝a1+b1t+c1t2+d1t3，t＝-2,-1,0,1,2，利用最小二乘原理求解三次拟合曲线的系数a1、b1、c1和d1，从而得到各线段的曲线方程；
s3.4.电离层tec时间序列的分离
将由极大值点构成的三次拟合曲线和由极小值点构成的三次拟合曲线求取平均值，得到电离层tec时间序列的子序列；
s3.5.重复步骤s3.1～s3.4，直到无法构成极值点线段，即无法再分离出子序列；
s4.探测各电离层tec时间子序列的背景值和上下限值，将每天同一时刻的tec子序列|tij|，按从小到大排列得到序列t′kj，然后分成四等分，其等分点依次表示为q1、q2和q3；
其中，i表示天数，i＝1,2,…,28，j表示时段，j＝1,2…,12；
k表示排序后的天数，k＝1,2,…,28；
公式如下：
式中，q1、q2、q3分别为第一、第二、第三四分位数；
q1表示由j时段排在第7位和第8位的tec值t′7j和t′8j计算的平均值；
q2表示由j时段排在第14位和第15位的tec值t′14j和t′15j计算的平均值；
q3表示由j时段排在第21位和第22位的tec值t′21j和t′22j计算的平均值；
iqr＝q3-q1(2)
iqr为四分位差，它是第三四分位数q3与第一四分位数q1的差；
则第i天第j时段的电离层tec异常的背景值为q2，其上下限值为：
εij＝3iqr(4)
其中，l1j和l2j表示第i天第j时段电离层tec异常的上下限；
εij表示以3倍的四分位差iqr作为误差限值；
s5.重构电离层tec的背景值和上下限值
步骤s4得到的是步骤s3分离后的各个电离层tec时间子序列的背景值和上下限值，而整个时间序列的背景值和上下限值需要将各个子序列的背景值和上下限值合并，即求和；
假设某电离层tec被分离成s个子时间序列，由各子序列的背景值求和得到总的背景值，即第i天第j时段的电离层tec的背景值为q′2(i,j)，计算公式如公式(5)所示：
式中，f表示被分离的各子时间序列的序号，s表示被分离的子时间序列的个数；
定义第i天第j时段的电离层tec的上下限值为l′1j和l′2j；
则l′1j和l′2j的计算公式分别如公式(6)和公式(7)所示：
当电离层tec值超过上限值或下限值，即表示电离层tec出现异常。
本发明具有如下优点：
1.解决观测资料局限性的问题：
由于gps接收机和电离层高频多普勒台站都位于陆地，而台风一般则在海洋深处生成，为更好地研究台风生成阶段时其对电离层tec的扰动，本发明利用了欧洲定轨中心根据igs站和各国gps站得到的以2个小时为间隔形成的50×2.50的tec格网图数据，以台风移动路径中风速最大时最近的tec网格点为研究对象，探测电离层tec异常情况。
2.解决探测方法局限性的问题：
由于电离层tec异常值探测的主要方法(包括平均值法、中位数法、滑动均值法、滑动四分位法)，存在着探测异常值不可靠或者计算的背景值的精度有限等缺点，本发明利用了滑动窗口思想解决了静态探测方法(平均值法和中位数法)未能剔除样本数据中异常数据的缺点，利用信号分离思想，将复杂信号分离成简单信号的集合，解决了滑动探测方法(滑动均值法和滑动四分位法)对非线性、非平稳信号异常值探测精度不高的问题。
附图说明
图1为本发明中台风中心电离层tec异常的滑动三次曲线探测方法流程图；
图2为风尼伯特移动轨迹及电离层tec格网点示意图；
图3为滑动均值法探测电离层tec上下限值示意图；
图4为滑动四分位法探测电离层tec上下限值示意图；
图5为滑动三次曲线分离法探测电离层tec上下限值示意图；
图6为滑动均值法、滑动四分位法以及本发明中的滑动三次曲线分离法在探测tec背景值相对精度对比示意图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
结合图1所示，台风中心电离层tec异常的滑动三次曲线探测方法，包括如下步骤：
s1.搜索台风中心电离层tec网格点
根据日本气象厅发布的台风观测数据，以及欧洲定轨中心(thecenterfororbitdeterminationineurope，简称code)提供的电离层tec网格点。
以2016年台风尼伯特(nepartak)事件为例，表1为台风“尼伯特”生成主要阶段相关信息。根据表1中的信息，台风尼伯特在2016年7月6日12时，其中心位于东经127度，北纬20.1度的太平洋上，此时台风中心风速最大，达到120节(knot)。
因此，可以搜索到欧洲定轨中心提供的tec网格数据中，距离该台风中心最近的tec网格点为东经125度，北纬20度，如图2所示。
表1台风“尼伯特”生成主要阶段相关信息
s2.构建电离层tec时间序列
以距离该台风中心最近的tec网格点为观测点，当探测该tec网格点第y天的tec背景值时，选择第y天之前d天和当天的tec数据构成一个序列。
其中，d被称为滑动时间窗口长，滑动时间窗口d为27天。
第y天和前27天电离层tec序列，一共构成28天的电离层tec序列{xz}，z＝1,2,…,28。
由于欧洲定轨中心提供的全球电离层地图时间分辨率是2小时，因此，构成的电离层tec时间序列总长度数目n＝12*28＝336。
s3.电离层tec时间序列的分离
s3.1.搜索电离层tec时间序列极值点
设电离层tec时间序列的第a个值为teca，在一个邻域t内的任意点，都有tecb<teca，则teca是电离层tec时间序列的一个极大值；同理，若在一个邻域t内的任意点，都有tecb>teca，则teca是电离层tec时间序列的一个极小值。
其中，t为tec时间序列的平均周期，tecb为teca一个邻域t内的任意值。
s3.2.极值点线段的分段
按照时间顺序分别将电离层tec极大值点和极小值点进行编号，同时分别将电离层tec的极大值和极小值连接起来，构成由多个极值点组成的直线段。
编号顺序为从1开始，直至最后一个极值点。
对数据进行分段时，按照由前向后的顺序，每个数据分为一段，直至将所有数据都分成不同的线段，若最后数据不足五个时，则与前一条线段组成线段。
s3.3.最小二乘三次曲线拟合
某极值点线段数据的三次拟合曲线函数为：wt＝a1+b1t+c1t2+d1t3，t＝-2,-1,0,1,2，利用最小二乘原理求解三次拟合曲线的系数a1、b1、c1和d1，从而得到各线段的曲线方程。
s3.4.电离层tec时间序列的分离
将由极大值点构成的三次拟合曲线和由极小值点构成的三次拟合曲线求取平均值，得到电离层tec时间序列的子序列。
s3.5.重复步骤s3.1～s3.4，直到无法构成极值点线段，即无法再分离出子序列。
s4.探测各电离层tec时间子序列的背景值和上下限值，将每天同一时刻的tec子序列|tij|，按从小到大排列得到序列t′kj，然后分成四等分，其等分点依次表示为q1、q2和q3。
其中，i表示天数，i＝1,2,…,28，j表示时段，j＝1,2…,12。
k表示排序后的天数，k＝1,2,…,28。
公式如下：
式中，q1、q2、q3分别为第一、第二、第三四分位数。
q1表示由j时段排在第7位和第8位的tec值t′7j和t′8j计算的平均值。
q2表示由j时段排在第14位和第15位的tec值t′14j和t′15j计算的平均值。
q3表示由j时段排在第21位和第22位的tec值t′21j和t′22j计算的平均值。
iqr＝q3-q1(2)
iqr为四分位差，它是第三四分位数q3与第一四分位数q1的差。
则第i天第j时段的电离层tec异常的背景值为q2，其上下限值为：
εij＝3iqr(4)
其中，l1j和l2j表示第i天第j时段电离层tec异常的上下限。
εij表示以3倍的四分位差(iqr)作为误差限值。
s5.重构电离层tec的背景值和上下限值
步骤s4得到的是步骤s3分离后的各个电离层tec时间子序列的背景值和上下限值，而整个时间序列的背景值和上下限值需要将各个子序列的背景值和上下限值合并，即求和。
假设某电离层tec被分离成s个子时间序列，由各子序列的背景值求和得到总的背景值，即第i天第j时段的电离层tec的背景值为q′2(i,j)，计算公式如公式(5)所示：
式中，f表示被分离的各子时间序列的序号，s表示被分离的子时间序列的个数。
定义第i天第j时段的电离层tec的上下限值为l′1j和l′2j。
则l′1j和l′2j的计算公式分别如公式(6)和公式(7)所示：
本发明通过将原始时间序列进行分离，得到各子时间序列的背景值和上下限值，然后再将各子时间序列的背景值以及上下限值分别求和，得到了原始时间序列的背景值和上下限值，即重构出电离层tec的背景值和上下限值。当电离层tec值超过上限值或下限值，即表示电离层tec出现异常，从而实现了电离层tec异常值的探测。
分别利用滑动均值法、滑动四分位法以及本发明中的滑动三次曲线分离法探测电离层tec异常上下限值及原始观测值，如图3～图5所示。
通过对比不难看出，本发明方法相比于滑动均值法、滑动四分位法，具有如下优点：
1.探测的电离层tec背景值平均相对精度更高；
2.探测的电离层tec背景值平均绝对精度更高。
具体表现在：
对搜索到的电离层tec网格点的电离层tec分别采用滑动均值法、滑动四分位法以及滑动三次曲线分离法来探测该点7天的tec背景值，并比较其精度。
采用相对精度prel和绝对精度pabs来进行评判，其定义如下：
式中，l为电离层tec时间序列值的序号，m为电离层tec时间序列值的个数；ibgv为tec背景值，iigs表示igs发布的tec观测值，σ为探测的背影值与tec观测值比较得到的中误差，δave表示探测的背影值与tec观测值比较得到的平均残差。
用以上三种方法探测电离层tec相对精度结果如下图6所示。
从图6可以看出，三种方法的相对精度变化趋势基本一致，说明三种方法的基本思路是一致的，只是在具体细节上有所不同。
从不同的时间段可以看出，滑动三次曲线分离法在总体上比另外两种方法的相对精度略高，仅仅在第2、3和6天的部分时间段比另外两种方法略低。
为了进一步比较三种方法探测电离层tec背景值的精度，分别对这三种方法探测tec背景值相对精度进行统计，结果如表2。
表2三种方法每天探测tec背景值相对精度平均值
从表2中可以看出，整体而言，三种方法每天的平均精度变化趋势基本一致。滑动三次曲线分离法相对精度除在第2天与另外两种方法精度相当，在第5天比另外两种方法略低以外，其它天数都高于另外两种方法。在第1天和第2天，滑动三次曲线分离法比另外两种方法的相对精度都高出约3％。从每天相对精度的平均值来看，滑动三次曲线分离法(87.39％)比另外两种方法(86.49％、85.94％)的相对精度分别高出约1％～2％。
三种方法探测tec背景值的绝对精度如表3所示。
表3三种方法每天探测tec背景值绝对精度比较
从表3中可以看出，三种方法在每天的背景值的绝对精度变化幅度都较大。三种方法，在不同的日期里，精度最高的是滑动四分位距法在第5天的绝对精度最高，为0.70tecu，而精度最低的仍然是滑动四分位距法在第7天的绝对精度最低，为14.91tecu，说明滑动四分位距法计算的背景值的精度变化幅度最大。滑动三次曲线分离法在第1、2和3天的绝对精度都最高，而在第4～6天的绝对精度都不是最低的，7天绝对精度的平均值来看，其平均精度是最高的，为5.82，而另外两种方法分别是6.41和7.18，分别高出10.14％和23.37％。
对三种方法探测tec背景值残差绝对值δ分四类，即：
δ≤1.5tecu、1.5tecu＜δ≤3tecu、3tecu＜δ≤5tecu和δ≥5tecu；
进行统计，结果如表4。
表4三种方法探测tec背景值残差绝对值分类百分比统计
从表4中可以看出，滑动均值法和滑动四分位距法探测tec背景值残差绝对值随着残差绝对值变大，所占的百分比呈现先减小后增大，而滑动三次曲线分离法的变化趋势是逐渐减小的。在δ≤1.5tecu这个区间段，滑动三次曲线分离法与滑动均值法所占的比例一样，而两种方法比滑动四分位距法略高2％。而在δ≥5tecu这个区间段，滑动三次曲线分离法比滑动均值法略低约1％，比滑动四分位距法略低约5％。
由此可以看出，本发明方法探测的背景值残差位于小误差区间(δ≤5tecu)的比例比滑动均值法和滑动四分位距法都高，因此，效果更好。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。