异地天线阵中基于预测模型的时延和时延率快速搜索方法

摘要

本发明公开了一种异地天线阵中基于预测模型的时延和时延率快速搜索方法，该方法包括以下步骤：预测值计算，预测模型拟合和残差搜索。本发明的时延和时延率快速搜索方法引入了时延和时延率的预测模型，对时延和时延率的搜索范围进行了有效限定，极大的提高了互相关处理的搜索效率，从而可以快速得到异地天线阵中各个天线的时延和时延率值，为后续的天线阵数据处理提供支持。

说明书

技术领域

本发明涉及深空探测技术领域，尤其涉及一种异地天线阵中基于预测模型的信号时延和时延率的快速搜索方法。

背景技术

随着人类深空探测活动的不断拓展，航天器的飞行距离越来越远，航天器发射的信号也随距离增加越来越弱。由于单面天线的口径和接收机系统的技术指标存在工程极限，需要借助天线组阵技术来达到接收信号所需要的最低信噪比。天线组阵的目的简要概括为：在相同的传输码率下，天线组阵技术能够增加地面对航天器的接收距离；在相同接收距离下，天线组阵技术能够提高地面与航天器之间的传输码速率。

异地天线阵指的利用现有的分布在不同区域的大口径天线同时观测航天器，然后通过数字信号处理技术将来自各个天线的信号进行记录、合成，以达到提高接收信号信噪比的目的。这种天线组阵技术的优势在于可以利用现有的天线及接收设备，只需要增加后端数据处理设备就可以实现天线组阵，研制周期短和成本低。

随着火星探测工程的立项，我国的深空探测活动也开始了由月球向火星的拓展。由于火星距离地球十分遥远，信号在空间传输过程中衰减巨大，相对于月球，火星距离增加了1000倍，信号空间衰减增加约60dB。地面接收到的探测器信号非常微弱，如采用现有的单天线接收，无法满足探测器的科学数据下行需求。如果新建口径更大的单天线来进行数据接收，由于受到单天线工程极限、研制进度、研制经费等诸多条件限制，使得整个项目存在极大的研制风险。

发明内容

(一)要解决的技术问题

有鉴于此，本发明的主要目的在于提供一种异地天线阵中基于预测模型的时延和时延率快速搜索方法，以快速得到异地天线阵中各个天线的时延和时延率值。

(二)技术方案

本发明提供一种基于预测模型的时延和时延率快速搜索方法，包括如下步骤：步骤1：根据天线位置和观测目标位置，初步计算出观测目标到天线之间的时延预测值；步骤2：对计算得到的时延预测值进行最小二乘模型拟合，得到时间连续的时延预测模型；步骤3：利用得到的时延预测模型对天线阵观测数据进行初步的时延和时延率补偿，然后使用互相关方法对天线数据进行互相关处理，对时延真值和时延预测值之间的残差进行搜索。

上述方案中，所述步骤1包括：步骤101：读取观测目标跟踪文件列表(t，ra，dec)和异地阵中各个天线地理坐标(x_i，y_i，z_i)，其中(t，ra，dec)表示在t时刻，观测目标的赤经和赤纬值；(x_i，y_i，z_i)表示地球空间直角坐标系中的每个天线位置，下标i用来区分不同的天线，阵中有N个天线，则有N组(x_i，y_i，z_i)，i取0，1...N-1；步骤102：根据经典天体时延计算方法计算观测目标到天线以及参考点的几何时延预测值。

上述方案中，步骤102中所述经典天体时延计算方法包括：步骤A1：输入测站坐标和探测器位置；步骤A2：修正由潮汐、板块运动等造成的测站本地坐标误差；步骤A3：将地心地球坐标系转换到地心天球坐标系；步骤A4：通过洛伦兹变换将地心天球坐标系变换到太阳系质心天球坐标系，在太阳系质心天球坐标系中计算几何时延，并修正太阳系重力吸引造成的信号传输路径弯曲误差；步骤A5：通过洛伦兹变换将太阳系质心天球坐标系变换回地心天球坐标系，得到所求的几何时延预测值。

上述方案中，所述步骤2包括：步骤201：将得到的每个天线关于地心参考点的几何时延预测值τ_i(t_j)读入，根据观测目标的运动特性，对上述预测值列表分成M分段，M为自然数；步骤202：使用最小二乘法，对上述一段预测值进行多项式拟合，得到多项式模型y；步骤203：评估多项式模型y拟合精度，拟合精度满足要求，得到时延预测模型。

上述方案中，步骤202中所述最小二乘法采用5阶多项式进行最小二乘拟合。

上述方案中，所述步骤203包括：根据公式计算得出一个标准精度，f为观测信号频率，单位为Hz；拟合误差是通过公式求得，式中y_j代表多项式模型y当X＝j时刻值，τ₀(j)，τ₁(j)对应几何时延计算结果表中密云站和昆明站j时刻的几何时延计算值，n表示每一分段预测值的个数；拟合误差RMS需要小于标准精度；如果精度满足要求，则该组多项式系数可以作为该天线的时延预测模型，如果拟合精度不满足要求，返回步骤201，调整分段或者增加多项式阶数，重新进行拟合，直到拟合精度满足要求。

上述方案中，步骤3中所述时延真值和时延预测值之间的残差包括残余时延和残余时延率。

上述方案中，所述步骤3包括：步骤301：选取天线阵中口径最大的天线作为参考天线，选取步骤203中得到天线阵的1个时延预测模型值减去参考天线的预测模型值，得到其他天线相对于参考天线的时延差值预测模型y(i，t_j)；步骤302：读入参考天线观测数据O_R(t_j)及阵中其他一个天线观测数据O_i(t_j)，去除观测数据的直流分量，得到O^a_R(t_j)，O^a_i(t_j)；步骤303：从O^a_R(t_j)数据中取2ⁿ个数据，利用步骤301得到的时延差值预测模型y(i，t_j)对观测数据O^a_i(t_j)进行时延补偿和时延率补偿，得到补偿后的2ⁿ个数据[O^a_i(t_j-y(i，t_j))]；步骤304：分别对O^a_R(t_j)，[O^a_i(t_j-y(i，t_j))]进行快速傅立叶变换，并进行互相关，得到互相关频谱fft(O_R^a(t_j))×fft^*([O^a_i(t_j-y(i，t_j))])；求互相关频谱fft(O_R^a(t_j))×fft^*([O^a_i(t_j-y(i，t_j))])的相位值，然后对相位值解卷绕，然后利用最小二乘方法对相位曲线进行一次线性函数拟合，得到这个2ⁿ个数据的残余时延Δτ¹；步骤305：取下2ⁿ个数据，同样进行上述304步骤计算，得到下2ⁿ个数据的残余时延Δτ²，计算得到时延率步骤306：将上述搜索得到的残余时延Δτ¹和时延率加到对应的时延差值预测模型y(i，t_j)中，即可得到可供后续天线组阵数据处理使用的真实时延模型。

(三)有益效果

本发明提供的时延和时延率快速搜索方法，所具有的积极效果在于：

引入了时延和时延率的预测模型，对时延和时延率的搜索范围了有效限定，极大地提高了互相关处理的搜索效率，从而可以快速得到异地天线阵中各个天线的时延和时延率值，为后续的天线阵数据处理提供支持。

附图说明

图1是异地天线组阵信号接收示意图；

图2是本发明异地天线阵中基于预测模型的时延和时延率快速搜索方法流程图；

图3是经典天体时延计算流程图；

图4是本发明的一种具体实施例的模型计算值与多项式拟合值对比；

图5是本发明的一种具体实施例的存在时延的互相关系数；

图6是本发明的一种具体实施例的存在残余时延率的互相关干涉条纹；

图7是本发明的一种具体实施例的残余时延补偿后的互相关系数；

图8是本发明的一种具体实施例的条纹旋转后的干涉条纹；

图9是本发明的一种具体实施例的分数时延调整前互相关相位谱；

图10是本发明的一种具体实施例的分数时延调整后的互相关相位谱。

具体实施方式

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

为了满足火星探测任务要求，在火星探测任务中拟采用异地组阵接收的方式来满足数据接收的需求。与本地组阵相比，由于接收探测器下行信号的路径不同，不同地点的天线接收到的信号会产生较大的时延和时延率，时延和时延率的修正对信号合成效率、合成信号的幅度有重要影响。鉴于此，本发明提供了一种异地天线阵中基于预测模型的时延和时延率快速搜索方法。

在本发明实施例中，异地天线阵由N个天线组成，最后可以简化为多个2天线阵进行数据处理，在此以一个2天线阵为示意性实例用于解释本发明，但不作为对本发明的限定，其中一个天线位于北京，另一个天线位于昆明，如图1所示，图1是异地天线组阵信号接收示意图。

图1中τ代表航天器到达两个天线之间的几何时延(时间差)，代表航天器到达两个天线之间的几何时延率(时间差变化速率)，g代表航天器到达两个天线之间的几何路程差，c代表光速。从图1中可以看出信号在空间传播的过程中，由于地球与航天器相对运动、大气折射、电离层扰动及接收信道等多方面因素，其到达两个地面站的时延值及时延率是不断变化的。为了能够将异地天线阵中的数据进行信号合成，必须要对时延值及时延率进行补偿。时延主要由几何时延、时钟引起的时延及大气引起的时延三部分构成。几何时延是整个时延中占比最大的一部分，主要由观测站与观测目标所处的几何位置引起。在知道观测站与观测目标所处的几何位置的情况下，可以通过计算得到几何时延的预测值。由于有了预测值的限定，后续的时延和时延率的残差搜索效率将大大提高。

图2是本发明异地天线阵中基于预测模型的时延和时延率快速搜索方法流程图，该方法包括以下步骤：预测值计算，预测模型拟合和残差搜索，其中：

步骤1：预测值计算是根据天线位置和观测目标位置，初步计算出观测目标到天线之间的几何时延预测值，该步骤具体包括：

步骤101：读取观测目标跟踪文件列表(t，ra，dec)和异地阵中各个天线地理坐标(x_i，y_i，z_i)，其中(t，ra，dec)表示在t时刻，观测目标的赤经和赤纬值；(x_i，y_i，z_i)表示地球空间直角坐标系中的每个天线位置，下标i用来区分不同的天线，阵中有N个天线，则有N组(x_i，y_i，z_i)，i取0，1...N-1。

所述观测目标跟踪文件是一组由观测目标位置数据组成列表。如果观测目标是卫星，那么观测目标跟踪文件就是由卫星轨道根数计算得到的一组卫星位置数据。为了方便理解观测目标跟踪文件数据格式，下面以嫦娥三号卫星跟踪文件为例给出观测目标跟踪文件数据格式：

YearMMUUHHMMSSDistanceRADEC201551000037460220.0849-15.1805201551001037460020.0856-15.188201551002037459720.0862-15.1771201551003037459520.0869-15.155201551004037459320.0875-15.1738201551005037459120.0882-15.1721201551006037458920.0888-15.1704201551007037458720.0895-15.1688201551008037458520.0901-15.1671201551009037458320.0908-15.1654

表格中各个参数的含义是：YearMMDDHHMMSS表示时间标识，格式为年月日时分秒；Distance表示航天器到地球地心的距离；RA(Right ascension)表示观测目标的赤经值；DEC(Declination)表示该观测目标的赤纬值。

读取异地天线阵中各个天线的地理坐标，此时，N＝2，有2个天线，分别位于，昆明和北京，如下表所示：

NameX(m)Y(m)Z(m)Km-1281152.9395640864.4072682653.403Bj-2201304.824324789.0454125367.718

表中Name异地天线阵中各个天线代号，Km表示昆明站40米天线，Bj表示北京密云站50米天线；X(m)Y(m)(m)分别表示每个天线在地心直角坐标系中的X轴，Y轴，Z轴数值，单位为米。

步骤102：根据经典天体时延计算方法计算观测目标到天线以及参考点的几何时延预测值。

在国际天球参考架(ICRF)定义的范畴内，以地心为参考点，计算观测目标到天线以及参考点的时延预测值列表τ_i(t_j)，其中t_j表示时间，下标i用来区分不同的天线，阵中有N个天线，则有N列几何时延预测值τ_i(t_j)，i取0，1...N-1。

在天球参考系中，天球参考系被近似为惯性坐标系。在该坐标系中，默认的原点是太阳系质心，赤道面定义为历元2000(J2000)赤道和春分点确定的平面。在地球坐标系中，采用国际习用原点(CIO，Conventional International Origin)确定的主轴方向和相应的赤道面。利用探测器的轨道数据和天线之间的基线矢量，进行几何时延的计算过程如图3所示。图3是经典天体时延计算流程图，首先输入测站坐标和探测器位置，然后修正由潮汐、板块运动等造成的测站本地坐标误差，将地心地球坐标系转换到地心天球坐标系(修正岁差、章动、摄动造成的误差)，随后通过洛伦兹变换将地心天球坐标系变换到太阳系质心天球坐标系，在太阳系质心天球坐标系中计算几何时延，并修正太阳系重力吸引造成的信号传输路径弯曲误差，最后通过洛伦兹变换将太阳系质心天球坐标系变换回地心天球坐标系，得到所求的几何时延。

根据上述计算，以地心为参考点，计算观测目标到天线以及参考点的几何时延预测值列表τ_i(t_j)，其中t_j表示时间，下标i用来区分不同的天线，阵中有2个天线，则有2列几何时延预测值τ_i(t_j)，i取0，1...N-1。

为了方便理解时延预测值的计算过程，下面给出了一组以嫦娥三号着陆器为观测目标(时间段为2015年5月10日00时2分0秒-59秒)，昆明站40米天线和密云站50米天线作为异地组阵天线，几何时延预测值计算结果(地心为参考点的几何时延值)。对几何时延预测值计算结果按照60秒进行分段，此为其中一段。

几何时延计算结果表

步骤2：预测模型拟合，是利用最小二乘方法，对上述计算得到的时延预测值进行最小二乘模型拟合，得到时间连续的几何时延预测模型；该步骤具体包括：

步骤201：将上述步骤1计算得到每个天线关于地心参考点的几何时延预测值列表τ_i(t_j)读入，根据观测目标的运动特性，对上述预测值列表分成M分段，M为自然数。对应到上述计算实例中，将表中密云站和昆明站时延预测值读入，并且以60秒时间进行分段。

步骤202：使用最小二乘法，对上述一段预测值进行多项式拟合，得到多项式模型y。

对应到上述计算实例中，结合探测器运动特性以及经过多次拟合精度评估，采用5阶多项式进行最小二乘拟合，能够满足后续计算要求。得到的多项式模型为：

Y＝6212.21513+0.110956058X+(-1.236351374e-05X²)+(-3.101422866e-10X³)+(3.286324942e-12X⁴)+(-1.306814583e-14X⁵)。

步骤203：评估多项式模型y拟合精度，拟合精度满足要求，得到时延预测模型。

拟合误差RMS需要小于f为观测信号频率，单位为Hz。如果精度满足要求，则该组多项式系数可以作为该天线的预测模型，如果拟合精度不满足要求，返回步骤201，调整分段或者增加多项式阶数，重新进行拟合，直到拟合精度满足要求，该M组多项式系数可以作为该天线的预测模型。然后依次对所有的天线进行上述模型拟合工作，得到整个异地天线阵的预测模型。

对应到上述计算实例中，多项式模型y＝6212.21513+0.110956058X+(-1.236351374e-05X²)+(-3.101422866e-10X³)+(3.286324942e-12X⁴)+(-1.306814583e-14X⁵)拟合误差RMS＝2.862e-07us，为均方根误差。该拟合误差是通过公式求得，式中y_j代表多项式模型y当X＝j时刻值，τ₀(j)，τ₁(j)对应几何时延计算结果表中密云站和昆明站j时刻的几何时延计算值。

如图4所示为本发明的一种具体实施例的模型计算值与多项式拟合值对比，横轴为时间，纵轴为时延，结合实施例说明了多项式拟合值与模型计算值误差很小，可以使用拟合得到多项式模型作为时延模型。

步骤3：残差搜索，是利用上述预测模型对天线阵观测数据进行初步的时延和时延率补偿，然后使用互相关方法对天线数据进行互相关处理，对时延真值和预测值之间的残差进行搜索；该步骤具体包括：

步骤301：选取天线阵中口径最大的天线作为参考天线，选取步骤203中得到天线阵的1个时延预测模型值减去参考天线的预测模型值，得到其他天线相对于参考天线的时延差值预测模型y(i，t_j)。

对应到上述计算实例中，选取密云站50米口径天线作为参考天线，

步骤302：读入参考天线观测数据O_R(t_j)及阵中其他一个天线观测数据O_i(t_j)，去除观测数据的直流分量，O^a_R(t_j)，O^a_i(t_j)得到。

对应到上述计算实例中，已参考天线为例，分别读取一段密云站天线的记录数据O_R(t_j)，求取这段数据的平均值(直流分量)然后将这段记录数据O_R(t_j)减去得到O^a_R(t_j)，其他天线的处理依次类推。

步骤303：从O^a_R(t_j)数据中取2ⁿ个数据，利用步骤301得到的时延差值预测模型y(i，t_j)对观测数据O^a_i(t_j)进行时延补偿和时延率补偿，得到补偿后的2ⁿ个数据[O^a_i(t_j-y(i，t_j))]。

对应到上述计算实例中，假设先从O^a_R(t_j)中取256个数据，然后取预测模型tj＝0，得到该时刻的时延预测值τ＝6.21221513e-03s，时延率预测值由于天线记录数据的采样速率是1e+07，这每个数据的采样时间Ts＝1e-07s，那个时延补偿数据位数＝τ/1e-7＝62122，对另一个天线数据进行时延补偿，就需要相同的记录时刻往后跳过时延补偿数据位数，读取数据256个数据，这个数据才与参考天线的的256个数据整数位对齐，对齐误差小于一个采样Ts。

然后，再对这段数据进行多普勒频移补偿。

步骤304：分别对O^a_R(t_j)，[O^a_i(t_j-y(i，t_j))]进行快速傅立叶变换，并进行互相关，得到互相关频谱fft(O_R^a(t_j))×fft^*([O^a_i(t_j-y(i，t_j))])；求互相关频谱fft(O_R^a(t_j))×fft^*([O^a_i(t_j-y(i，t_j))])的相位值，然后对相位值解卷绕，然后利用最小二乘方法对相位曲线进行一次线性函数拟合，得到这个2ⁿ个数据的残余时延Δτ¹。

步骤305：取下2ⁿ个数据，同样进行上述304步骤计算，得到下2n个数据的残余时延Δτ²，计算得到时延率

步骤306：将上述搜索得到的残余时延Δτ¹和时延率加到对应的几何时延预测y(i，t_j)中，即可得到可供后续天线组阵数据处理使用的真实时延模型。

对应到本实施例中Δτ¹＝5e-7s，那么真实时延模型

对剩下M-1组数据重复步骤201至步骤306的运算，即可得到该天线所有观测数据对应的真实时延模型。

对天线阵中剩余的其他天线进行上述步骤101至步骤306所有操作，就可以得到所有天线观测数据对应的真实时延模型。

对应到上述计算实例中，对两站数据进行初步时延和时延率校正，校正完成后对两站数据进行互相关处理，得到互相关系数、相关条纹如图5和图6所示。

如图5所示为本发明的一种具体实施例的存在时延的互相关系数图，横轴为时间，纵轴为互相关系数幅度，结合实施例说明了由于残余时延存在，互相关系数的最大值没有出现时间零点，存在一定偏移。

如图6所示为本发明的一种具体实施例的存在残余时延率的互相关干涉条纹图，横轴为时间，纵轴为互相关干涉条纹相位，结合实施例说明了由于残余时延率存在，互相关干涉条纹相位斜率不为零。

从互相关结果可以看出，由于模型的估算值存在一定的误差，相关系数和相关频谱的相位与理想的自相关谱存在一定的偏差，这个偏差就是模型时延值和真实时延之间的残差。通过进一步的时延和时延率估计得到残余时延为5e-7s，时延率4.11e-11s/s。

如图7所示为本发明的一种具体实施例的残余时延补偿后的互相关系数图，横轴为时间，纵轴为互相关系数幅度，结合实施例说明了由于对残余时延进行了补偿，互相关系数的最大值出现在时间零点。

如图8所示为残余时延率补偿的互相关干涉条纹图，横轴为时问，纵轴为互相关干涉条纹相位，结合实施例说明了由于对残余时延率进行了补偿，互相关干涉条纹相位斜率基本为零。

如图9所示为分数时延调整前互相关相位频谱图，横轴为频率，纵轴为相位，结合实施例说明了由于分数位时延的存在，互相关相位频谱图中相位斜率不为零。

如图10所示为分数时延调整后的互相关相位频谱图，横轴为频率，纵轴为相位，结合实施例说明了由于对分数位时延进行补偿，互相关相位频谱图中相位斜率基本为零。

从互相关结果可以看出，互相关系数、条纹和相位谱与理想的自相关谱相符合。

本发明的时延和时延率快速搜索方法引入了时延和时延率的预测模型，对时延和时延率的搜索范围进行了有效限定，极大的提高了互相关处理的搜索效率，从而可以快速得到异地天线阵中各个天线的时延和时延率值，为后续的天线阵数据处理提供支持。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。