技术领域
[0001] 本
发明属于空间物理领域，涉及极光物理，具体涉及一种基于行星际和地磁参数的极光卵边界位置预测方法。
背景技术
[0002] 极光卵是极光粒子沉降以磁极为中心在地球南北极地区形成的椭圆带状区域，是地球上来自太阳的
能量粒子的影响区域。极光卵能从全局尺度上对空间天气进行表征和预测，所以一个好的极光卵模型能有效地帮助理解空间天气并预测其对极区和次极区的影响。
[0003] 现有的对极光卵大小和位置的研究工作分为两大类：一种是定性分析，如讨论行星际
磁场(IMF)的By、Bz分量对极光卵边界位置的影响；另一种是单变量回归分析，回归量取某个地磁指数，如极光电集流指数(AE)，地磁扰动指数(KP)等。
[0004] 可事实上，一方面，极光活动是受所有行星际参数共同影响的，所以传统的一元回归分析或定量分析的方法不能充分表征极光卵的变化。另一方面，从物理源头来说，这些地磁活动指数(AE，Kp等)与极光是同一个层面的东西：它们都受
太阳风-
磁层-电离层耦合作用的影响。换句话说，这些地磁活动指数仅仅能从一个侧面来反映极光卵的属性，因为它们都受同一个“源”的影响，但并非因果关系。此外，以往研究还有一个明显缺点是将地磁指数按大小划分为几个区间，每个区间里只取其均值来做回归分析。这样就会出现一些非常相近的值被划分到了两个区间里，导致回归结果完全不同；此外，仅仅一个均值也不足以表征整个区间。
[0005] 卫星图像是同时确定整个(或绝大部分)极光卵边界的唯一途径。可是从海量的图像中人工确定极光卵边界是非常繁琐的，几乎不太可能完成。
发明内容
[0006] 本发明的目的在于提供一种基于行星际和地磁参数的UVI图像极光卵边界位置建模和预测方法，以克服上述
现有技术存在的
缺陷，本发明可用于对极光卵边界(包括高纬极向边界和低纬赤道向边界)位置预测，可以有效地帮助理解空间天气并预测其对极区和次极区的影响。
[0007] 为达到上述目的，本发明采用如下技术方案：
[0008] 基于行星际和地磁参数的极光卵边界位置预测方法，包括以下步骤：
[0009] 步骤1：选取Polar卫星拍摄到的UVI极光图像，并对UVI极光图像预处理，以消除噪声点、移除无效背景
像素；
[0010] 步骤2：对预处理后的UVI极光图像中的极光卵边界进行提取；
[0011] 步骤3：利用NASA OMNI数据，并根据每幅极光图像发生的UT时间，得到每个极光卵边界点对应的行星际和地磁参数；
[0012] 步骤4：对步骤2获得的极光卵边界数据结合步骤3获得的行星际地磁参数值进行筛选，去除无效边界点，构建有效数据集；
[0013] 步骤5：根据步骤4得到的有效数据集，用一元线性回归模型，对赤道向和极向极光卵边界位置进行建模，因变量取赤道向或极向极光卵边界的地磁纬度值，自变量取行星际和地磁参数；
[0014] 步骤6：根据步骤5建立的极光卵边界位置模型，对任意给定的行星际和地磁条件，能够预测出此条件下的赤道向和极向极光卵边界位置。
[0015] 进一步地，步骤1中对UVI极光图像预处理的方法为：
[0016] 步骤1.1：移除UVI极光图像中纬度低于50度的区域；
[0017] 步骤1.2：掩膜：根据UVI极光图像外边界呈椭圆状的特点，用直径为228像素的圆去截取经过步骤1.1处理后的UVI极光图像，并从中间把宽度限制到200像素，最后图像大小为228*200；
[0018] 步骤1.3：清除经步骤1.2处理后的UVI极光图像中像素值为负值的点；
[0019] 步骤1.4：将经过步骤1.3处理后的UVI极光图像中面积小于50个像素，极光强度大于整幅图像灰度均值加3倍方差的区域的像素值重置为整幅图像像素的灰度均值。
[0020] 进一步地，步骤2中对UVI极光图像中极光卵边界进行提取的方法为：
[0021] 步骤2.1：利用包含空间信息的模糊c均值聚类方法将UVI极光图像聚成6类；
[0022] 步骤2.2：根据聚类结果把UVI极光
图像分割为极光卵区和背景区域；
[0023] 步骤2.3：对聚类后出现的与周围区域类别不同的区域，用其周围邻域的类别信息对其进行平滑处理；
[0024] 步骤2.4：获取极光卵边界并将其映射到磁地方时-地磁纬度
坐标系下。
[0025] 进一步地，步骤2.2中根据聚类结果把UVI极光图像分割为极光卵区和背景区域的准则为：将强度最小的2个或3个聚类簇视为背景区域，其他则为极光卵区。
[0026] 进一步地，步骤3中所述的行星际参数包括行星际磁场三分量Bx、By、Bz和
太阳风动压Pdyn、太阳风速度Vp、太阳风
密度Np；地磁参数包括极光电集流指数AE、极盖指数PC。
[0027] 进一步地，步骤5中建立的模型为：
[0028] MLAT＝a0+a1·Bx+a2·By+a3·Bz+a4·Pdyn+a5·Vp+a6·Np
[0029] 其中，a0-a6表示此模型的回归系数，Bx,By,Bz为行星际磁场的三分量，Pdyn,Vp,Np分别为太阳风动压、速度和密度。
[0030] 进一步地，步骤5中建立的模型为：
[0031] MLAT＝b0+b1·Bx+b2·By+b3·Bz+b4·Pdyn+b5·Vp+b6·Np+b7·AE+b8·PC+b9·Ekl[0032] 其中，b0-b9表示此模型的回归系数，Bx,By,Bz为行星际磁场的三分量，Pdyn,Vp,Np分别为太阳风动压、速度和密度，AE为极光电集流指数，PC为极盖指数，Ekl为Kan-Lee
电场，它与行星际参数之间的关系为：
[0033] 其中θ取值方法为：
[0034] 若Bz>0,By>＝0，θ＝arctan(|By|/Bz)；
[0035] 若Bz<＝0,By>0，θ＝arctan(|Bz|/By+π/2)；
[0036] 若Bz<0,By<＝0，θ＝arctan(|By|/|Bz|+π)；
[0037] 若Bz>＝0,By<0，θ＝arctan(Bz/|By|+1.5π)；
[0038] 若Bz＝0,By＝0，θ＝0。
[0039] 进一步地，步骤5中建立的模型为：
[0040] MLAT＝c0+c1·VB2+c2·VPdyn+c3·Ekl+c4·Ez+c5·AE
[0041] 其中，c0-c5表示此模型的回归系数，VB2为行星际磁场磁能通量，VPdyn为太阳风2
动能通量，Ekl为Kan-Lee电场，Ez为南北电场，AE为极光电集流指数；VB 、VPdyn、Ez与行星际参数之间的关系为：
[0042]
[0043] VPdyn＝Vp·Pdyn
[0044] Ez＝By·Vp。
[0045] 进一步地，步骤5中建立的模型为：
[0046] MLAT＝d0+d1·Bx+d2·By+d3·Bz+d4·Vp+d5·Np+d6·AE
[0047] 其中，d0-d6表示此模型的回归系数，Bx,By,Bz为IMF的三分量，Vp,Np分别为太阳风速度和密度，AE为极光电集流指数。
[0048] 进一步地，步骤5中建立的模型按24个磁地方时分别做回归分析，步骤6中用于预测的极光卵边界数据与步骤5中建模所用的边界数据均属于步骤4得到的有效数据集且不重叠；步骤6中采用相同地磁环境下用模型预测的极光卵边界纬度与实际极光卵边界之间的平均绝对误差用来评价模型好坏。
[0049] 与现有技术相比，本发明具有以下有益的技术效果：
[0050] 第一，本发明充分利用Polar卫星返回的海量UVI图像数据，借助计算机
图像处理和
模式识别技术自动获取极光卵边界，为统计分析提供了海量的边界点数据。大
数据库本身可以作为一个查询表，用于对给定的地磁环境来求极光卵的位置。第二，多变量回归分析，首次定量地从根源上来分析极光卵随行星际和太阳风
等离子体变化。第三，回归模型评2
价指标及结果。尽管有少数指标，如确定性系数(R)，F统计量和t统计量可以用来对一个回归模型进行评估，可实际上，对海量的数据量而言，每个系数都很可能是有统计意义的。
为了考察计算模型的实际
质量(而非这些统计指标)，本发明把
训练数据和测试数据分开对待，即用80％的数据点来拟合模型参数，剩余20％当做测试集，用训练好的模型来对其进行预测。相同地磁环境下，用模型估计的边界纬度与实际边界之间的平均绝对误差(MAD)用来评价模型好坏。仅采用最简单的一元线性回归模型，就获得了极向边界平均MAD值大约为1.7个地磁纬度，赤道向边界平均MAD值约为1.8个地磁纬度，这说明本发明提出的方法是非常有效的。本发明对每个磁地方时处的极光卵赤道向&极向边界位置都进行了多元回归建模，这样得到的回归模型可以用来对所有磁地方时处的赤道向/极向边界的地磁纬度进行预测，因而可以用于空间天气预测整个极光卵的位置。
附图说明
[0052] 图2是本发明极光卵边界自动提取示意图；其中(a)为极光原始图像，(b)为预处理后的极光图像，(c)为用SFCM
算法聚类得到的结果，(d)为根据聚类结果把极光分割为2类(即极光卵区域和背景区域)，(e)为对分割结果做进一步平滑处理，(f)为获得的极光卵极向边界(内边界)和赤道向边界(外边界)；
[0053] 图3是本发明的4个模型预测结果与极光卵实际边界的平均绝对误差图，其中(a)为模型1和2的极向边界与赤道向边界预测误差，(b)为模型3和4的极向边界与赤道向边界预测误差；
[0054] 图4是本发明的模型3预测边界与实际极光卵边界、SFCM算法自动分割出的边界三种边界对比图，其中极光卵的真实边界在图中表现为白色阴影区域，SFCM自动分割得出的边界点用白色圆点表示，模型3预测边界用菱形黑点表示；图(a)-(f)所示的极光图像分别来自1997年12月25日的05:34:48，1997年1月2日的11:12:00，1997年12月20日的21:20:39，1997年1月23日的08:43:21，1997年1月13日的11:30:20，1998年1月29日的13:44:18。
具体实施方式
[0055] 下面结合附图对本发明作进一步详细描述：
[0056] 参见图1，一种基于行星际和地磁参数的UVI图像极光卵边界位置建模和预测方法，本发明包括4种基于行星际和地磁参数的UVI极光卵边界模型，通过对UVI极光图像预处理模
块、极光卵边界自动分割模块、行星际及地磁参数获取模块、数据集构建模块、回归建模模块、预测评估模块完成以下具体步骤：
[0057] 步骤1：UVI极光图像预处理
[0058] 在对极光卵边界提取与建模之前，需要对高噪声的Polar卫星UVI极光图像进行预处理操作，消除噪声点，移除无效
背景像素，具体如下：
[0059] 步骤1.1：因为极光发生在高纬地区，移除大范围的背景(如低纬区域)；
[0060] 步骤1.2：掩膜：根据UVI(紫外)极光图像外边界呈椭圆状的特点，用直径为228像素的圆去截取原图，并从中间把宽度限制到200像素，最后图像大小为228*200；
[0061] 步骤1.3：负值点清零：清除因仪器噪声等导致的图像像素值为负值的点；
[0062] 步骤1.4：剔除噪声点：噪声特点是小范围内的大强度像素，所以把那些面积非常小(经验值，如50个像素)、极光强度非常强(大于图像灰度均值+3倍方差)的区域的像素值重置为整幅图像像素的灰度均值；
[0063] 步骤2：极光卵边界自动分割
[0064] 步骤2.1：利用包含空间信息的模糊c均值聚类方法(SFCM)，把极光图像聚成6类(图2(c))；
[0065] 步骤2.2：根据聚类结果，把极光图像分割为前景(极光卵)和背景区域(图2(d))；基于SFCM聚类结果把极光图像分割为前景(极光卵)和背景区域的准则为：强度最小的2个或3个聚类簇被视为背景，其他则为极光卵区。默认是将强度最小的2个聚类簇视为背景区域，但如果此时极光卵区域形状很不规则(与椭圆相比)，则将强度最小的3个聚类簇视为背景区域；
[0066] 步骤2.3：对聚类后出现的一些小块与周围区域类别不同的情况，用其周围邻域的类别信息来对其进行平滑处理(图2(e))；
[0067] 步骤2.4：边界获取并将其映射到磁地方时(MLT)-地磁纬度(MLAT)坐标系下。图2(f)外边界点表示赤道向边界，内边界点表示极向边界。
[0068] 步骤3：行星际、地磁参数获取
[0069] 利用美国国家航空航天局(NASA)提供的OMNI数据库，根据每幅极光图像发生的世界时(UT)时间，可以求得每个极光边界点对应的行星际磁场(IMF)、太阳风(SWD)和极光电集流指数(AE)、极盖指数(PC)的值。同一幅图像中所有边界点的世界时(UT)时间相同，这些参数值相等。
[0070] 步骤4：数据集构建
[0071] 对获得的边界数据结合行星际地磁参数值做进一步筛选。具体地，仅关注北半球区域的边界点、删除那些有一个或多个无效参数值的边界点、对每个参数忽略那些不经常出现的情况。
[0072] 步骤5：回归建模
[0073] 用一元线性回归模型，利用行星际和地磁参数，对赤道向和极向极光卵边界位置进行建模。因变量取赤道向或极向极光卵边界的地磁纬度值，自变量取行星际和地磁参数。
[0074] 模型1：旨在从源头来研究极光卵的变化。首先考虑行星际磁场(IMF)三分量。IMF尤其是IMF Bz分量能控制极光卵的位置已被广泛认识，不同的IMF条件会使得极光卵转向不同的方向。第二，除了IMF外，太阳风(SWD)，特别是太阳风动压(Pdyn)，能显著地影响地球磁场，从而影响极光。尽管很多之前的研究都把太阳风动压作为一个整体来考虑，可太阳风密度(Np)和速度(Vp)与磁层的交互作用是不同的，因此在这个回归模型中，太阳风动压、速度、密度都被考虑进来。即模型参量包括行星际6个单独量(IMF Bx，By，Bz；Pdyn，Vp，Np)。
[0075] 模型1的回归方程为：
[0076] MLAT＝a0+a1·Bx+a2·By+a3·Bz+a4·Pdyn+a5·Vp+a6·Np
[0077] 其中，a0-a6表示此模型的回归系数，Bx,By,Bz为行星际磁场(IMF)的三分量，Pdyn,Vp,Np分别为太阳风动压、速度和密度三分量。
[0078] 模型2：为了达到更好的预测效果，在模型2中，我们除了模型1中包含的6个参量，还加入了三个与极光活动紧密相关的地磁指数。这三个参数分别为：1)极盖(PC)指数，它能监控极盖区的地磁扰动并反映极光卵的大小；2)极光电集流指数(AE)，它与大的地磁活动比如磁暴或亚暴紧密相关，因而能反映夜侧极光卵的边界情况；3)Kan-Lee电场(Ekl)，这是一个太阳风-磁层耦合函数，用来反映日侧重联中从太阳风到磁层的能
力，特别是成分重联。因此，在此模型中，模型参量包括行星际6个单独量(IMF Bx,By,Bz；Pdyn,Vp,Np)以及地磁参数AE、PC和Ekl，
[0079] 模型2的回归方程为：
[0080] MLAT＝b0+b1·Bx+b2·By+b3·Bz+b4·Pdyn+b5·Vp+b6·Np+b7·AE+b8·PC+b9·Ekl[0081] 其中，b0-b9表示此模型的回归系数，Bx,By,Bz为IMF的三分量，Pdyn,Vp,Np分别为太阳风动压、速度和密度三分量，AE为极光电集流指数，PC为极盖指数，Ekl为Kan-Lee电场，它与行星际参数之间的关系为：
[0082] 其中Ekl与基本的行星际参数指数之间的关系为：
[0083] 其中θ取值方法为：
[0084] 若Bz>0,By>＝0，θ＝arctan(|By|/Bz)；
[0085] 若Bz<＝0,By>0，θ＝arctan(|Bz|/By+π/2)；
[0086] 若Bz<0,By<＝0，θ＝arctan(|By|/|Bz|+π)；
[0087] 若Bz>＝0,By<0，θ＝arctan(Bz/|By|+1.5π)；
[0088] 若Bz＝0,By＝0，θ＝0。
[0089] 模型3：极光卵分布由行星际参数和地磁环境决定，但这些参数是通过许多不同的磁层动力学过程影响磁层的。所以考虑用动力学参数来建模。模型参数包括AE指数、2
Kan-Lee电场(Ekl)、IMF磁能通量(VB)、太阳风动能通量(VPdyn)、南北电场(Ez)，[0090] 模型3的回归方程为：
[0091] MLAT＝c0+c1·VB2+c2·VPdyn+c3·Ekl+c4·Ez+c5·AE
[0092] 其中，c0-c5表示此模型的回归系数，VB2为IMF磁能通量，VPdyn为太阳风动能通2
量，Ekl为Kan-Lee电场，Ez为南北电场，AE为极光电集流指数；VB 、VPdyn、Ez与行星际参数之间的关系为：
[0093] 除AE、Ekl外，另外3个动力学参数与基本的行星际参数之间的关系为：
[0094]
[0095] VPdyn＝Vp·Pdyn
[0096] Ez＝By·Vp
[0097] 模型4：独立参量模型。不考虑参量之间的耦合作用，选择完全独立的6个影响极光演变的参量，包括IMF Bx、By、Bz，太阳风速度(Vp)、密度(Np)及AE指数。
[0098] 模型4的回归方程为：
[0099] MLAT＝d0+d1·Bx+d2·By+d3·Bz+d4·Vp+d5·Np+d6·AE
[0100] 其中，d0-d6表示此模型的回归系数，Bx,By,Bz为IMF的三分量，Vp,Np分别为太阳风速度和密度，AE为极光电集流指数。
[0101] 对上述四个模型来说，因为在不同的磁地方时(共24小时)，每一个参数与极光卵边界的相互关系都是不同的，因此建模的时候按24个MLT分别做回归分析。
[0102] 步骤6：极光卵边界位置预测
[0103] 根据步骤5获得极光卵边界模型，对任意给定的行星际和地磁条件，即可计算出相应的赤道向和极向极光卵边界位置。
[0104] 步骤7：预测结果评价
[0105] 极光卵边界位置预测包含定性和定量两方面的评价。其中定性评价为：对比UVI图像极光卵的真实边界、SFCM
聚类算法自动分割得到的边界点以及运用回归模型得到的预测边界，视觉比较三种边界，评价算法及模型的效果；定量评价为：根据行星际和地磁参数用回归模型估计极光卵边界位置，计算其与真实边界的平均绝对误差值(MAD)，给出各模型的效果。
[0106] 下面结合具体
实施例对本发明做进一步描述：
[0107] 本发明基于行星际和地磁参数，提出4种的UVI图像极光卵边界位置的回归模型，用来对极光卵边界位置进行预测。为了使实验结果更可信，本发明将利用Polar卫星拍摄到的1996年12月～1999年1月三个冬季，共6万多幅UVI极光图像构成的数据库来研究极光卵边界在不同地磁条件下的变化；且所有的空间参数都采用他们的真实值来做回归拟合。
[0108] 步骤一：对高噪声的Polar UVI极光图像进行预处理操作(附图2(b)所示)。具体操作包括：1)移除大范围的背景(如低纬区域)；2)掩膜：根据UVI极光图像外边界呈椭圆状的特点，用直径为228像素的圆去截取原图，并从中间把宽度限制到200像素，最后图像大小为228*200；3)负值点清零：清除因仪器噪声等导致的图像像素值为负值的点；4)剔除噪声点：噪声特点是小范围内的大强度像素，所以把那些面积非常小(经验值，如50)、极光强度非常强(大于图像均值+3倍方差)的区域像素值重置为其周围邻域图像像素的均值。
[0109] 步骤二：运用SFCM聚类算法自动提取极光卵边界。在分析了大量UVI图像后，我们最终决定将每幅图像聚成6类(附图2(c))。根据这些聚类结果，强度最小的两个或三个聚类簇被视为背景，其他则为极光卵区(附图2(d))。具体地，默认将强度最小的2个聚类簇视为背景区域，但如果此时极光卵区域形状很不规则(与椭圆相比)，则将强度最小的3个聚类簇视为背景区域。因为聚类过程是基于像素强度进行的，所以经常会出现图像中一些小块与周围区域类别不同的情况，这不符合实际情况，所以我们用其周围邻域的类别信息来对这些小块进行平滑处理(附图2(e))。附图2(f)给出了获得的赤道向边界(外边界)和极向边界(内边界)结果。为了能更好的进行后续回归分析，我们对一些不好的分割结果进行了人工剔除，如一些“θ极光”。
[0110] 步骤三：行星际和地磁参数获取。根据每幅图像对应的UT时间，我们都可以由NASA OMNI数据获得每个边界点的MLT、MLAT值，以及其对应的行星际和地磁环境(包括行星际磁场IMF三分量Bx,By,Bz值、太阳风三参量Pdyn,Vp,Np值、及AE、PC参数值)。据此，我们就可以把极光卵边界映射到[MLT，MLAT]地磁坐标系下来讨论。而且根据这8个基本参量，我们还可以求其他如Ekl、Ez、VPdyn等动力学参数值。
[0111] 步骤四：数据集构建。经过上述处理，我们共得到了极向边界分割很好的61210幅图像和赤道向边界分割很好的60180幅图像。为了研究不同参数对极光卵边界的影响，获得的边界点被进行了进一步的处理。首先，我们仅仅关注北半球区域的边界点。其次，删除那些有一个或多个无效参数值的边界点。第三，对每个参数，包括地磁纬度，我们忽略那些不经常出现的情况，即，那些值特别大或特别小的数据(在20个bin的直方图里，那些出现
频率小于0.02的bin里的数据视为不经常出现的数据)，因为这些点不是统计上显著的。第四，我们认为5min间隔内极光卵边界位置变化不会超过5个地磁纬度。这样处理之后，超过380.5万个极向边界点和121.5万个赤道向边界点被提取出来。
[0112] 步骤五：回归建模。根据上面得到的数据集，我们可以看出，在不同的磁地方时，任何一个参数与极光卵边界的相互关系是不同的。因此我们需要对每个磁地方时分别做回归分析。第二，极光卵的活动是一个综合反映，没有哪个空间参数可以单独将其完全表征。所以，需要多元回归分析。因为目前导致极光卵变化的具体原因尚不清楚，我们不知道回归分析的具体形式，所以为了简单起见，本发明中选用最简单、常用的一元线性回归分析方法。回归模型的因变量取赤道向或极向极光卵边界的地磁纬度值，自变量取行星际和地磁参数。本发明选用了4种回归模型，每一个模型的回归参量(自变量)都不同。对每一种模型，建模的时候都是按24个MLT分别做回归分析。
[0113] 步骤六：边界预测与评价。(1)定量评价：为了考察计算模型的实际质量(而非这些统计指标)，我们把训练数据和测试数据分开对待，即用80％的数据点来拟合模型参数，剩余20％当做测试集，用训练好的模型来对其进行预测。预测值与实际值之间的平均绝对误差(MAD)用来评价模型好坏。这个过程在整个数据库上重复操作五次，每次都用不同的20％数据部分作为测试集。这就是所谓的五重交叉验证方法。为了增加实验结果的可靠性，我们进行了5轮上述操作，让每次操作的训练集和测试集数据都不一样。这样得到的25个MAD值取平均后即为最终结果(附图3)。MAD值越小越好。可以看出，最好的是回归模型2的结果，极向边界最佳预测出现在15MLT处，误差为1.31MLAT，平均误差为1.69MlAT，而对赤道向边界，最佳预测出现在24MLT处，为1.29MLAT，平均误差为1.77MLAT。
[0114] (2)定性评价：除了上述定量评价之外，附图4还给出了极坐标系下6幅UVI图像，读者可以从视觉上比较极光卵的真实边界(白色阴影区域)、SFCM自动分割得到的边界点(白色圆点)及运用回归模型3得到的预测边界(菱形黑点)三者之间的差异，由此定性评价本发明的方法效果。其中回归模型3的回归方程为：
[0115] MLAT＝c0+c1·VB2+c2·VPdyn+c3·Ekl+c4·Ez+c5·AE
[0116] 回归建模后得到的极向和赤道向边界如下：(其中VB2、VPdyn、Ekl、Ez由于参数的单位原因，回归建模前参数值除以了100)，详见表1和表2。
[0117] 表1 0100-2400MLT时刻回归模型3的极向边界回归模型系数
[0118]
[0119] 表2 0100-2400MLT时刻回归模型3的赤道向边界回归模型系数
[0120]
[0121]
