CN111580061B

CN111580061B - 基于clean算法的电离层电子密度反演方法

Info

Publication number: CN111580061B
Application number: CN202010430957.5A
Authority: CN
Inventors: 郝红连; 赵必强; 乐新安; 丁锋; 曾令旗
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2020-05-20
Filing date: 2020-05-20
Publication date: 2020-10-27
Anticipated expiration: 2040-05-20
Also published as: CN111580061A

Abstract

本发明属于信号与信息处理领域，具体涉及一种基于CLEAN算法的电离层电子密度反演方法、系统、装置，旨在解决非相干散射等离子体谱线计算量大以及电子密度剖面反演精度不够的问题。本系统方法包括：获取IQ数字信号；提取并删除IQ数字信号中脉冲发射期间采集的信号部分，构建二维矩阵；通过频域FFT算法解码计算各距离门高度的信号功率谱；迭代获取各距离门高度的功率谱数据进行累加，并通过系统频率响应函数进行去噪；基于CLEAN算法对去噪后的功率谱数据解卷积；通过样条插值法获取整个高度剖面的等离子体谱线，并利用朗缪尔色散关系拟合得到电离层电子密度剖面。本发明提高了等离子体谱线计算的实时性和电子密度剖面反演的精度。

Description

基于CLEAN算法的电离层电子密度反演方法

技术领域

本发明属于信号与信息处理领域，具体涉及一种基于CLEAN算法的电离层电子密度反演方法、系统、装置。

背景技术

在电离层研究中使用超大天线增益的大功率甚高频/超高频雷达系统叫做非相干散射雷达。非相干散射的机理是发射的雷达脉冲入射到电离层中每个自由电子上时，会引起电子振荡并加速电子，从而引起电磁能的再辐射，这种过程称作汤姆逊散射或非相干散射。因此非相干散射信号谱是非相干散射雷达对电离层电子密度的热起伏引起电磁散射的统计测量。

非相干散射雷达发射电磁波信号后散射回来的回波信号功率谱由具有双峰结构的窄带离子谱线和较弱更窄的上移、下移等离子体谱线组成，它们分别是与离子声波和朗缪尔波的色散关系有关。离子谱线中包含有等离子体漂移速度、电子密度、电子/离子温度和离子成分等信息。通过离子声波散射的离子谱线来获取的测量精度受到雷达系统的状态、积累时间、脉冲编码方式等多种因素的影响，使得某些电离层参量如电子密度时空分辨率变得比较差，还需要进一步的标校和定标。那么如果能准确获取朗缪尔波散射的等离子体谱线，然后提取等离子体谱线每个高度的频偏，利用朗缪尔色散关系，就可以反演得到精度更高的电子密度剖面，这样得到的电子密度信息也可以作为离子谱线散射模型的先验信息。但是由于等离子体谱线通常位于信号功率谱中心频带外的1～15MHz处，带宽非常窄并且很微弱，而且散射信号在被接收采样后带来的距离模糊使得更加难以从回波信号中分离出等离子体谱线。因此等离子体谱线的测量需要非相干散射雷达接收系统带宽要大于等于两倍等离子体频率以及较好的频谱分辨率。同时要在整个高度剖面范围内进行数据采样并获得所需的频谱分辨率，尤其是对于实时计算来说，就需要在数据处理方面也要具备较强的数据计算能力。

CLEAN算法是一种在去除特定频率信号的同时将信号副瓣也可以同时去除的解卷积技术，在改善图像质量方面有广泛应用。在常规硬目标雷达回波处理中常用来抑制杂噪信号。在非相干散射雷达信号处理领域，将CLEAN算法应用于等离子体谱线的提取问题目前还没有研究。实际上，计算实测的非相干散射回波信号的自相关过程就是某个高度的等离子体自相关与对应模糊函数的卷积过程，因此对实测回波信号自相关解模糊的过程其实就是对实测回波信号的解卷积。在频域，就是实测回波信号的功率谱等于等离子体谱线与模糊函数做傅里叶变换后的乘积。因此，本发明提出了一种基于CLEAN算法的电离层电子密度反演方法、系统、装置。

发明内容

为了解决现有技术中的上述问题，即为了解决现有的非相干散射等离子体谱线计算量大以及电子密度剖面反演精度不够的问题，本发明第一方面，提出了一种基于CLEAN算法的电离层电子密度反演方法，该方法包括：

步骤S100，获取一个脉冲重复周期内电离层散射的回波信号并进行下变频，得到IQ数字信号；所述回波信号包括脉冲发射期间，雷达接收机采样的回波信号；

步骤S200，提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各门的FFT点数，通过预设的第一方法构建二维矩阵；

步骤S300，基于所述二维矩阵中各距离门对应的信号数据，通过频域FFT算法解码计算各距离门高度的信号功率谱，作为第一功率谱数据；

步骤S400，循环执行步骤S100-步骤S300，获取W个脉冲重复周期内各距离门高度的第一功率谱数据并进行累加，作为第二功率谱数据；通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪，得到第三功率谱数据；其中，W为正整数；

步骤S500，基于CLEAN算法，结合脉冲发射期间采样的回波信号，对第三功率谱数据进行解卷积，得到第四功率谱数据；

步骤S600，获取第四功率谱数据中每个高度的等离子体频率，并通过样条插值法获取整个高度剖面的等离子体谱线；基于所述等离子体谱线，利用朗缪尔色散关系拟合得到电离层的电子密度剖面。

在一些优选的实施方式中，步骤S200中“提取并删除IQ数字信号中脉冲发射期间采集的信号部分”，其方法为：

计算IQ数字信号各采样点的瞬时功率值，作为第一功率值；

结合多个采样点的第一功率值求平均，得到IQ数字信号各采样点对应的第二功率值；

计算各采样点第二功率值之间的差，将差大于设定的门限值的数据点的位置作为脉冲发射期间，雷达采样信号的起始点位置，将差小于设定的门限值的数据点的位置作为脉冲发射期间，雷达采样信号的截止点位置；

计算所述起始点位置与所述截止点位置的差的绝对值，作为脉冲发射期间，雷达采样的总的信号数据点个数。

在一些优选的实施方式中，步骤S200中“通过预设的第一方法构建二维矩阵”，其方法为：

步骤S210，设脉冲发射期间，雷达采样的总的信号数据点个数为M；

步骤S211，将M+1个数据点至M+M个数据点的信号作为二维矩阵第一行元素的前M个元素，剩余元素置零，组成第一个距离门要进行FFT计算的数据；

步骤S212，令M＝M+1，根据预设的距离门个数，迭代执行步骤S211，构建二维矩阵。

在一些优选的实施方式中，步骤S300中“通过频域FFT算法解码计算各距离门高度的信号功率谱”，其方法为：

其中，

表示各距离门对应的信号功率谱，X_N(e^jω)为各距离门对应的频谱，N表示各距离门对应的信号数据点个数，x_N(n)表示各距离门对应的有限长数据序列，n表示各数据点下标。

在一些优选的实施方式中，若步骤S400中“通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪”其方法为：

其中，|H(ω)|²表示雷达接收机系统的频率响应，

表示去噪后的功率谱数据。

在一些优选的实施方式中，步骤S500中“对第三功率谱数据进行解卷积”，其方法为；

步骤S510，将发射期间采样的回波信号与移位后的回波信号共轭相乘，构建时域数据矩阵；

步骤S520，对所述时域数据矩阵中的数据进行傅里叶变换，并进行归一化，将归一化后的矩阵作为第一矩阵；

步骤S530，构建与第三功率谱数据矩阵大小相同的矩阵，作为第二矩阵，并将该矩阵中的各元素初始化为0；构建与第一矩阵大小相同的矩阵，作为第三矩阵，其中该矩阵元素构成的频谱中心主瓣通过高斯拟合得到，其余元素初始化为0；

步骤S540，对所述第三功率谱数据进行归一化后，获取该数据矩阵中幅值最大的数据点的幅值以及该数据点的位置；将第一矩阵的中心位置对准归一化后第三功率谱数据中该最大值点位置，找出二者数据矩阵重叠的部分，并与设定的循环比例因子相乘后从第三功率谱数据矩阵中减去，同时将第三矩阵中对应位置范围的数据与所述最大幅值、设定的循环比例因子相乘作为第二矩阵中对应位置处的数据；

步骤S550，循环执行步骤S540，直至第三功率谱数据剩余的部分接近背景噪声或迭代次数达到最大阈值，则终止循环，并将第二矩阵与第三功率谱数据剩余部分相加得到第四功率谱数据。

在一些优选的实施方式中，所述朗缪尔色散关系其计算过程为:

其中，ω_L表示朗缪尔波频率，ω_pe表示电子等离子体频率，ω_ce表示电子回旋频率，v_Te表示电子热速度，θ表示雷达视线和地磁场之间的夹角，k表示朗缪尔波的波数。

本发明的第二方面，提出了一种基于CLEAN算法的电离层电子密度反演系统，该系统包括信号获取模块、构建矩阵模块、频域解码模块、去噪模块、解卷积模块、拟合模块；

所述信号获取模块，配置为获取一个脉冲重复周期内电离层散射的回波信号并进行下变频，得到IQ数字信号；所述回波信号包括脉冲发射期间，雷达接收机采样的回波信号；

所述构建矩阵模块，配置为提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各距离门的FFT点数，通过预设的第一方法构建二维矩阵；

所述频域解码模块，配置为基于所述二维矩阵中各距离门对应的信号数据，通过频域FFT算法解码计算各距离门高度的信号功率谱，作为第一功率谱数据；

所述去噪模块，配置为循环执行信号获取模块-频域解码模块，获取W个脉冲重复周期内各距离门高度的第一功率谱数据并进行累加，作为第二功率谱数据；通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪，得到第三功率谱数据；其中，W为正整数；

所述解卷积模块，配置为基于CLEAN算法，结合脉冲发射期间采样的回波信号，对第三功率谱数据进行解卷积，得到第四功率谱数据；

所述拟合模块，配置为获取第四功率谱数据中每个高度的等离子体频率，并通过样条插值法获取整个高度剖面的等离子体谱线；基于所述等离子体谱线，利用朗缪尔色散关系拟合得到电离层的电子密度剖面。

本发明的第三方面，提出了一种存储装置，其中存储有多条程序，所述程序应用由处理器加载并执行以实现上述的基于CLEAN算法的电离层电子密度反演方法。

本发明的第四方面，提出了一种处理装置，包括处理器、存储装置；处理器，适用于执行各条程序；存储装置，适用于存储多条程序；所述程序适用于由处理器加载并执行以实现上述的基于CLEAN算法的电离层电子密度反演方法。

本发明的有益效果：

本发明提高了等离子体谱线计算的实时性和电子密度剖面反演的精度。本发明采用频域FFT解码计算等离子体谱线方法，大大降低了等离子体谱线计算的时间和空间复杂度，并基于GPU-CUDA平台的并行处理，加快了实测的非相干散射等离子体谱线的计算速度。

同时，利用CLEAN算法以及采样发射期间的回波信号对信噪比很差的功率谱进行了解卷积处理，从卷积失真中恢复得到的真实的等离子体谱线图中准确提取每个高度的等离子体频率，从而提高了电子密度反演的精度。

附图说明

通过阅读参照以下附图所做的对非限制性实施例所做的详细描述，本申请的其他特征、目的和优点将会变得更明显。

图1是本发明一种实施例的基于CLEAN算法的电离层电子密度反演方法的流程示意图；

图2是本发明一种实施例的基于CLEAN算法的电离层电子密度反演系统的框架示意图；

图3是本发明一种实施例的基于CLEAN算法的电离层电子密度反演方法的详细流程示意图；

图4是本发明一种实施例的非相干散射回波信号的时域数据矩阵重组的示意图；

图5是本发明一种实施例的雷达发射和接收的示意图；

图6是本发明一种实施例的基于CLEAN算法进行解模糊的流程示意图；

图7是本发明一种实施例的去背景噪声前的功率谱效果示意图；

图8是本发明一种实施例的去背景噪声后的功率谱效果示意图；

图9是本发明一种实施例的基于CLEAN算法对上移等离子体谱线解模糊前后的对比示意图；

图10是本发明一种实施例的基于CLEAN算法对下移等离子体谱线解模糊前后的对比示意图；

图11是本发明一种实施例的电子密度剖面反演结果的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

本发明的基于CLEAN算法的电离层电子密度反演方法，如图1和图3所示，包括以下步骤：

步骤S200，提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各距离门的FFT点数，通过预设的第一方法构建二维矩阵；

为了更清晰地对本发明基于CLEAN算法的电离层电子密度反演方法进行说明，下面结合附图对本发明方法一种实施例中各步骤进行展开详述。

步骤S100，获取一个脉冲重复周期内电离层散射的回波信号并进行下变频，得到IQ数字信号；所述回波信号包括脉冲发射期间，雷达接收机采样的回波信号。

在本实施例中，天线阵面接收电离层散射的回波信号，对回波信号进行下变频得到中频模拟信号，对中频模拟信号进行AD采样，变为中频数字信号。对中频数字信号进行数字正交下变频，得到IQ数字信号。

步骤S200，提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各距离门的FFT点数，通过预设的第一方法构建二维矩阵。

雷达接收机系统为了能够完整地采样发射期间的回波信号，接收机开始采样的时间通常要比脉冲开始发射时间提前，这样就会导致接收机刚开始采样的信号都是系统的噪声信号。因此需要从每个IPP(脉冲重复周期)内的回波数据流中检测出真正开始采样发射期间回波信号的数据点位置，可以结合平方律包络检波原理，通过记录包络检波检测到的发射期间回波信号包络的有效上升沿和下降沿到达的时间，从而找到采样发射期间回波信号的起始点和截止点。记录发射期间采样的回波信号长度可以从回波信号数据流中提取出发射期间的回波信号，用来后续的解卷积处理，同时也可以验证发射信号编码调制的正确性。检测并提取脉冲发射期间的信号的具体方法如下：

步骤S201，计算IQ数字信号各采样点的瞬时幅值的平方，即瞬时功率值，作为第一功率值，如公式(1)所示：

X(n)＝I(n)²+Q(n)²,n＝1,2,.....N (1)

其中，I(n)²表示IQ数字信号中I路的瞬时功率值，Q(n)²表示IQ数字信号中Q路的瞬时功率值，n表示采样点下标，N表示采样点的数量，X(n)表示IQ数字信号的瞬时功率值。

步骤S202，由于噪声的影响，若只利用单个采样点功率值来检测难免会有偏差。我们可以利用下式累加多个点的数据值并求平均，这样得到误差较小的瞬时功率值，从而降低对信噪比的要求，即对多个IQ数字信号采样点的第一功率值求平均，得到IQ数字信号各采样点对应的第二功率值，如公式(2)所示：

其中，m表示预设的参数，为正整数，X′(n′)为求平均后的瞬时功率值，即第二功率值。

步骤S203，计算各数据点第二功率值之间的差，寻找功率差大于设定门限值的数据点位置，该点位置即为采样的发射期间回波信号的起始点位置，同样寻找功率差小于设定门限值的数据点位置，该点位置即为采样的发射期间回波信号的截止点位置，也是雷达开始采样距离盲区外的回波信号起始点位置。具体如公式(3)所示：

其中，Threshold1为预设的第一门限值，Threshold2为预设的第二门限值。

步骤S204，计算所述起始点位置与所述截止点位置的差的绝对值，作为脉冲发射期间，雷达采样的总的信号数据点个数。

由于采样的IQ回波信号中包含有发射期间采样的回波信号，属于雷达盲距内的回波信号采样，因此通常在实际距离门分析时，要先剔除掉雷达盲距内采样的回波数据点，再进行每个距离门的数据重组。

具体方法如下：

步骤S212，令M＝M+1，根据预设的距离门个数，迭代执行,步骤S211，构建二维矩阵(矩阵大小为距离门个数×每个距离门的FFT点数)。

以本发明实施例为例，假设发射脉冲宽度为500us，采样间隔为40ns，IPP(T)为20ms，那么在一个IPP内采样的总的数据点个数为500000个复IQ原始数据，其中前12500个为发射期间采样的数据点，因此要剔除掉前12500个数据点，从剩下的数据点中根据预设的距离门个数以及FFT点数重新排列数据为二维矩阵。数据重组过程如图4所示，先将从第V12501到第V25000个复数据放到二维矩阵的第一行，接下来跳过一个数据采样点，将后12500个复数据即V12502到V25001个复数据放到二维矩阵的第二行，依此类推完成数据重排列过程。

步骤S300，基于所述二维矩阵中各距离门对应的信号数据，通过频域FFT算法解码计算各距离门高度的信号功率谱，作为第一功率谱数据。

通常我们假设接收的非相干散射回波信号都是宽平稳随机信号，那么就可以运用周期图法估计真实的功率谱，在本实施例中，由获得的N点数据构成的有限长序列x_N(n)直接求傅里叶变换，得频谱；取频谱幅度的平方，并处以N，以此作为对x(n)真实功率谱的估计。具体如公式(4)(5)所示：

其中，

表示各距离门对应的信号功率谱，X_N(e^jω)为各距离门对应的频谱，N表示各距离门对应的信号数据点个数，x_N(n)表示各距离门对应的有限长数据序列，n表示各数据点下标。上述公式中用FFT来完成傅里叶变换。

CUDA是一种并行编程模型和软件环境，采用Python语言等标准语言进行操作。该技术封装了GPU的硬件细节，支持数据的并行化和线程并行化处理，在可接受的时间内，使得由非相干原始数据到高分辨率宽频带功率谱的计算以及谱线的准确提取成为可能。在本发明实施例中优选使用PyCUDA实现GPU-CUDA并行化计算，PyCUDA为我们提供了一个简单的使用Nvidia CUDA并行计算的接口API，在同一个代码块中允许有许多线程，与同一块中的其他线程相比每个线程都有它ID，它们具有唯一标识符，每个代码块都有它拥有唯一的ID，以便每个代码块中的线程不会被意外地视为线程另一个代码块。对于每个块的大小，它可以显示每个块中有多少线程，线程被放置到块中以存储线程跟踪信息，而不需要使用大量内存。因此基于CUDA的并行多任务处理方案，可以高效实现大规模的FFT计算，并可以同时完成多个距离门的功率谱计算，提高了整个反演方法的实时性。

步骤S400，循环执行步骤S100-步骤S300，获取W个脉冲重复周期内各距离门高度的第一功率谱数据并进行累加，作为第二功率谱数据；通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪，得到第三功率谱数据；其中，W为正整数。

在本实施例中，根据预设的积累周期个数，循环执行步骤S100-步骤S300获取设定周期内对应高度的功率谱数据并进行多周期积累，得到对应距离门的功率谱数据，从而提高信号的信噪比。

另外在接收机接收采样的过程中，回波信号会受到宇宙噪声的影响以及进入后检测滤波后通常会受到接收机系统本身的电磁振荡导致的热噪声的影响，一般将宇宙噪声和接收机系统热噪声都统称为系统噪声，这些噪声都可以看作是白噪声，因此要对实测功率谱进行系统噪声去除。根据白噪声功率谱特性，其功率谱在整个频域内是均匀分布的，但是随着散射信号经过接收机后检测滤波系统后，噪声功率谱形状近似等于接收机系统的频率响应。通过对最高高度处一些距离门对应的功率谱求平均得到噪声功率谱，也就是接收机系统幅频响应|H(ω)|²，那么去噪后的实测功率谱获取过程如公式(6)所示：

其中，

表示去噪后的功率谱数据。

去噪前后的功率谱如图7、8所示，可以看出采用本发明的去背景噪声方法效果明显。

步骤S500，基于CLEAN算法，结合脉冲发射期间采样的回波信号，对第三功率谱数据进行解卷积，得到第四功率谱数据。

根据维纳辛钦定理，要计算非相干散射信号的功率谱，可以先计算非相干散射回波信号的自相关。那么计算非相干散射回波信号自相关的过程化简后可以看作是等离子体自相关与不同时延处对应模糊函数的卷积过程，如式(7)所示，因此可以通过对时延剖面解卷积消除掉每个距离门处的距离模糊。

其中，

为模糊函数，

为距离门步进值，c为光速，env为发射脉冲包络，h为接收机滤波器的脉冲响应，*表示卷积操作，

n(iΔr,t)为第i个距离门、t时刻的电子密度扰动，k为散射波矢，z表示接收的散射信号，j表示距离门索引，x表示电子扰动的位移大小，τ表示时延。

基于以上所述，那么模糊函数可以看作是一个滤波系统，等离子体自相关函数为系统输入，实测回波信号的自相关函数为系统输出，将发射的宽脉冲产生的这种距离模糊影响可以看作是由于非相干散射信号经过此滤波系统后多个高度加权平均的结果。如图5所示，假设雷达发射信号为4位的编码信号，分别为a₀,a₁,a₂,a₃，接收机采样的散射信号为V₀,V₁,V₂,V₃，那么散射回波信号V_n不仅包括高度h₀的散射信号，还包括来自于高度h₁到h₃和高度h_-1到h_-3的散射信号，高度h₀位于倾斜的4×4的正方形中心。为了解决只提取出高度h₀处的散射信号，本发明采用CLEAN算法对实测信号的时延剖面进行解卷积处理从而消除掉距离模糊，同时使用解卷积处理方法，不需要对反演参数强加不切实际的假设，例如在积累时间内假设电离层参量在整个距离门内认为是常量等。

在本实施例中，基于CLEAN算法进行解卷积(解模糊)的流程如图6所示，具体步骤如下：

步骤S510，将发射期间采样的回波信号与移位(位置平移)后的回波信号共轭相乘，构建时域数据矩阵；

构造卷积过程中的模糊函数时域数据矩阵，即CLEAN算法中的时域“脏”光束。如图5中的倾斜的正方形所示，通过每次从中心到矩阵的顶部和底部，将采样的发射期间回波信号与移位后的回波信号共轭相乘，以至于使得乘积的总长度依次递减1，形成2M×M的时域数据矩阵P，例如一共10个数，移位1位，那么与交叠部分共轭相乘，即第1到9个数与第2到8个数共轭相乘，移位2位，就是第1到8个数与第3到第7个数共轭相乘。

对矩阵P中的每一行时域数据做傅里叶变换，得到的频域数据矩阵即为CLEAN算法中的频域“脏”光束，并进行归一化处理得到归一化的“脏”光束B。

步骤S530，构建与第三功率谱数据矩阵大小相等的矩阵，作为第二矩阵，并将该矩阵中的各元素初始化为0；构建与第一矩阵大小相同的矩阵，作为第三矩阵，其中该矩阵元素构成的频谱中心主瓣通过高斯拟合得到，其余元素初始化为0。

所述第二矩阵,即CLEAN算法中的“清洁”图S，所述第三矩阵是与“脏”光束B大小相同的矩阵，即理想“清洁”光束B_c。

步骤S540，对所述第三功率谱数据进行归一化后，获取该数据矩阵中幅值最大的数据点的幅值以及该数据点的位置，然后将第一矩阵的中心位置对准归一化后第三功率谱数据中该最大值点位置，找出二者数据矩阵重叠的部分，并与设定的循环比例因子相乘后从第三功率谱数据矩阵中减去，同时将第三矩阵中对应位置范围的数据与所述最大幅值、设定的循环比例因子相乘作为第二矩阵中对应位置处的数据。

上述计算得到的第三功率谱数据，即CLEAN算法中的“脏”图，进行归一化处理得到归一化的初始“脏”图即初始“脏”图S₀，从初始“脏”图S₀(迭代后为S_q，q为迭代次数)中筛选出最大幅值的数据点，并记录最大幅值

和位置(i,j)。将“脏”束的中心对准“脏”图S₀中的最大幅值位置，找出“脏”图中要去除的频率范围B_i,j乘以循环步长比例因子γ，然后从“脏”图中减去，即S_q-γB_i,j。同时将“清洁”光束B_c中对应的频率范围乘以最大幅值

与循环步长比例因子γ后放进“清洁”图S的对应位置。

图9为CLEAN算法对上移等离子体谱线的解模糊前后的对比图；图10为CLEAN算法对下移等离子体谱线的解模糊前后的对比图，由此可以看出采用CLEAN算法可以有效地去除距离模糊对等离子体谱线带来的影响。

步骤S550，循环执行步骤S540，直至第三功率谱数据剩余的部分接近背景噪声或迭代次数达到最大阈值，则终止循环，并将第二数据矩阵与第三功率谱数据剩余部分相加得到第四功率谱数据。

所述步骤S550中迭代循环后第三功率谱数据剩余的部分，即CLEAN算法中的“残”图S_res将“清洁”图S与“残”图S_res相加重构得到经过CLEAN算法处理干净的等离子体谱线图。

将真实的等离子体谱图从卷积失真中恢复后，可以用Argmax函数逼近功率谱，即提取每个高度上等离子体谱线中功率最大值所对应的频率，这个频率就叫等离子体频率。

在本实施例中，获取第四功率谱数据中每个高度的等离子体频率，并通过样条插值法获取整个高度剖面的等离子体谱线。基于所述等离子体谱线，利用朗缪尔色散关系拟合得到电离层的电子密度剖面。

在磁化等离子体中，朗缪尔波的色散关系如公式(8)所示：

其中，ω_L表示朗缪尔波频率，ω_pe表示电子等离子体频率，

ω_ce表示电子回旋频率，ω_ce＝eB₀/m_e，v_Te表示电子热速度，

θ表示雷达视线和地磁场之间的夹角，k表示朗缪尔波的波数，B₀表示地磁场强度，m_e,n_e,e,T_e分别表示电子质量，电子密度，电子电量以及电子温度，ε₀表示自由空间的电导率，k_b表示玻尔兹曼常量。式(8)进一步推导得出电子密度，如公式(9)所示：

然而，在给定的积累时间内，等离子体谱线在整个高度剖面上应该是连续的。本发明实施例表明，当每个高度独立拟合功率谱时，沿曲线的连续点之间存在振荡。因此，为了保证等离子体谱线的连续性和平滑性，可以进一步采用样条插值，使拟合得到的等离子体谱线更加平滑。通过从等离子体谱线的每个高度获得等离子体频率后再根据式(8)反演得到电子密度剖面。

同时，本发明的效果通过Arecibo雷达的长脉冲实验进行了仿真验证。其中参数为脉冲宽度为500us的长脉冲，IPP为20ms，积累时间10s。

图11为拟合得到的电子密度剖面，从图中可以看出电子密度清晰平滑，并且具有较高的时空分辨率。

本发明第二实施例的一种基于CLEAN算法的电离层电子密度反演系统，如图2所示，包括：信号获取模块100、构建矩阵模块200、频域解码模块300、去噪模块400、解卷积模块500、拟合模块600；

所述信号获取模块100，配置为获取一个脉冲重复周期内电离层散射的回波信号并进行下变频，得到IQ数字信号；所述回波信号包括脉冲发射期间，雷达接收机采样的回波信号；

所述构建矩阵模块200，配置为提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各门的FFT点数，通过预设的第一方法构建二维矩阵；

所述频域解码模块300，配置为基于所述二维矩阵中各距离门对应的信号数据，通过频域FFT算法解码计算各距离门高度的信号功率谱，作为第一功率谱数据；

所述去噪模块400，配置为循环执行信号获取模块-频域解码模块，获取W个脉冲重复周期内各距离门高度的第一功率谱数据并进行累加，作为第二功率谱数据；通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪，得到第三功率谱数据；其中，W为正整数；

所述解卷积模块500，配置为基于CLEAN算法，结合脉冲发射期间采样的回波信号，对第三功率谱数据进行解卷积，得到第四功率谱数据；

所述拟合模块600，配置为获取第四功率谱数据中每个高度的等离子体频率，并通过样条插值法获取整个高度剖面的等离子体谱线；基于所述等离子体谱线，利用朗缪尔色散关系拟合得到电离层的电子密度剖面。

所述技术领域的技术人员可以清楚的了解到，为描述的方便和简洁，上述描述的系统的具体的工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。

需要说明的是，上述实施例提供的基于CLEAN算法的电离层电子密度反演系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。

本发明第三实施例的一种存储装置，其中存储有多条程序，所述程序适用于由处理器加载并实现上述的基于CLEAN算法的电离层电子密度反演方法。

本发明第四实施例的一种处理装置，包括处理器、存储装置；处理器，适于执行各条程序；存储装置，适于存储多条程序；所述程序适于由处理器加载并执行以实现上述的基于CLEAN算法的电离层电子密度反演方法。

所述技术领域的技术人员可以清楚的了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实例中的对应过程，在此不再赘述。

本领域技术人员应该能够意识到，结合本文中所公开的实施例描述的各示例的模块、方法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，软件模块、方法步骤对应的程序可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。为了清楚地说明电子硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以电子硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。本领域技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。

术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。

至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims

1.一种基于CLEAN算法的电离层电子密度反演方法，其特征在于，该方法包括：

步骤S500，基于CLEAN算法，结合脉冲发射期间采样的回波信号，对各第三功率谱数据进行解卷积，得到第四功率谱数据；

2.根据权利要求1所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，步骤S200中“提取并删除IQ数字信号中脉冲发射期间采集的信号部分”，其方法为：

计算IQ数字信号各采样点的瞬时功率值，作为第一功率值；

3.根据权利要求1所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，步骤S200中“通过预设的第一方法构建二维矩阵”，其方法为：

4.根据权利要求1所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，步骤S300中“通过频域FFT算法解码计算各距离门高度的信号功率谱”，其方法为：

其中，

5.根据权利要求4所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，若步骤S400中“通过雷达接收机系统的频率响应函数对第二功率谱数据进行去噪”其方法为：

其中，|H(ω)|²表示雷达接收机系统的频率响应，

表示去噪后的功率谱数据。

6.根据权利要求1所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，步骤S500中“对第三功率谱数据进行解卷积”，其方法为：

步骤S530，构建与第三功率谱数据矩阵大小相等的矩阵，作为第二矩阵，并将该矩阵中的各元素初始化为0；构建与第一矩阵大小相同的矩阵，作为第三矩阵，其中该矩阵元素构成的频谱中心主瓣通过高斯拟合得到，其余元素初始化为0；

7.根据权利要求1所述的基于CLEAN算法的电离层电子密度反演方法，其特征在于，所述朗缪尔色散关系其计算方法为：

8.一种基于CLEAN算法的电离层电子密度反演系统，其特征在于，该系统包括：信号获取模块、构建矩阵模块、频域解码模块、去噪模块、解卷积模块、拟合模块；

所述构建矩阵模块，配置为提取并删除IQ数字信号中脉冲发射期间采集的信号部分；基于剩余的信号数据，结合预设的距离门个数及各门的FFT点数，通过预设的第一方法构建二维矩阵；

所述拟合模块，配置为获取各第四功率谱数据中每个高度的等离子体频率，并通过样条插值法获取整个高度剖面的等离子体谱线；基于所述等离子体谱线，利用朗缪尔色散关系拟合得到电离层的电子密度剖面。

9.一种存储装置，其中存储有多条程序，其特征在于，所述程序应用由处理器加载并执行以实现权利要求1-7任一项所述的基于CLEAN算法的电离层电子密度反演方法。

10.一种处理装置，包括处理器、存储装置；处理器，适用于执行各条程序；存储装置，适用于存储多条程序；其特征在于，所述程序适用于由处理器加载并执行以实现权利要求1-7任一项所述的基于CLEAN算法的电离层电子密度反演方法。