CN102930562A - 电离层的压缩层析成像方法 - Google Patents

Abstract

本发明公开了一种电离层的压缩层析成像方法(Computerized Ionosphere Tomography，CIT)，包括：步骤一：计算投影矩阵；步骤二：基于GPS原始数据计算TEC；步骤三：选择稀疏变换；步骤四：构造稀疏约束的优化目标函数；步骤五：求解步骤四所得的优化函数，完成电离层层析成像。本发明方法解决了传统电离层层析成像由于有限的观测数据导致的成像结果不确定的问题。该技术基于电离层电子密度分布的可压缩性(或可稀疏性)，实现基于有限观测数据的高精度电离层成像。

Description

电离层的压缩层析成像方法

技术领域

本发明涉及雷达成像技术领域，是一种电离层的压缩层析成像方法。

背景技术

地球电离层是一种色散媒质，其特性是时间和空间的函数；会引起穿过其中的低频波段(例如，UHF/VHF、P、L等)无线电信号的相位、幅度及极化产生重要影响；在通信、导航、遥感、隐身/反隐身、全球气候预测与监控、等众多军民用领域有重要的研究和应用价值。目前为止，已发展了众多电离层不规则体的探测与成像的技术和方法，例如，测高仪、顶部探测仪、非相干探测，基于GPS(Global PositionSystem)等观测数据的电离层层析成像，等。如图1和图2所示，根据所处理的电离层不规则体尺度的不同，这些方法可分为确定性和统计两类。统计成像方法主要针对小尺度(米量级)的电离层不规则体；当感兴趣的电离层区域包含大量随时间和空间变化的小尺度不规则体时，这些不规则体的统计参数的重建显得更有意义。确定性方法又可分为射线层析成像方法和衍射层析成像方法；其中以电离层层析成像技术为代表的射线层析成像方法主要针对大尺度的不规则体(数十公里至数百公里)，然而衍射层析成像方法针对中等尺度(数百米至数公里量级)的不规则体。

电离层层析成像基于GPS的信标观测，利用若干空间上分开的TEC(Total Electron Content)监测站，通过断层扫描反演技术获得电离层结构的方法。1986年，美国伊利诺伊大学Austen等首次提出联合应用低轨卫星无线电信标测量和计算机层析成像(ComputerizedTomography，CT)技术反演电离层电子密度二维分布的设想。CIT为传统的电离层无线电探测技术注入了新的活力，使原来仅能提取沿传播路径积分电子密度信息的卫星信标差分多普勒测量数据，可以用于反演电离层电子密度沿水平和垂直方向的二维分布信息。众多科研工作者为CIT的推广做出了努力，其中大部分研究工作是基于GPS信标的CIT技术，一部分研究工作是基于星载SAR(Synthetic ApertureRadar)信号提取电离层TEC信息，然后对电离层结构进行CT成像。CIT已经成为目前最流行的电离层成像方法之一。

通常，电离层层析成像方法基于GPS信标，通过测量电离层观测剖面内多条路径的TEC，获取探测区域内大量彼此交叉路径电子密度的积分数据，然后利用合适的重建算法，由线积分数据反演出剖面内电子密度的分布。以二维电子密度分布为例，CIT探测的几何示意如图3所示。沿某路径的TEC为沿该路径的电子密度积分，根据Austen方法将该电子密度用一组基函数展开。将电离层用立方体进行剖分，选择Pixel基函数：

将GPS与地面之间的电离层探测区域沿垂直高度方向和纬度方向均匀剖分，将连续的电子密度图像离散化，并假定每个网格内部的电子密度相同，离散值密度值为X(i，j)，其中，i＝1，…，n，j＝1…，N。将待反演区域内的n×N个网格按照一定的顺序排列，相应地，离散的电子密度值排列成一个列向量x_k，k＝1，…，n×N，CIT问题转化成求解：

$\stackrel{\‾}{y}=\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}+\stackrel{\‾}{n}---\left(4\right)$

但由于GPS、接收站、电离层区域之间特定的位置关系：(1)GPS轨道很高，它与地面接收机之间的射线主要集中在接近垂直的方向上；(2)地面接收机数量有限且间隔不均。这导致CIT是一个典型的有限视角成像问题，另一方面为降低地面跟踪系统的成本，接收站无法做到细密均匀分布，观测射线分布不均匀，因此有限观测视角CIT的测量数据量远远小于未知数的个数，属于不适定或病态问题。

因此在实际电离层层析成像问题中，方程组(4)通常是病态的，已发展了一些正则化方法来处理该不适定问题，例如迭代方法，典型的有代数重建方法(Algebraic reconstruction technique，ART)，乘法代数重建方法(the Multiplicative algebraic reconstructiontechnique，MART)等；又如最小2-范数或最小能量约束的正则化，尽管最小能量正则化后的不适定问题存在闭式结果，但在许多CIT情况中最小能量的假设并不合理，且该方法依赖于电离层先验模型的选取。

如何利用有限测量数据获得高精度的电离层图像是实现有限视角CIT成像的关键。理论和经验结果表明：自然界的大多数信号和图像在小波等变换域下具有稀疏表示。鉴于此，Donoho、Candès及华裔科学家Tao等人近年来提出了压缩感知理论，该理论表明：对于可稀疏信号(包括信号本身稀疏或者在某组基函数下具有稀疏表示形式)可通过较少的测量数据精确重建原始信号。

发明内容

本发明的目的是公开一种电离层的压缩层析成像方法，将稀疏约束的优化的理论和方法应用到电离层层析成像上，解决了传统基于GPS等卫星信标测得的稀疏电离层TEC观测数据进行电离层层析成像问题中成像分辨率低和成像结果高度倚赖初值选择的瓶颈问题，实现了基于有限观测数据的高精度电离层成像。

为达到上述目的，本发明的技术解决方案是：

一种电离层的压缩层析成像方法，其包括步骤：

步骤一：计算投影矩阵；

步骤二：从GPS原始数据计算TEC；

步骤三：选择稀疏变换；

步骤四：构造稀疏约束的优化函数；

步骤五：求解步骤四所得的优化函数，实现电离层电子密度分布的重构，完成电离层层析成像。

所述的电离层的压缩层析成像方法，其所述步骤三中，包括输出一种稀疏变换

在该变换下，电离层电子密度分布得到很好的压缩；或经过该变换仅存有极少数重要的非零元素。

所述的电离层的压缩层析成像方法，其所述稀疏变换

是经典的离散余弦或小波的正交变换，或是经验上观测到的或通过机器学习训练得到的欠完备或超完备的线性变换。

所述的电离层的压缩层析成像方法，其所述通过机器学习，为通过主分量分析、样本聚类、K-SVD的方法学习。

所述的电离层的压缩层析成像方法，其所述步骤五中求解稀疏约束优化问题，是基于稀疏约束的优化方法对电离层电子密度进行重建，依据电离层电子密度分布的可稀疏性是一个正则化原病态问题的非常重要的先验信息；电离层TEC数据组成的n维观测数据

n×N的投影矩阵

以及二维或三维电离层电子密度组成的N维向量x等三者之间关系由式

表征，其中

观测误差。该病态方程组的求解表征为如下两类稀疏约束的最优化问题：

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1---\left(5\right)$

$s.t.,\stackrel{\‾}{y}=\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}+\stackrel{\‾}{\mathrm{\ϵ}}$

或

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|\stackrel{\‾}{y}-\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}\left|\right|}_2^2---\left(6\right)$

其中，

是稀疏变换设计模块输出的稀疏变换矩阵；λ是平衡数据失配项和稀疏正则化项

的正则化参数，

和

分别为1-范数和2-范数，即，

和

其中a_i是

的第i个元素。

所述的电离层的压缩层析成像方法为求解(5)式或(6)式的方法，包括贪婪算法、迭代硬/软门限、稀疏贝耶斯、线性规划的优化算法。

本发明方法将稀疏信号处理原理和方法应用到有限观测视角的CIT问题中，采用L1范数作为稀疏性度量，将有限观测视角的CIT问题转化为稀疏约束的最优化问题，通过求解最优化问题可以用远远小于未知数个数的测量数据量对电离层图像实现较好地重建。

附图说明

图1为电离层成像方法的分类示意图；

图2为电离层无线电成像方法分类示意图；

图3为电离层层析成像测量结构示意图；

图4为不同方法对含有等离子体泡结构的电离层电子密度分布在110经度处的重建的电离层电子密度随高度变化的重建结果的比较图；其中：

图4(a)为TSVD方法的解；

图4(b)为ART方法的解；

图4(c)为最小二乘最小二范数方法的解；

图4(d)为稀疏约束优化问题(5)的解；

图4(e)为稀疏约束优化问题(5)的解；在这些图中，黑色实线表示重建的结果，圆圈标注的黑色实线表示准确的结果。

图5为不同方法对含有等离子体泡结构的电离层电子密度分布350公里高度处的重建的电离层电子密度随经度变化的结果比较；其中：

图5(a)为TSVD方法的解；

图5(b)为ART方法的解；

图5(c)为最小二乘最小二范数方法的解；

图5(d)为稀疏约束优化问题(5)的解；

图5(e)为稀疏约束优化问题(6)的解；在这些图中，黑色实线表示重建的结果，圆圈标注的黑色实线表示准确的结果。

图6为本发明的一种电离层的压缩层析成像方法流程框图。

具体实施方式

本发明是一种电离层的压缩层析成像方法，其工作原理是电离层电子密度图像在选定的基函数下可呈现较好的稀疏性。以32×32的查普曼背景电离层电子密度为例，仅利用该二维电子密度分布的2％个最大Haar小波系数就可得到精确的重建。电离层电子密度图像的稀疏表示为实现有限视角电离层CIT成像提供了可能和依据。压缩感知理论表明：可通过较少的测量数据实现稀疏图像(或在某个已知基函数集上可稀疏的图像)的精确重建。

对于含有等离子体泡、电离槽、行扰等不规则体结构的电离层的分析表明，电离层电子密度图像的水平方向分布在Harr小波变换域通常呈现良好的稀疏性，垂直分布在离散余弦变换域(DCT)上呈现良好的稀疏性。构造如下两类稀疏约束的优化问题：

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1---\left(7\right)$

$s.t.,\stackrel{\‾}{y}=\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}+\stackrel{\‾}{\mathrm{\ϵ}}$

或

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|\stackrel{\‾}{y}-\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}\left|\right|}_2^2---\left(8\right)$

在公式(7)和(8)中，‖Φx‖₁取为

其中和D₂为分别针对水平方向的Harr小波变换和针对垂直方向的DCT变换，α₁和α₂为正则化参数，

为先验的参考图像。本发明采用共轭梯度方法优化。价格函数中的正则化参数采用广义交叉验证(GeneralizedCross Validation，GCV)方法来优化。

图4给出了以上反演方法对含有等离子体泡结构的电离层电子密度分布在110经度处的重建的电离层电子密度随高度变化的重建结果的比较图，图4(a)TSVD(Truncated SVD)方法的解，图4(b)ART方法的解，图4(c)最小二乘最小二范数方法的解，图4(d)为稀疏约束优化问题(5)的解，图4(e)为稀疏约束优化问题(6)的解；在这些图中，黑色实线表示重建的结果，圆圈标注的黑色实线表示准确的结果。仿真采用峰状电子密度结构，背景电离层采用Chapman模型，信噪比为27.7dB。图5不同方法对含有等离子体泡结构的电离层电子密度分布在350公里高度处的重建的电离层电子密度随经度变化的结果比较，图5(a)TSVD方法的解，图5(b)ART方法的解，图5(c)最小二乘最小二范数方法的解，图5(d)为稀疏约束优化问题(5)的解，图5(e)为稀疏约束优化问题(6)的解；在这些图中，黑色实线表示重建的结果，圆圈标注的黑色实线表示准确的结果。利用本发明的两种优化方法，如图5d、图5e所示，较好的重建出“泡”状结构的形状、位置，电子密度在垂直和水平方向上的分布，与原始值吻合较好。

本发明的一种电离层的压缩层析成像方法，具体执行步骤：

步骤一：确定投影矩阵。

步骤二：测量TEC数据。

步骤三：选择一组合适的基函数将电离层电子密度分布图像稀疏化。

步骤四：构造稀疏约束的优化函数。

步骤五：利用共轭梯度方法和硬门限贪婪算法求解稀疏约束优化问题，反演电离层电子密度分布。

如图6所示，为本发明的一种电离层的压缩层析成像方法流程框图。与传统的电离层层析成像类似，本发明方法流程也包含GPS卫星数据读取、电离层TEC提取、成像区域确定、投影矩阵计算及病态矩阵求解五个基本模块；除病态矩阵求解模块外，其它模块的工作原理和流程与传统的电离层层析成像完全相同。区别与传统电离层层析成像，本发明的方法存在如下两个重要的不同：

1)稀疏变换设计模块

该模块是本发明特有的，它的目的是输出一种稀疏变换Φ。在该变换下，电离层电子密度分布可得到很好的压缩；或者说，经过该变换仅有极少数重要的非零元素。稀疏变换

可以是经典的离散余弦和小波等正交变换，也可以是经验上观测到的或通过机器学习(例如，主分量分析、样本聚类、K-SVD，等方法)训练得到的欠完备或超完备的线性变换。

2)稀疏约束的优化问题求解模块

传统的电离层层析成像中病态矩阵求解模块采用TSVD，ART及其变形，最小二乘-最小二范数求解，等方法；这些方法均高度依赖于初值的选择，否则无法得到满意的重建结果。本发明采用基于稀疏约束的优化方法对电离层电子密度进行重建，其动机是电离层电子密度分布的可稀疏性是一个正则化原病态问题的非常重要的先验信息。电离层TEC数据组成的n维观测数据

n×N的投影矩阵

以及二维或三维电离层电子密度组成的N维向量x等三者之间关系由式

表征。该病态方程组的求解表征为如下两类稀疏约束的最优化问题：

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1---\left(9\right)$

$s.t.,\stackrel{\‾}{y}=\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}+\stackrel{\‾}{\mathrm{\ϵ}}$

或

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|\stackrel{\‾}{y}-\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}\left|\right|}_2^2---\left(10\right)$

其中，

是稀疏变换设计模块输出的稀疏变换矩阵；λ是平衡数据失配项

和稀疏正则化项

的正则化参数，

和

分别为1-范数和2-范数，即，

和

其中a_i是

的第i个元素。

已发展了众多优秀的方法求解优化问题(9)式或(10)式，例如，贪婪算法、迭代硬/软门限、稀疏贝耶斯、线性规划优化方法。

流程图注解：

基于微分多普勒方法电离层TEC估计的工作原理是GPS卫星与地面接受站之间距离的改变或电离层电子浓度的改变将会导致多普勒频移。设频率为f的电磁波在折射率分布为n的介质中传播了距离S后，其相位变化为：

$\mathrm{\θ}=\frac{2\mathrm{\πf}}{c}{\∫}_{s}n\·\mathrm{dS}---\left(11\right)$

其中c为真空中的光速。设接收点距地心的距离为R_E，卫星到地心的距离为r_T。距离r_T的变化和电子密度的变化均会导致(11)式的相位θ产生变化。对式(11)求导可得多普勒频移：

$\mathrm{\Δf}=-\frac{f}{c}\frac{d}{\mathrm{dt}}{\∫}_{R_E}^{r_T}\mathrm{ndS}=-\frac{f}{c}[n{V}_{r_T}+{\∫}_{R_E}^{r_T}\frac{d}{\mathrm{dt}}\mathrm{ndS}]---\left(12\right)$

其中表示卫星在r_T位置时相对于接收机的速度。

GPS卫星发射两个频率f₁和f₂的信号，其中f₂＝mf₁，m为大于1的正整数。由于Doppler效应，地面上接受到的这两个信标的频率差拍是：

$f_b\≈\frac{{\mathrm{mf}}_1}{c}{\∫}_{R_E}^{r_T}\frac{d}{\mathrm{dt}}(n_2-n_1)\·\mathrm{dS}\≈\frac{40.3m}{{\mathrm{cf}}_1}(1-\frac{1}{m^2})\frac{d}{\mathrm{dt}}\mathrm{TEC}---\left(13\right)$

式(13)表明利用测量得到的频率拍频f_b可求得TEC随时间的变化率。

基于法拉第旋转的电离层TEC估计的工作原理总结如下。任一线极化高频电波通过各向异性电离层后极化面总的旋转角为：

$\mathrm{\Ω}\≈\frac{K}{f^2}\stackrel{\‾}{B\cos \mathrm{\θ}\sec \mathrm{\χ}}{\∫}_{R_E}^{r_T}\mathrm{Ndr}=\frac{K}{f^2}\stackrel{\‾}{B\cos \mathrm{\θ}\sec \mathrm{\χ}}\·\mathrm{TEC}---\left(14\right)$

其中，θ是射线方向与地磁场的夹角，χ是射线的天顶角，dS＝secχ·dr，B是射线路径上质心高度处的地磁场高度。式(14)中的K是常数，

e为电子电量，m为电子质量，c为光速，ε₀为真空中的介电常数。

M是Bcosθsecχ在路径上的加权平均值，它依赖于射线路径相对磁场的方向以及电子浓度剖面，对于给定的磁场模式是可计算的。式(14)表明只要测得法拉第旋转角Ω，就可以求得TEC。

基于群/相时延的电离层TEC工作原理如下。电磁波在电离层中传播的群速度v_g和相速度v_p分别为：

$v_g=\frac{\mathrm{d\ω}}{\mathrm{dk}}=\mathrm{cn}=c\sqrt{1-{({\mathrm{\ω}}_p/\mathrm{\ω})}^2}\≈c(1-40.3\frac{N_e}{f^2})---\left(15\right)$

$v_p=\frac{\mathrm{\ω}}{k}=\frac{c}{n}=\frac{c}{\sqrt{1-{({\mathrm{\ω}}_p/\mathrm{\ω})}^2}}\≈c(1+40.3\frac{N_e}{f^2})---\left(16\right)$

对应GPS卫星两个不同频率f₁和f₂的伪距(L＝∫_sv_gdt)的差值ΔS与电离层TEC之间的关系是：

$\mathrm{TEC}=\frac{\mathrm{\ΔS}}{40.3(\frac{1}{{f_1}^2}-\frac{1}{{f_2}^2})}---\left(17\right)$

对应GPS卫星两个不同频率f₁和f₂的相位差值与电离层TEC之间的关系是：

$\mathrm{TEC}=\frac{c}{40.3}\frac{{f_1}^2{f_2}^2}{{f_1}^2-{f_2}^2}(\frac{{\mathrm{\Φ}}_1+2\mathrm{\π}{m}_1}{f_1}-\frac{{\mathrm{\Φ}}_2+2\mathrm{\π}{m}_2}{f_2})---\left(18\right)$

其中m₁和m₂是相位测量的模糊度。

Claims (9)

1.一种电离层的压缩层析成像方法，其特征在于通过稀疏的电离层TEC观测数据实现电离层电子密度的高分辨率成像，包括步骤：

步骤一：计算投影矩阵；

步骤二：基于GPS原始数据计算TEC；

步骤三：选择稀疏变换；

步骤四：构造稀疏约束的优化目标函数；

步骤五：求解步骤四所得的优化函数，完成电离层层析成像。

2.如权利要求1所述的电离层的压缩层析成像方法，其特征在于，所述步骤一的投影矩阵的计算流程如下：首先确定反演区域，然后将待反演区域内连续分布的电子密度离散化，并假设在每一个像素网格内的电子密度是相同的；每一次TEC测量事件都对应投影矩阵中的一行，投影矩阵的行数对应GPS卫星位置的采样数与地面接收台站个数的乘积；投影矩阵每一行元素值对应与该行相应的射线穿过待反演区域被网格区域所截的线段长度；因此，电离层CT成像转化成求解线性方程组

其中

表示观测误差，

表示投影矩阵，该矩阵元素的值等于各个网格所截射线线段的长度值；

代表原始图像的一维列向量，表示待重建区域电子密度分布的离散网格值；

表示沿GPS-地面接收站之间射线路径进行线积分的TEC数据，向量

的第i个元素对应第i条射线路径的TEC。

3.如权利要求1所述的电离层的压缩层析成像方法，其特征在于，所述TEC计算方法包括微分多普勒方法、群/相时延方法或法拉第旋转方法。

4.如权利要求1所述的电离层的压缩层析成像方法，其特征在于，所述步骤三中，包括输出一种稀疏变换

在该变换下，电离层电子密度分布可得到很好的压缩；或经过该变换后仅有极少数重要的非零元素。

5.如权利要求4所述的电离层的压缩层析成像方法，其特征在于，所述稀疏变换

是经典的离散余弦或小波的正交变换，或是经验上观测到的或通过机器学习训练得到的欠完备或超完备的线性变换。

6.如权利要求5所述的电离层的压缩层析成像方法，其特征在于，所述通过机器学习，为通过主分量分析、样本聚类或K-SVD的方法学习。

7.如权利要求1所述的电离层的压缩层析成像方法，其特征在于，所述步骤五中求解稀疏约束优化问题，是基于稀疏约束的优化方法对电离层电子密度进行重建，其依据是电离层电子密度分布的可稀疏性是一个正则化原病态问题的非常重要的先验信息。

8.如权利要求1所述的电离层的压缩层析成像方法，其特征在于，所述电离层TEC数据组成的n维观测数据

n×N的投影矩阵

以及二维或三维电离层电子密度组成的N维向量x等三者之间关系由

表征，其中

观测误差。该病态方程组的求解表征为如下两类稀疏约束的最优化问题：

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1---\left(1\right)$

$s.t.,\stackrel{\‾}{y}=\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}+\stackrel{\‾}{\mathrm{\ϵ}}$

或

${\min }_x{\left|\right|\stackrel{\‾}{\stackrel{\‾}{\mathrm{\Φ}}}\stackrel{\‾}{x}\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|\stackrel{\‾}{y}-\stackrel{\‾}{\stackrel{\‾}{A}}\stackrel{\‾}{x}\left|\right|}_2^2---\left(2\right)$

其中，

是稀疏变换设计模块输出的稀疏变换矩阵；λ是平衡数据失配项

和稀疏正则化项

的正则化参数，

和

分别为1-范数和2-范数，即，

和

其中a_i是

的第i个元素。

9.如权利要求7所述的电离层的压缩层析成像方法，其特征在于，所述求解式(1)或式(2)的优化方法包括贪婪算法、迭代硬/软门限、稀疏贝耶斯或线性规划算法。

Images