CN105844662A - 基于流体力学的极光运动方向确定方法 - Google Patents

Abstract

本发明公开了一种基于流体力学的极光运动方向确定方法，创新性的在极光动态过程的研究中引入了流体力学的理论，主要解决现有极光运动方向自动确定方法效率低的问题。其实现步骤为：1.输入待分析的极光序列；2.对序列中的每帧全天空极光图像进行预处理；3.应用流体力学中的离散格子玻尔兹曼方法计算整个极光动态过程中沿各方向移动的粒子数目；4.根据沿各方向移动的粒子数目，确定极光运动方向。本发明通过结合流体粒子运动，对极光运动方向进行确定。克服了极光作为非刚体的这一特性，具有计算简单、运算量小、处理速度快的优点，提高了确定极光运动方向的效率，可用于极光动态过程分析。

Description

基于流体力学的极光运动方向确定方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种极光运动方向的计算机确定方法，可用于非刚体物体动态过程分析。

背景技术

极光是太阳风通过日侧极隙区注入到地球磁层时，沉降粒子沿磁力线与地球高层大气相互作用而产生的绚丽光辉。从物理的角度来讲，极光是太阳中的高能带电粒子在地球磁场作用下，与极区高层大气中的原子和分子相互作用而产生，也就是说，太阳风、地球磁场和极区高层大气是形成极光的必要。因此，极光的发生反映了太阳与地磁活动之间的动力学关系，有助于人们了解太阳对地球的影响方式和程度。其次，发生极光时辐射出的某些无线电波，会直接影响地球上的无线电通信、导航、定位以及线路传输等。然而极光发生时，在地球大气层中所爆发的能量，几乎可以达到全世界各国发电厂所产生的电容量总和。因此，如何使用极光所产生的巨大能量造福人类，成为当今科学领域的重要研究课题。研究和事实表明，极光现象是太阳系中有磁星体所具有的普遍现象。哈勃太空望远镜已经能很清楚的看到木星和土星这两颗行星上的极光。因此对地球上极光的研究将有助于人类研究其他行星上的极光现象。

极光的动态过程研究属于空间物理学和天文学的范畴，能集中反应外磁层高能等离子体的变化过程，对极光的活动情况进行长期系统地观测，可以帮助人们进一步探究磁层边界和磁层结构的相关动力学方程以及太阳风与地球磁层之间的相互作用关系。

中国北极黄河站的全天空成像系统(All-sky Camera)同时对极光的三个典型谱段427.8nm，557.7nm和630.0nm进行连续观测，产生数以万计的极光序列，数据量巨大。但是，目前关于极光的动态过程分析大多是通过人工标记对少数极光序列进行案例分析，由于所能涉及的极光事件非常有限，案例分析的结果难以推广。2006年，Blixt等人将鲁棒光流场用于极光的全局二维运动场估计，并且研究了光流场方法用于极光数据的可行性。但该方法是基于亮度不变假设的，然而极光的形态在动态过程中会不断地发生变化，因此，假设极光的动态过程满足亮度不变性，这显然是不严谨的。2012年王倩等人在“日侧全天空极光图像分类及动态过程分析方法研究”中提出基于局部运动向量时空统计量的极光事件监测方法。该方法虽说克服了光流场亮度不变假设的约束，但却存在着计算复杂，运算量大，处理速度慢的不足。

发明内容

本发明的目的在于针对上述现有技术存在的不足，引入流体力学中的离散格子玻尔兹曼方法，提出一种基于流体力学的极光运动方向确定方法，以提高极光运动方向自动确定的速度，取代繁杂的人工标记。

为实现上述目的本发明的技术方案包括如下：

(1)输入待分析的极光序列；

(2)对极光序列中的每帧全天空极光图像进行预处理：

(2a)采用半径为220个像素的圆形掩膜去掉全天空极光图像圆周上的边缘噪声；

(2b)使用模糊C-均值聚类方法，对掩膜处理后的全天空极光图像进行粗分割，得到全天空极光图像的弧状区域；

(2c)在极光弧区域内，计算所有像素的平均值，并将极光弧区域内所有像素的灰度值设为该平均值，得到预处理后的全天空极光图像；

(3)采用流体力学理论中的离散格子玻尔兹曼方法进行极光运动方向的确定：

(3a)按照离散格子玻尔兹曼方法中的二维九方向正方形网格元胞结构，把经过预处理的全天空极光图像上的每一个像素点映射为流体粒子，并以相邻全天空极光图像间对应位置上像素点的灰度值差作为外力项，每个粒子在外力的影响下向其周围八个方向的元胞移动或者静止不动，并与移动方向上的元胞发生碰撞；

(3b)根据有外力项的离散格子玻尔兹曼方程，逐帧计算沿各个方向移动的粒子数目；

(3c)判断外力项是否接近于0，即相邻的两帧全天空极光图像是否接近一致，如果外力项接近0，则极光动态过程结束，否则返回(3b)；

(3d)计算整个极光过程中沿各个方向移动的粒子数目之和，以粒子移动方向标号为横坐标，沿该方向移动的粒子数目之和为纵坐标，绘制这些移动粒子数目在各方向上的分布直方图，该直方图峰值所处的方向即为极光运动方向。

本发明与现有技术相比具有以下优点：

1.本发明引入流体力学的思想进行极光运动方向分析，克服了极光作为非刚体的这一特性；

2.取代以往人工分析的手段，完全使用计算机自主对输入序列进行分析；

3.计算简单、运算量小、具有很高的可并行性。

附图说明

图1是本发明的实现流程图；

图2是本发明实验中待处理极光序列中的3幅极光图像；

图3是本发明实验中对待处理极光序列中全天空极光图像经过掩膜处理后的结果；

图4是本发明实验中在掩膜处理后的全天空极光图像中提取出的极光弧区域；

图5是本发明实验中经全部预处理后的全天空极光图像；

图6是离散格子玻尔兹曼模型的元胞空间示意图；

图7是本发明设定的运动方向的示意图；

图8是地磁南北极标示图；

图9是本发明实验得到的沿不同方向移动的粒子数目的统计直方图。

具体实施方式

下面结合附图对本发明做进一步的详细描述。

步骤1，从全天空极光图像数据库中选取一个有明显运动的弧状极光序列作为待分析极光序列，如图2所示。

步骤2，对序列中的每幅全天空极光图像进行预处理。

(2a)构造一幅与输入全天空极光图像尺寸一致的数值矩阵B，将B中以矩阵中心为圆心，半径为220的圆内的点的值设为1，其余点的值设为0，利用如下公式对输入的极光图像进行掩膜处理，得到掩膜处理后的全天空极光图像：O₁＝I*B，如图3所示，其中I为原始输入的全天空极光图像；

(2b)使用模糊C-均值聚类方法，对掩膜处理后的全天空极光图像进行粗分割，得到全天空极光图像的弧状区域：

(2b1)将掩膜处理后的全天空极光图像O₁转换为数值矩阵U，输入到模糊C-均值聚类算法中，并设定聚类数为3；

(2b2)运行模糊C-均值聚类算法，得到与U尺寸一致的数值矩阵V，V中元素的取值范围为：{1,2,3}；

(2b3)将V中的1和3全部替换为0，得到极光弧区域分割图，如图4所示；

(2c)生成预处理后的全天空极光图像：

(2c1)构造与V尺寸相同的数值矩阵V₁：

V₁＝O₁*V，

(2c2)计算V₁中所有非零数的平均值m；

(2c3)生成预处理后的全天空极光图像O，如图5所示：

O＝M*V+O₁*(E-V)，

其中，M为与V尺寸相同且值全部为m的数值矩阵，E为与V尺寸相同且值全部为1的数值矩阵。

步骤3，根据步骤1中得到的预处理后的全天空极光图像O，采用流体力学理论中的离散格子玻尔兹曼方法计算沿每个方向移动的粒子密度分布函数：

(3a)按照离散格子玻尔兹曼方法中的二维九方向正方形网格元胞结构，把经过预处理的全天空极光图像上的每一个像素点映射为流体粒子，并构造元胞空间；

(3b)令为t时刻全天空极光图像上处沿方向i移动的粒子的密度分布函数，为t时刻全天空极光图像上处粒子平衡分布函数，为t时刻全天空极光图像上处的粒子密度，即：

$\mathrm{\ρ}(\stackrel{\→}{r},t)=\underset{i=0}{\overset{8}{\Σ}}{f}_i(\stackrel{\→}{r},t),$

将粒子密度分布函数和平衡分布函数的初始值设置为：

$f_i(\stackrel{\→}{r},0)=f_i^{\left(0\right)}(\stackrel{\→}{r},0)=\left\{\begin{array}{c}\frac{M\mathrm{\ρ}(\stackrel{\→}{r},0)}{9},i\≠0\\ 1-\frac{8M\mathrm{\ρ}(\stackrel{\→}{r},0)}{9},i=0\end{array}\right.,$

其中，M∈[0,1]；

(3c)令F为相邻全天空极光图像间对应位置上像素点的灰度值差，根据下式计算外力项g_i(x)：

$g_i\left(x\right)=\left\{\begin{array}{c}0,i=0\\ \frac{C_i}{3}\·F,i=1,2,3,4\\ \frac{C_i}{12}\·F,i=5,6,7,8\end{array}\right.,$

其中：

$C_i=\left\{\begin{array}{c}(0,0),i=0\\ (\cos {\mathrm{\θ}}_i,sin{\mathrm{\θ}}_i),{\mathrm{\θ}}_i=\frac{(i-1)\mathrm{\π}}{2},i=1,2,3,4\\ \sqrt{2}(\cos {\mathrm{\θ}}_i,sin{\mathrm{\θ}}_i),{\mathrm{\θ}}_i=\frac{(i-5)\mathrm{\π}}{2},i=5,6,7,8\end{array}\right.;$

(3d)根据有外力项的离散格子玻尔兹曼方程，计算各时刻沿每个方向移动的粒子密度分布函数：

$f_i(\stackrel{\→}{r}+\mathrm{\Δ}h\·C_i,t+\mathrm{\Δ}t)=f_i(\stackrel{\→}{r},t)+\frac{1}{\mathrm{\ζ}}\[f_i^{\left(0\right)}(\stackrel{\→}{r},t)-f_i(\stackrel{\→}{r},t)\]+g_i\left(x\right),$

其中，Δt和Δh分别为时间步长和空间步长；ζ为松弛时间，即元胞内粒子密度趋向于平衡态的时间；

(3e)计算t时刻沿每个方向运动的粒子数目：

$N_i\left(t\right)=\∫f_i(\stackrel{\→}{r},t)dS.$

步骤4，确定极光运动方向：

(4a)计算整个极光过程中沿各个方向移动的粒子数目之和：

N_i＝∑N_i(t)；

(4b)以粒子移动方向标号为横坐标，沿该方向移动的粒子数目之和为纵坐标，绘制这些移动粒子数目在各方向上的分布直方图，该直方图峰值所处的方向即为极光运动方向。

本发明的效果通过以下实验进一步说明：

1、实验条件

实验硬件设备：Intel Core i3(2.93GHz)、4GB RAM

实验软件平台：Windows7 64bit操作系统下的MATLAB R2012a；

实验数据：对我国北极黄河站的全天空极光数据中2003年12月到2004年G波段极光图像进行手工标记，去除天气等因素造成的无效数据，标记出102个弧状极光序列，序列长度为15帧到35帧之间。

2、实验内容与结果

选取2003.12.25的12:39-12:54中连续的24幅全天空极光图像作为图像序列，确定其运动方向，其中的三幅全天空极光图像如图2所示，对这三幅图像经过掩膜处理后的结果如图3所示，使用模糊C-均值聚类方法在掩膜处理后的全天空极光图像中提取出的极光弧区域如图4所示，经全部预处理后的全天空极光图像如图5所示

对比图7中本发明设定的运动方向的示意图和图8的地磁南北极标示，粒子运动应该主要集中在运动方向1区域，其次集中在5区域。即各运动方向粒子数的统计直方图应主要在1方向上存在峰值，其次在5方向上存在较高的值，其余方向上数值较低。

根据步骤3计算沿每个方向移动的粒子密度分布函数，计算整个极光过程中沿各个方向移动的粒子数目之和，以粒子移动方向标号为横坐标，沿该方向移动的粒子数目之和为纵坐标，绘制这些移动粒子数目在各方向上的分布直方图得到图9。

图9中的分布直方图主要在1方向上存在峰值，其次在5方向上存在较高的值，其余方向上数值较低，其结果与理论预测的结果一致，说明通过该方法成功实现了极光运动方向的自动确定。

Claims (3)

1.基于流体力学的极光运动方向确定方法，包括：

(1)输入待分析的极光序列；

(2)对极光序列中的每帧全天空极光图像进行预处理：

(2a)采用半径为220个像素的圆形掩膜去掉全天空极光图像圆周上的边缘噪声；

(2b)使用模糊C-均值聚类方法，对掩膜处理后的图像进行粗分割，得到全天空极光图像的弧状区域；

(2c)在极光弧区域内，计算所有像素的平均值，并将极光弧区域内所有像素的灰度值设为该平均值，得到预处理后的全天空极光图像；

(3)采用流体力学理论中的离散格子玻尔兹曼方法进行极光运动方向分析：

(3a)按照离散格子玻尔兹曼方法中的二维九方向正方形网格元胞结构，把经过预处理的全天空极光图像上的每一个像素点映射为流体粒子，并以相邻全天空极光图像间对应位置上像素点的灰度值差作为外力项，每个粒子在外力的影响下向其周围八个方向的元胞移动或者静止不动，并与移动方向上的元胞发生碰撞；

(3b)根据有外力项的离散格子玻尔兹曼方程，逐帧计算沿各个方向移动的粒子数目；

(3c)判断外力项是否接近于0，即相邻的两帧全天空极光图像是否接近一致，如果外力项接近0，则极光动态过程结束，否则返回(3b)；

(4)计算整个极光过程中沿各个方向移动的粒子数目之和，以粒子移动方向标号为横坐标，沿该方向移动的粒子数目之和为纵坐标，绘制这些移动粒子数目在各方向上的分布直方图，该直方图峰值所处的方向即为极光运动方向。

2.根据权利要求1所述的方法，其中步骤(2b)中全天空极光图像弧状区域的获取，按如下步骤进行：

(2b1)设定聚类数为3及迭代停止阈值ε，构造隶属矩阵U＝(u_ij)_2*N，用值在0和1之间的随机数初始化该隶属度矩阵，其中N为单幅全天空极光图像上的像素点数目，u_ij表示第i个类型下第j个像素点的隶属度；

(2b2)根据下式计算各聚类中心c_i,i＝1,2，3：

$c_i=\frac{\underset{j=1}{\overset{N}{\Σ}}{\left(u_{ij}\right)}^m{x}_j}{\underset{j=1}{\overset{N}{\Σ}}{\left(u_{ij}\right)}^m}$

其中x_j为第j个像素点的位置，m∈[1,∞)是一个加权指数，取值为2；

(2b3)根据下式计算目标函数：

$J(U,c_1,c_2,c_3)=\underset{i=1}{\overset{3}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(u_{ij}\right)}^m{d_{ij}}^2,$

其中，d_ij＝||c_i-x_j||为第i个聚类中心与第j个数据点间的欧几里何距离；

(2b4)将目标函数值与停止阈值ε进行比较，如果目标函数值小于停止阈值ε，则此时的隶属度矩阵为最优的隶属度矩阵，否则更新第i个类型下第j个像素点的隶属度：然后返回(2b2)对隶属矩阵U进行更新；

(2b5)使用最优隶属度矩阵实现去模糊化，令C_j表示第j个数据点所属的类别：用所有C_j＝2的数据点构成的区域即为极光弧区域，其中分别表示‘所有i’和‘所有j’。

3.根据权利要求1所述的方法，其中步骤(3b)中计算沿各个方向移动的粒子数目，按如下步骤进行：

(3b1)令为t时刻全天空极光图像上处沿方向i移动的粒子的密度分布函数，为t时刻全天空极光图像上处粒子平衡分布函数，为t时刻全天空极光图像上处的粒子密度，即：

$\mathrm{\ρ}(\stackrel{\→}{r},t)=\underset{i=0}{\overset{8}{\Σ}}{f}_i(\stackrel{\→}{r},t),$

将粒子密度分布函数和平衡分布函数的初始值设置为：

$f_i(\stackrel{\→}{r},0)={f_i}^{\left(0\right)}(\stackrel{\→}{r},0)=\left\{\begin{array}{c}\frac{M\mathrm{\ρ}(\stackrel{\→}{r},0)}{9},i\≠0\\ 1-\frac{8M\mathrm{\ρ}(\stackrel{\→}{r},0)}{9},i=0\end{array}\right.$

其中，M∈[0,1]；

(3b2)令F为相邻全天空极光图像间对应位置上像素点的灰度值差，根据下式计算外力项g_i(x)：

$g_i\left(x\right)=\left\{\begin{array}{c}0,i=0\\ \frac{C_i}{3}\·F,i=1,2,3,4\\ \frac{C_i}{12}\·F,i=5,6,7,8\end{array}\right.$

其中，

$C_i=\left\{\begin{array}{c}(0,0),i=0\\ ({\mathrm{cos\θ}}_i,{\mathrm{sin\θ}}_i),{\mathrm{\θ}}_i=\frac{(i-1)\mathrm{\π}}{2},i=1,2,3,4\\ \sqrt{2}({\mathrm{cos\θ}}_i,{\mathrm{sin\θ}}_i),{\mathrm{\θ}}_i=\frac{(i-5)\mathrm{\π}}{2},i=5,6,7,8\end{array}\right.;$

(3b3)根据有外力项的离散格子玻尔兹曼方程，计算各时刻沿每个方向移动的粒子密度分布函数：

$f_i(\stackrel{\→}{r}+\mathrm{\Δ}h\·C_i,t+\mathrm{\Δ}t)=f_i(\stackrel{\→}{r},t)+\frac{1}{\mathrm{\ζ}}\[{f_i}^{\left(0\right)}(\stackrel{\→}{r},t)-f_i(\stackrel{\→}{r},t)\]+g_i\left(x\right),$

其中，Δt和Δh分别为时间步长和空间步长；ζ为松弛时间，即元胞内粒子密度趋向于平衡态的时间；

(3b4)计算t时刻沿每个方向运动的粒子数目：

$N_i\left(t\right)=\∫f_i(\stackrel{\→}{r},t)dS.$