一种地外天体探测器信号的频谱检测估计方法与流程

1.本发明属于航天器测控技术领域，具体涉及一种地外天体探测器信号的频谱检测估计方法。


背景技术：

2.地外天体探测器飞行的距离越来越远，其下行信号的强度越来越弱。目前，火星探测器到地面的最远距离达到4亿公里，而木星、土星及天王星等行星探测器到地面的距离更远，探测器到地面距离的平方与信号的衰减程度成正比，因此，随着距离的增大，探测器下行信号的强度显著减弱。此外，对于上述的行星探测任务，探测器将经历行星捕获的过程，部分探测器还将实现在行星上的着陆，探测器的速度变化快，由此造成下行信号的动态非常大。
3.传统方式中，地面站对探测器信号采用锁相环技术进行闭环的信号估计，完成对探测器信号的检测。当探测器信号强度相对较强、信号动态相对不大，可以有效的完成探测器信号的闭环检测。然而，对于上述探测器大动态、微弱信号，闭环信号估计方法不能够很好的适应，探测器信号检测估计难度大，信号失锁的概率显著增大。


技术实现要素：

4.针对现有技术存在的缺陷，本发明提供一种地外天体探测器信号的频谱检测估计方法，可有效解决上述问题。
5.本发明采用的技术方案如下：
6.本发明提供一种地外天体探测器信号的频谱检测估计方法，包括以下步骤：
7.步骤1，预先估计地面站需采集的探测器信号的频率范围和频率变化率范围；
8.步骤2，根据步骤1预先估计的地面站需采集的探测器信号的频率范围和频率变化率范围，对地面站的信号采集参数进行配置；
9.步骤3，地面站采用配置的信号采集参数，对探测器信号进行采集，得到原始采集信号；
10.步骤4，将原始采集信号，以时间长度t为周期，划分为多段采集信号单元；
11.对于每段采集信号单元d
t
，均执行步骤5-步骤8，搜索采集信号单元d
t
中的目标信号：
12.步骤5，将采集信号单元d
t
的频带宽度等间距连续划分为m个子频带，第m个子频带的采集信号表示为：d
t
(m),m＝1，2，3，
…
，m；
13.步骤6，将步骤1估计得到的频率变化率范围等间距连续划分为n个频率变化率估计值子区间；
14.对第m个子频带的采集信号d
t
(m)，分别计算其在每个频率变化率估计值子区间的频率与频率变化率二维补偿值，由此得到n个二维补偿后的采集信号g
t
(m,n)，其中，n＝1,2,3...,n；
15.步骤7，将第m个子频带的频带宽度进一步等间距连续划分为q个子频率区间；
16.对第m个子频带在第n个频率变化率估计值子区间的二维补偿后的采集信号g
t
(m,n)，分别计算其在每个子频率区间的傅里叶变换结果，由此得到q个傅里叶变换后的信号，即为三维采集信号g
t
(m,n,q)；
17.步骤8，因此，对于采集信号单元d
t
，在m个子频带、n个频率变化率估计值子区间和q个子频率区间进行精细搜索计算，得到m*n*q个三维采集信号g
t
(m,n,q)；
18.比较m*n*q个三维采集信号g
t
(m,n,q)的功率，功率峰值对应的子频带序号、频率变化率估计值子区间序号和子频率区间序号，即为搜索到的目标信号的位置，进而搜索到目标信号；
19.步骤9，根据搜索到的目标信号，估计其频谱检测估计值。
20.优选的，步骤1中，探测器信号的频率范围表示为[f
min
,f
max
]，频率变化率范围表示为[f
′
min
,f
′
max
]：
[0021]
其中：
[0022]fmin
为需采集的探测器信号的频率最小值估计值；
[0023]fmax
为需采集的探测器信号的频率最大值估计值；
[0024]f′
min
为需采集的探测器信号的频率变化率最小值估计值；
[0025]f′
max
为需采集的探测器信号的频率变化率最大值估计值；
[0026]
步骤2中，配置的地面站的信号采集参数包括地面站的频带宽度b0和中心频率fb，通过以下公式得到：
[0027]
b0＝f
max-f
min
[0028][0029]
优选的，步骤5中，采集信号单元d
t
的频带宽度，与地面站配置的频带宽度b0相等。
[0030]
优选的，步骤5中，每个子频带的采集信号d
t
(m)表示为：
[0031][0032]
其中：n0为噪声，s(f(t))为目标信号；
[0033]
含义为：
[0034]
在m个子频带中，只有一个子频带具有目标信号，假设在第z个子频带具有目标信号，此处z为未知值；其他m-1个子频带均为噪声；
[0035]
步骤6具体为：
[0036]
将步骤1估计得到的频率变化率范围[f
′
min
,f
′
max
]等间距连续划分为n个频率变化率估计值子区间，表示为：e(n)，n＝1,2,3...,n；因此，每一个频率变化率估计值子区间的宽度δf
′
＝e(n)/n；
[0037]
对第m个子频带的采集信号d
t
(m)，采用以下方式，进行频率与频率变化率二维补偿，得到二维补偿后的采集信号g
t
(m,n)：
[0038][0039]
其中：
[0040]
j代表虚数单位；
[0041]
fc(m)代表第m个子频带的下边界频率补偿值；
[0042]f′c(n)代表第n个频率变化率估计值子区间的下边界频率变化率补偿值；
[0043]
其中，fc(m)和f
′c(n)采用下式得到：
[0044]
fc(m)＝(m-1)δf+f
min
[0045]f′c(n)＝(n-1)δf
′
+f
′
min
[0046]
其中：δf＝b0/m，代表每一个子频带的频带宽度。
[0047]
优选的，步骤7具体为：
[0048]
对第m个子频带在第n个频率变化率估计值子区间的二维补偿后的采集信号g
t
(m,n)，采用以下公式，转化为三维采集信号g
t
(m,n,q)：
[0049][0050]
其中：
[0051]
fft表示傅里叶变换；
[0052]
含义为：
[0053]
将第m个子频带的频带宽度δf等间距连续划分为q个子频率区间，表示为：fr(q)，q＝1,2,...,q；对每个子频带的二维补偿后的采集信号g
t
(m,n)在每个子频率区间进行傅里叶变换，得到的傅里叶变换后的信号，即为三维采集信号g
t
(m,n,q)。
[0054]
优选的，步骤8中，三维采集信号g
t
(m,n,q)的功率表示为p
t
(m,n,q)，通过以下方式获得：
[0055][0056]
其中：
[0057]
含义为噪声部分傅里叶变换后结果；
[0058]
p[n0]代表n0的功率，含义为：噪声部分傅里叶变换后的功率；
[0059]
含义为目标信号s(f(t))部分傅里叶变换后结果；
[0060]
p[s(m,n,q)]代表s(m,n,q)的功率，含义为：目标信号s(f(t))部分傅里叶变换后的功率。
[0061]
优选的，步骤8中，功率峰值表示为p
t
(m
t
,n
t
,q
t
)；m
t
,n
t
,q
t
含义分别为：功率峰值所在的子频带序号、频率变化率估计值子区间序号和子频率区间序号，表示为：
[0062]
(m
t
,n
t
,q
t
)＝max[p
t
(m,n,q)]|
m∈[1,m],n∈[1,n],q∈[1,q]
[0063]
步骤9具体为：
[0064]
目标信号的频谱检测估计值，包括：目标信号的频率精确估计值频率变化率精确估计值以及信号功率噪声谱密度比精确估计值p/n0；
[0065]
采用下式，得到采集信号单元d
t
中目标信号的频率精确估计值和频率变化率精确估计值
[0066][0067][0068]
采用下式，得到采集信号单元d
t
中目标信号的信号功率噪声谱密度比精确估计值p/n0：
[0069]
p/n0＝p
t
(m
t
,n
t
,q
t
)/p[n0]
[0070]
p[n0]＝n/(δf/q)
[0071]
获得频谱曲线g
t
(m
t
,n
t
,q)，其中，q＝1,2,...,q，包括幅频曲线和相频曲线；
[0072]
其中：
[0073]
频谱曲线g
t
(m
t
,n
t
,q)代表横坐标为fr(q)，纵坐标为g
t
(m
t
,n
t
,q)形成的曲线。
[0074]
本发明提供的一种地外天体探测器信号的频谱检测估计方法具有以下优点：
[0075]
本发明将地面天线接收到的地外天体采集信号，按照频带范围划分多个子频带、多个子频率区间和多个频率变化率估计值子区间，从而开展快速的信号的分析处理，降低计算资源，实现信号快速分析处理。
附图说明
[0076]
图1为本发明提供的一种地外天体探测器信号的频谱检测估计方法的流程示意图；
[0077]
图2为本发明提供的频率和频率变化率二维补偿网格的示意图；
[0078]
图3为本发明提供的火星探测器信号频谱结果图。
具体实施方式
[0079]
为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0080]
本发明提供一种地外天体探测器信号的频谱检测估计方法，将接收到的采集信号频带划分多个子频带，每个子频带划分多个子频率区间；同时，将预估的频率变化率范围划分多个频率变化率估计值子区间，由此在子频带、子频率区间和频率变化率估计值子区间的三维空间中，开展快速的信号捕获及跟踪，恢复目标信号的频谱，并完成目标信号频率、频率变化率和信号功率噪声谱密度比的估计，实现地外天体探测器微弱信号频谱监测。本发明利用多个子频带、多个子频率区间和频率变化率估计值子区间的处理方法，可以实现极微弱信号的捕获及监测，特别适合于地外天体探测器信号频谱监测估计。
[0081]
本发明提供一种地外天体探测器信号的频谱检测估计方法，包括以下步骤：
[0082]
步骤1，预先估计地面站需采集的探测器信号的频率范围和频率变化率范围；
[0083]
本步骤中，探测器信号的频率范围表示为[f
min
,f
max
]，频率变化率范围表示为[f
′
min
,f
′
max
]；
[0084]
其中：
[0085]fmin
为需采集的探测器信号的频率最小值估计值；
[0086]fmax
为需采集的探测器信号的频率最大值估计值；
[0087]f′
min
为需采集的探测器信号的频率变化率最小值估计值；
[0088]f′
max
为需采集的探测器信号的频率变化率最大值估计值。
[0089]
由于地外天体目标的位置并不能精确已知，但可以根据先验信息进行初略估计，因此可初略确定地外天体目标(即探测器)和地面站之间的相对运动，从而实现地面站需采集的探测器信号的频率及频率变化率的初略估计。
[0090]
步骤2，根据步骤1预先估计的地面站需采集的探测器信号的频率范围和频率变化率范围，对地面站的信号采集参数进行配置；
[0091]
本步骤中，配置的地面站的信号采集参数包括两个，分别为地面站的频带宽度b0和中心频率fb，通过以下公式得到：
[0092]
b0＝f
max-f
min
[0093][0094]
步骤3，地面站采用配置的信号采集参数，对探测器信号进行采集，得到原始采集信号；
[0095]
步骤4，将原始采集信号，以时间长度t为周期，划分为多段采集信号单元；
[0096]
对于每段采集信号单元d
t
，均执行步骤5-步骤8，搜索采集信号单元d
t
中的目标信号：
[0097]
步骤5，将采集信号单元d
t
的频带宽度等间距连续划分为m个子频带，第m个子频带的采集信号表示为：d
t
(m),m＝1，2，3，
…
，m；其中，采集信号单元d
t
的频带宽度，与地面站配置的频带宽度b0相等。
[0098]
每个子频带的采集信号d
t
(m)表示为：
[0099][0100]
其中：n0为噪声，s(f(t))为目标信号；
[0101]
含义为：
[0102]
在m个子频带中，只有一个子频带具有目标信号，假设在第z个子频带具有目标信号，此处z为未知值；其他m-1个子频带均为噪声。
[0103]
步骤6，将步骤1估计得到的频率变化率范围等间距连续划分为n个频率变化率估计值子区间；
[0104]
具体的，将步骤1估计得到的频率变化率范围[f
′
min
,f
′
max
]等间距连续划分为n个频率变化率估计值子区间，表示为：e(n)，n＝1,2,3...,n；因此，每一个频率变化率估计值子区间的宽度δf
′
＝e(n)/n。
[0105]
对第m个子频带的采集信号d
t
(m)，分别计算其在每个频率变化率估计值子区间的频率与频率变化率二维补偿值，由此得到n个二维补偿后的采集信号g
t
(m,n)，其中，n＝1,2,3...,n；
[0106]
作为一种具体实现方式，对第m个子频带的采集信号d
t
(m)，采用以下方式，进行频率与频率变化率二维补偿，得到二维补偿后的采集信号g
t
(m,n)：
[0107][0108]
其中：
[0109]
j代表虚数单位；
[0110]
fc(m)代表第m个子频带的下边界频率补偿值；
[0111]f′c(n)代表第n个频率变化率估计值子区间的下边界频率变化率补偿值；
[0112]
其中，fc(m)和f
′c(n)采用下式得到：
[0113]
fc(m)＝(m-1)δf+f
min
[0114]f′c(n)＝(n-1)δf
′
+f
′
min
[0115]
其中：δf＝b0/m，代表每一个子频带的频带宽度。
[0116]
作为一种具体实现方式，可制作得到图2所示的频率与频率变化率二维补偿网格，每个网格预先存储对应的fc(m)和f
′c(n)的值，具体的，对于横轴第m个点(m＝1,2,3...,m)，纵轴第n个点(n＝1,2,3...,n)，频率及频率变化率的二维补偿矩阵为c(m,n)，c(m,n)＝(fc(m),fc′
(n))。
[0117]
例如，在图2的c(1,3)网格，存储fc(1)和f
′c(3)的值。因此，可直接查找图2，快速得到所需的fc(m)和f
′c(n)的值。
[0118]
步骤7，将第m个子频带的频带宽度进一步等间距连续划分为q个子频率区间；
[0119]
对第m个子频带在第n个频率变化率估计值子区间的二维补偿后的采集信号g
t
(m,n)，分别计算其在每个子频率区间的傅里叶变换结果，由此得到q个傅里叶变换后的信号，即为三维采集信号g
t
(m,n,q)；
[0120]
本步骤中，对第m个子频带在第n个频率变化率估计值子区间的二维补偿后的采集信号g
t
(m,n)，采用以下公式，转化为三维采集信号g
t
(m,n,q)：
[0121][0122]
其中：
[0123]
fft表示傅里叶变换；
[0124]
含义为：
[0125]
将第m个子频带的频带宽度δf等间距连续划分为q个子频率区间，表示为：fr(q)，q＝1,2,...,q；对每个子频带的二维补偿后的采集信号g
t
(m,n)在每个子频率区间进行傅里叶变换，得到的傅里叶变换后的信号，即为三维采集信号g
t
(m,n,q)。
[0126]
步骤8，因此，对于采集信号单元d
t
，在m个子频带、n个频率变化率估计值子区间和q个子频率区间进行精细搜索计算，得到m*n*q个三维采集信号g
t
(m,n,q)；
[0127]
比较m*n*q个三维采集信号g
t
(m,n,q)的功率，功率峰值对应的子频带序号、频率变化率估计值子区间序号和子频率区间序号，即为搜索到的目标信号的位置，进而搜索到目标信号；
[0128]
本步骤中，三维采集信号g
t
(m,n,q)的功率表示为p
t
(m,n,q)，通过以下方式获得：
[0129][0130]
其中：
[0131]
含义为噪声部分傅里叶变换后结果；
[0132]
p[n0]代表n0的功率，含义为：噪声部分傅里叶变换后的功率；
[0133]
含义为目标信号s(f(t))部分傅里叶变换后结果；
[0134]
p[s(m,n,q)]代表s(m,n,q)的功率，含义为：目标信号s(f(t))部分傅里叶变换后的功率。
[0135]
本步骤中，功率峰值表示为p
t
(m
t
,n
t
,q
t
)；m
t
,n
t
,q
t
含义分别为：功率峰值所在的子频带序号、频率变化率估计值子区间序号和子频率区间序号，表示为：
[0136]
(m
t
,n
t
,q
t
)＝max[p
t
(m,n,q)]|
m∈[1,m],n∈[1,n],q∈[1,q]
[0137]
步骤9，根据搜索到的目标信号，估计其频谱检测估计值。
[0138]
本步骤中，目标信号的频谱检测估计值，包括：目标信号的频率精确估计值频
率变化率精确估计值以及信号功率噪声谱密度比精确估计值p/n0；
[0139]
采用下式，得到采集信号单元d
t
中目标信号的频率精确估计值和频率变化率精确估计值
[0140][0141][0142]
采用下式，得到采集信号单元d
t
中目标信号的信号功率噪声谱密度比精确估计值p/n0：
[0143]
p/n0＝p
t
(m
t
,n
t
,q
t
)/p[n0]
[0144]
p[n0]＝n/(δf/q)
[0145]
获得频谱曲线g
t
(m
t
,n
t
,q)，其中，q＝1,2,...,q，包括幅频曲线和相频曲线；
[0146]
其中：
[0147]
频谱曲线g
t
(m
t
,n
t
,q)代表横坐标为fr(q)，纵坐标为g
t
(m
t
,n
t
,q)形成的曲线。
[0148]
以下对本发明的计算复杂度进行分析：
[0149]
地面站配置的频带宽度为b0，计算的复杂度可以表示为o1：
[0150]
o1∝
(2b0)2×m×n[0151]
共划分为q个子频率区间后，计算的复杂度可以表示为o2：
[0152]
o2∝
[(2b0/m)2]
×m×
n＝(2b0)2×
n/m
[0153]
可以看出，o2＝o1/m2，采用的子频带数目(即m)越多，则本发明计算复杂度越低，需要的时间有效的缩短。
[0154]
仿真分析试验
[0155]
条件：
[0156]
1、地外天体目标模拟信号按照如下三种实际工况，假设标称下行频率f0＝8431mhz：
[0157]
工况一，频率f0+20khz、频率变化率-100hz/s，信号功率噪声谱密度比p/n0＝20dbhz。
[0158]
工况二，频率f0+20khz、频率变化率100hz/s，信号功率噪声谱密度比p/n0＝20dbhz。
[0159]
工况三，频率f0+20khz、频率变化率-10hz/s，信号功率噪声谱密度比p/n0＝20dbhz。
[0160]
2、数据生成：地外天体目标模拟信号经采集转换和记录设备(参数配置：采集记录带宽2mhz，量化位数16bit)，转换为原始数据，通过文件形式加载。
[0161]
3、信号电平、多普勒及动态预报文件在上述1的基础上增加误差，频率误差
±
5khz、频率变化率误差
±
20hz/s，信号功率噪声谱密度比误差
±
5dbhz。
[0162]
4、子通道数目(即q)为20，频率和频率变化率补偿网格数目取50，10。
[0163]
采用本发明方法，进行频率、频率变化率以及信号功率噪声谱密度比的估计，结果如下：
[0164]
工况一，频率f0+20.001khz、频率变化率-99.017hz/s，p/n0＝18.02dbhz。
[0165]
工况二，频率f0+20.000khz、频率勒变化率98.68hz/s，p/n0＝18.10dbhz。
[0166]
工况三，频率f0+19.999khz、频率勒变化率-9.7hz/s，p/n0＝17.73dbhz。
[0167]
由此可见，利用本发明方法可以有效的对地外天体探测器信号频谱进行监测估计，根据上述三种工况的结果，频率的估计精度达到1hz，频率变化率估计精度达到1hz/s，信号功率噪声谱密度比估计精度达到1.5dbhz。
[0168]
计算复杂度分析：
[0169]
传统的分析方法，o1∝
(2
×2×
106)2×
20
×
10＝3.2
×
10
15
，本发明方法，o2＝(2
×2×
106/20)2×
20
×
10＝8
×
10
12
，可以看出，本发明的计算复杂度，远小于传统方法。
[0170]
实际应用分析
[0171]
条件：
[0172]
1、火星着陆器信号理论值：频率f
0-362.1khz、频率变化率1hz/s，信号功率噪声谱密度比p/n0＝30dbhz。
[0173]
2、数据生成：火星着陆器信号由地面天线接收，经采集转换和记录设备(参数配置：采集记录带宽2mhz，量化位数16bit)，转换为原始数据，通过文件形式加载。
[0174]
3、采用本发明方法，进行频率、频率变化率以及信号功率噪声谱密度比的估计，结果如下：
[0175]
火星着陆器信号估计结果：频率f
0-362.103khz、频率变化率为0.95hz/s，p/n0＝31.29dbhz。估计到的火星着陆器信号频谱见图3。
[0176]
因此，应用本发明方法，成功监测并估计火星着陆器信号频谱，信号的频率与理论值的偏差为3hz，频率变化率与理论值的偏差为0.05hz/s，信号功率噪声谱密度比与基带测试结果的偏差为1.29dbhz。由此证明本发明估计值的精度非常高。
[0177]
本发明提供的一种地外天体探测器信号的频谱检测估计方法，具有以下优点：
[0178]
1、本发明将地面天线接收到的地外天体采集信号，按照频带范围划分多个子频带、多个子频率区间和多个频率变化率估计值子区间，从而开展快速的信号的分析处理，降低计算资源，实现信号快速分析处理。
[0179]
2、本发明可实现信号高精度估计，目标信号频率的估计精度达到1hz、频率变化率的估计精度达到1hz/s，、信号噪声谱密度比估计精度达到1.5dbhz。
[0180]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。