一种电离层层析方法及装置与流程

本发明涉及一种电离层层析方法及装置，属于大地测量及大气探测技术领域。

背景技术：

电离层是从大约50到1000公里延伸到等离子体层的部分电离的大气高层，电离层中含有足够的自由电子来影响无线电波传播。电离层电子密度分布受到太阳辐射、粒子沉淀等诸多自然因素的影响，导致自由电子在时间和空间上的电子分布都有很大的不同。因此，从卫星应用以及了解电离层动力学的角度，电离层电子密度分布的实时确定尤为重要。

电离层层析是利用实测电离层的电子含量反演所在路径上的电子密度，在电离层层析过程中，考虑到实际的情况，将反演的电离层空间离散为一个个的立体格网，称为像素，以像素中心点的电子密度值代表该像素的电子密度。这样射线穿过每个像素的电子含量值可以表示为射线在该像素内的截距与电子密度的乘积，该路径上的总电子含量为各个像素电子含量的总和。

mart算法(传统乘法代数重构算法)作为一种常用的电离层层析算法，目前被广泛应用于电离层电子密度的反演过程中，例如申请公布号为cn104933737a的中国发明专利申请文件公开了一种基于共轭梯度法的电离层层析成像混合反演方法，主要通过迭代方式调整电子密度。

如图1所示，现有mart算法利用iri2012经验模型获取电子密度初值，然后采用迭代的方式逐步改善待重构电离层中的电子密度值，其中mart算法的迭代方程为：

其中，为第j个像素中电子密度的第k+1步迭代值；为第j个像素中电子密度的第k步迭代值；yi为第i条射线的实测电子含量；x^k为第k步迭代时第i条射线上各个像素的电子密度组成的向量；ai为第i条射线上各个像素的截距值组成的向量，aij为第i条射线穿过第j个像素的截距值；u为迭代系数；λk为松弛因子，一般取0<λk<2；n为一条射线包含的像素数目；为迭代方程中的组合量。

从图1可以得到，当电子含量的实测值与反演值(<x^k,ai>)的误差在设定范围时则迭代结束，否则对电子密度继续进行迭代。

可以看出，mart算法中将电子含量的实测值与反演值的误差以比值的形式给出，并以此比值作为迭代方程的组合量，但是这样就会使得组合量的形式较为单一，对影响迭代的因素考虑不全面，导致mart算法的迭代精度不高。

技术实现要素：

本发明的目的是提供一种电离层层析方法及装置，用于解决现有mart算法中构建组合量时对影响因素的考虑不全面导致mart算法反演结果存在误差的问题。

为解决上述技术问题，本发明提供了一种电离层层析方法，包括以下步骤：

根据mart算法中的组合量和截距值构建新的组合量，以新的组合量确定电离层电子密度的迭代方程，根据该迭代方程对电离层进行层析反演得到电子密度；所述迭代方程为：

其中，为第i条射线上第j个像素中电子密度的第k+1步迭代值；为第i条射线上第j个像素中电子密度的第k步迭代值；yi为第i条射线的实测电子含量；x^i,k为第k步迭代时第i条射线上各个像素的电子密度组成的行向量；ai为第i条射线上各个像素的截距值组成的列向量；<x^i,k,ai>表示求行向量x^i,k和列向量ai的乘积；aij为第i条射线穿过第j个像素的截距值；u为迭代系数；为第i条射线对应的松弛因子；n为一条射线包含的像素数目；为新的组合量。

为解决上述技术问题，本发明还提供了一种电离层层析装置，包括处理器和存储器，所述处理器用于处理存储在所述存储器中的指令以实现如下方法：

根据mart算法中的组合量和截距值构建新的组合量，以新的组合量确定电离层电子密度的迭代方程，根据该迭代方程对电离层进行层析反演得到电子密度；所述迭代方程为：

其中，为第i条射线上第j个像素中电子密度的第k+1步迭代值；为第i条射线上第j个像素中电子密度的第k步迭代值；yi为第i条射线的实测电子含量；x^i,k为第k步迭代时第i条射线上各个像素的电子密度组成的行向量；ai为第i条射线上各个像素的截距值组成的列向量；<x^i,k,ai>表示求行向量x^i,k和列向量ai的乘积；aij为第i条射线穿过第j个像素的截距值；u为迭代系数；为第i条射线对应的松弛因子；n为一条射线包含的像素数目；为新的组合量。

本发明的有益效果是：将通过像素内的截距值与电子含量作为组合量，构造新的组合量，由于综合考虑了电子含量和截距值，而截距值直接影响着迭代精度，有效提高了反演精度，使迭代结果更加准确。

作为方法和装置的进一步改进，为了削弱噪声对反演结果的影响，克服由测量误差导致的不稳定性问题，所述松弛因子为：

其中，为第i条射线上第j+1个像素中电子密度的第k步迭代值；m为迭代次数。

附图说明

图1是现有技术mart算法的流程图；

图2是本发明电离层层析方法的流程图；

图3(a)是现有技术mart算法和iri2012经验模型重构的电子密度对比图；

图3(b)是本发明imart算法和iri2012经验模型重构的电子密度对比图；

图4(a)-图4(d)为本发明imart算法重构的电子密度在经度面111.4°e上的周日变化图；

图5(a)为本发明2016年8月9日在00:30ut时刻的电子密度轮廓对比图；

图5(b)为本发明2016年8月9日12:30ut时刻的电子密度轮廓对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施例对本发明进行进一步详细说明。

电离层层析方法实施例：

本实施例提供了一种电离层层析方法(简称为imart算法)，其对应的流程图如图2所示，具体包括以下步骤：

(1)获取观测数据：确定目标区域范围，选取目标区域范围内卫星导航系统观测数据和电离层测高仪数据。

也就是，在步骤(1)中，需要确定重构区域的经纬度范围、高度范围和观测时段，获取反演区域的全球导航卫星系统(gnss)观测数据和电离层测高仪数据。其中，在本实施例中，重构区域的经纬度范围分别为108.5°-114.5°e和24.5°-30.5°n，高度范围为70-1000km，电离层测高仪数据为(111.4°e,27.2°n)的实测数据。使用2016年8月9日四个时段(04:30-05:00ut，08:00-08:30ut，16:00-16:30ut，20:30-21:00ut)的观测数据，采用经度0.5°×纬度0.5°×高度15km的像素尺度对电离层区域进行剖分。

(2)使用基于像素基的电离层层析模型，将待反演的电离层空间离散化，该基于像素基的电离层层析模型的数学表达式如下：

其中，stec为电子总含量，ne是射线ρ上某一位置的电子密度，是该位置在时间t的位置向量。

为了获得电子密度，电离层被剖分成大量的小像素，且假设每个像素内的电子密度ied在一个短时期内是一个常数，因此每个像素中的各条射线对应的电子密度是一样的。此时，公式(1)中沿着第i条射线的电子总含量steci可近似表示为：

其中，aij为第i条射线穿过第j个像素的截距值，xj为射线穿过第j个像素对应的电子密度，εi为误差值，n为一条射线包含的像素数目。

公式(2)可以用简单的矩阵符号写成：

ym×1＝am×nxn×1+εm×1(3)

其中，n为一条射线包含的像素数目，m为射线总数，ym×1为由m个已知的电子总含量stec值组成的列向量，am×n为m条射线的截距组成的矩阵，xn×1为所有像素内电子密度组成的向量，εm×1是由离散误差和仪器误差组成的向量。

(3)计算电子含量的实测值。

具体的，步骤(3)利用全球导航卫星系统双频观测数据，使用如下的公式求解电子总含量stec，该公式如下：

其中，p4,sm为使用载波相位平滑后的伪距观测值，f1和f2为全球导航卫星系统双频观测数据的频率，bcdⁱ为全球导航卫星系统卫星的仪器偏差，dcbj为全球导航卫星系统接收机的仪器偏差。

利用全球导航卫星系统双频观测数据，通过平滑，可以获取更高精度的dcb值，在此基础上，可获得可靠的电子含量的实测值。由于根据全球导航卫星系统双频观测数据计算电子含量是本领域的现有技术，在此不再展开描述。

(4)根据mart算法中的组合量和截距值构建新的组合量，以新的组合量确定电离层电子密度的迭代方程，根据该迭代方程对电离层进行层析反演得到电子密度；迭代方程为：

其中，为第i条射线上第j个像素中电子密度的第k+1步迭代值；为第i条射线上第j个像素中电子密度的第k步迭代值；yi为第i条射线的实测电子含量；x^i,k为第k步迭代时第i条射线上各个像素的电子密度组成的行向量；ai为第i条射线上各个像素的截距值组成的列向量；<x^i,k,ai>表示求行向量x^i,k和列向量ai的乘积；aij为第i条射线穿过第j个像素的截距值；u为迭代系数；为第i条射线对应的松弛因子；n为一条射线包含的像素数目；为新的组合量。

松弛因子在层析算法中起着重要作用，它有助于克服由测量误差导致的不稳定性问题。当使用的数据中含有噪声时，即使降低收敛速度，通过适当选择松弛因子也可以改善重建质量。松弛因子的取值对反演结果的影响也较为明显，当松弛因子取较大值时，则反演电子密度较为平滑，电离层局部变化特征将会被掩盖；当松弛因子取较小值时，又会受到噪声的干扰。

因此，为了有效克服传播噪声对电子密度反演的影响，使松弛因子自适应地和电子密度相关联，在本实施例中，松弛因子的表达式为：

其中，为第i条射线上第j个像素中电子密度的第k步迭代值；为第i条射线上第j+1个像素中电子密度的第k步迭代值；m为迭代次数，其取值视实际情况而定。

需要说明的是，在忽略松弛因子对电子密度的影响的情况下，作为其他的实施方式，上述步骤(4)中迭代方程中的松弛因子也可以使用现有技术中的取值，例如可以取为0～2。

(5)解算反演区域电离层电子密度。

其中，由于步骤(5)中解算反演区域电离层电子密度的具体过程属于现有技术，此处不再赘述。

为了验证上述的电离层层析方法的可靠性和优越性，针对步骤(1)中确定的目标区域范围以及选取的目标区域范围内卫星导航系统观测数据和电离层测高仪数据，分别运行现有技术中的iri2012经验模型和mart算法以及上述的imart算法，所获取的电子总含量stec重构误差结果如表1所示。表1中包括四个时段内stec的均方根误差rmsereal和绝对误差ereal，其误差计算方法如下：

其中，m为射线总数，是第i条射线上由反演电子密度计算得到的stec，是第i条射线上由实测全球导航卫星系统观测数据提取的stec。

表1

由表1可知，imart算法相较于iri2012经验模型和mart算法，其反演性能更加突出。另外，与mart算法和imart算法相比，iri2012经验模型的性能最差，例如在04:30-05:00ut时段，mart算法、imart算法的均方根误差rmsereal分别为1.7622、1.5493，其绝对误差ereal分别为1.0729、1.0574，然而由iri2012经验模型得到的rmsereal和绝对误差ereal却非常的高。总的来说，imart算法、mart算法和iri2012经验模型获取的平均rmsereal分别为1.5024、1.8321和4.2130，由此可知，imart算法比mart算法的反演精度提高了18.0％、比iri2012经验模型提高了64.34％。

图3(a)为2016年8月9日04:30-05:00ut历元时刻mart算法和iri2012经验模型重构的电子密度比较情况，图3(b)为2016年8月9日04:30-05:00ut历元时刻imart算法和iri2012经验模型重构的电子密度比较情况。把经验模型作为一个重要的参考，则根据图3(a)和图3(b)可知，由imart算法获取的电子密度在整体上更加与iri2012经验模型一致，这也能够说明imart算法在层析过程中反映出的优良性能。

图4(a)-图4(d)为imart算法重构的2016年8月9日的电子密度在经度面111.4°e上的周日变化图，其中，图4(a)对应的时段为04:30-05:00ut，图4(b)对应的时段为08:00-08:30ut，图4(c)对应的时段为16:00-16:30ut，

图4(d)对应的时段为20:30-21:00ut，单位为el/m³，高度和纬度范围分别为70-1000km和24.5-30.5°n。由图4(a)-图4(d)可知，不同区域的电子密度分布特征具有较大差别，从整体上来看，北方的电子密度要小于南方的，因此这也揭示了电子密度分布和纬度之间存在很强的相关性。另外，比较图4(a)-图4(d)可知，电子密度的峰值高度随时间的变化而变化，如从04:30到08:30ut，峰值高度逐渐增加，而从16:30ut开始，峰值高度又开始减小，直到20:30ut时，峰值高度降到205km。因此，上述这些特征反映了电子密度在一天之内垂直高度上的变化情况。

图5(a)和图5(b)分别为2016年8月9日在00:30ut时刻和12:30ut时刻的电子密度轮廓，其电子密度分别来自(111.4°e,27.2°n)电离层测高仪(ionosonde)以及iri2012经验模型、mart算法和imart算法。由图5(a)和图5(b)可知，在不同时刻，由imart算法重构的电子密度轮廓更加接近实测的电子密度轮廓。因此，基于实测数据进一步验证了imart算法的性能。

上述电离层层析方法一方面通过将像素内的截距值与电子含量作为组合量，构造了新的组合量；另一方面改变了松弛因子的获取方式，削弱了噪声对反演结果的影响。并利用全球导航卫星系统实测数据，对该方法在地震监测系统中应用的可靠性和优越性进行了仔细验证，根据统计分析和定量比较的结果，发现该方法表现出显著的重构精度且相较于其他传统层析方法有18％的提高，从而验证了上述电离层层析方法的可靠性及其相对于其他传统层析方法的优越性。

电离层层析装置实施例：

本实施例提供了一种电离层层析装置，包括处理器和存储器，该处理器用于处理存储在存储器中的指令，以实现上述的电离层层析方法。对于本领域的技术人员来讲，可以根据该电离层层析方法，生成相应的计算机指令，以获取该电离层层析装置。由于该电离层层析方法已经在上述的电离层层析方法实施例中进行了详细介绍，此处不再赘述。

最后应当说明的是，以上实施例仅用于说明本发明的技术方案而非对其保护范围的限制，尽管参照上述实施例对本申请进行了详细的说明，所属领域的普通技术人员应当理解，本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或者等同替换，但这些变更、修改或者等同替换，均在本发明的权利要求保护范围之内。