CN105158777B - 用于无源测向定位的数据源生成方法 - Google Patents

Abstract

本发明提供了一种用于无源测向定位的数据源生成方法，其针对卫星无源测向定位试验成本高的特点，基于无源测向定位理论，包括初始设置、地球卫星运动学环境递推、观测数据生成等环节。本发明具有参数设置灵活、计算速度快、输出观测量准确性高的特点。本发明解决了星上定位算法设计过程中校验难和校验效率低的工程实际问题，能够为定位算法的验证工作提供大批量可靠数据源，以此得到的定位结果可作为测向定位算法准确性的考察依据。

Description

用于无源测向定位的数据源生成方法

技术领域

本发明涉及属卫星无源定位技术领域，具体地，涉及用于无源测向定位的数据源生成方法。

背景技术

卫星无源定位通常可分为单星体制和多星体制。单星体制是采用瞬时测向法对目标定位，基本方案主要有幅度法测向、相位法测向和时差法测向三种。无论使用哪一种测向原理，其目的都是为了得到目标在卫星本体坐标系下的方向。其中，相位法技术成熟，在工程上便于实现，容易被工程采用。

目前，测向定位算法的验证方法主要是在STK软件中设置仿真时间和卫星轨道和姿态参数，并假定一个目标位置，从而由STK计算出目标方向在卫星本体系下的表示，即β_x(目标方向与本体系X轴的夹角)和β_y(目标方向与本体系Y轴的夹角)，并生成报表。将STK中得出的“偏航角”和“俯仰角”代入定位算法得出目标位置经纬度，与STK中的输入值(目标位置经纬度)进行比较。

STK软件是在航空航天行业内被普遍认可的仿真软件，它具有图形界面和计算准确等优点。但对于测向定位算法的验证工作，STK软件亦存在一些局限性，可归纳为以下三点：

1、在一些研究中，根据相位测向原理，接收机在测角后输出载波相位的小数部分，再由硬件或软件解模糊以得到完整的载波相位，反算测角结果，继而进行定位解算，STK无法直接输入或输出载波相位信息，因此也无法分析载波频率对定位精度带来的影响；

2、在一些研究中，往往希望数据源产生形如卫星位置、速度、天线安装位置、载波相位误差等各种随机误差或按照误差模型产生误差，STK无法直接产生这些误差，即无法得到这些误差对测向误差的影响，继而难以系统地研究这些误差对定位误差的分布情况的影响；

3、在实际工程中，往往希望大批量处理一些输入，例如用于在一些给定工况下比较不同的定位算法的优劣，或是批量设置某一参数用以研究这一参数的变化对定位结果的影响分布，反复在STK中设置场景并生成报表的方法效率较低。

发明内容

针对现有技术中的缺陷，为了解决星上定位算法设计过程中校验难和校验效率低的工程实际问题，本发明提出了一种用于研究无源测向定位算法的数据源生成方法。

根据本发明提供的一种用于无源测向定位的数据源生成方法，包括如下步骤：

步骤1：初始设置环节；

步骤2：运动学环境构建与递推环节；

建立运动学环境，通过引入地球岁差、章动、自转和极移模型，计算惯性系和地固系的转换矩阵，继而得到目标随地球转动在惯性空间中的位置；通过轨道模型(二体轨道模型、J2轨道模型、甚至更加精确的轨道模型)根据时间递推卫星轨道，从而计算卫星及卫星载荷天线在空间中的位置以及卫星姿态；

步骤3：观测数据生成环节；

将目标相对于各载荷天线的位置转换至卫星本体坐标系，并以目标相对于各载荷天线在卫星本体坐标系中的位置计算目标至各天线的距离，继而得到无线电波从目标至各天线所需要的载波相位和各条基线相位差，并依据误差模型加入误差。

优选地，所述步骤1包括如下任一个或任多个步骤：

设置计时器时间步进、仿真时长、仿真开始时间；

设置地球参数，地球参数包括地球引力常数、赤道半径、极地半径、地球扁率；

设置卫星历元时刻的轨道六要素、基线长度并根据基线长度设置三个天线的安装位置；

设置载荷下视角大小；

设置卫星的属性，卫星的属性包括载波信号波长、基线长度、天线安装位置、载荷波束角；

设置目标的地理经纬度。

优选地，在所述步骤2中：

地球的转动由转换矩阵表达如下

HG＝[EP]·[ER]·[NR]·[PR]其中，为地固系下坐标，为J2000惯性系下坐标，[·]表示坐标转换矩阵，其中，HG为所要计算的自转矩阵，先后包括岁差矩阵PR、章动矩阵NR、地球自转矩阵ER和极移矩阵EP，分别由下列各式表达：

[PR]＝R_z(-z_A)R_y(θ_A)R_z(-ξ_A)

[NR]＝R_x(-Δε)R_y(Δθ)R_z(-Δμ)

[EP]＝R_y(-x_p)R_x(-y_p)

其中，R为绕指定轴的旋转矩阵，R_x表示绕x轴旋转的旋转矩阵，R_y表示绕y轴旋转的旋转矩阵，R_z表示绕z轴旋转的旋转矩阵，x_p表示极移x分量，y_p表示极移y分量；

岁差常数Z_A,θ_A,ξ_A由下列公式计算

章动量Δμ,Δθ,Δε，计算公式如下：

Δψ＝(A₀+A₁T)sin(α(T))

Δε＝(B₀+B₁T)cos(α(T))

A₀＝-17.1996″,A₁＝-0.0174.2″,B₀＝9.2025″,B₁＝0.00089″

α表示幅角，T表示儒略世纪数，Δψ表示黄经章动；

幅角计算公式为：

α＝125°02′40″.280

其中，Δμ和Δθ为

Δμ＝Δψcos(23°26′21″.448)

Δθ＝Δψsin(23°26′21″.448)

中的是格林威治平恒星时，

S_G＝18^h.697374558+879000^h.051336907T其中，h表示时角，T表示儒略世纪数，由下式计算；

JD(t)是当前UT1时刻的儒略日形式，UT1表示世界时；

平近点角由下式表示。

M(t)＝M(0)+n·t

式中，n表示卫星的轨道角速度，t表示仿真时间，M(t)表示t时刻的平近点角，M(0)表示0时刻的平近点角。

优选地，在所述步骤3的观测数据生成环节中，假设以卫星本体坐标系与轨道坐标系重合，其中，Z轴指向地心，Y轴沿轨道面负法向，X轴约束沿速方向，以卫星本体系原点建立球坐标系重合于轨道坐标系，辐射源目标T方向与Z轴的夹角为θ，方向角为其中，定义X轴正向为0，绕Z轴旋转为正，轨道坐标系原点O为地球的球心，G为星下点，点S₀，S₁，S₂为三副天线组成的等腰直角三角形天线阵，分别位于原点、X轴和Y轴上，设S₀S₁，S₀S₂为定位基线，基线的长度S₀S₁＝S₀S₂＝d，辐射源目标T的坐标为卫星高度为h，地球半径为R_e；其中，中的分别表示目标在球坐标下的距离、方位角和高度角的余角；

设天线S₀和天线S₁同目标的距离的差值为r-r₁，其中，r表示天线S₀到目标的距离，r₁表示天线S₁到目标的距离，那么，天线S₀和天线S₁的相位差由下式求得：

φ₁＝2π(r-r₁)/λ

式中，λ为载波信号波长；

设天线S₀和天线S₂同目标的距离的差值为r-r₂，其中，r表示天线S₀到目标的距离，r₂表示天线S₂到目标的距离，那么，天线S₀和天线S₂的相位差由下式求得：

φ₂＝2π(r-r₂)/λ

其中，φ₁表示天线S₀和天线S₁的相位差，φ₂表示天线S₀和天线S₂的相位差，r₂表示天线S₂到目标的距离。

优选地，在所述步骤3的观测数据生成环节中，还需将仿真生成的相位差信息转换至测角信息；

具体地，已知相位差并由相位差得到距离，有

其中，r表示天线S₀到目标的距离，r₁表示天线S₁到目标的距离，d表示基线长度，θ表示目标的高度角的余角，表示目标的方位角；

天线S₀和天线S₂的距离差r-r2为

天线S₀和天线S₁组成定位基线的相差测向方程为

天线S₀和天线S₂组成定位基线的相差测向方程为

式中，λ为载波信号波长，φ₁为天线S₀和天线S₁之间的相差，φ₂为天线S₀和天线S₂之间的相差；

解相差测向方程

由φ和θ求得目标在卫星本体系下的方向矢量(β_x,β_y,β_z)，β_x,β_y,β_z分别表示该方向矢量(β_x,β_y,β_z)与X，Y，Z轴的夹角，方向矢量(β_x,,β_y,β_z)和卫星速度、位置、姿态一同作为无源测向定位的数据源。

与现有技术相比，本发明具有如下的有益效果：

1、利用本发明能够为定位算法的验证工作提供大批量可靠数据源；

2、本发明具有参数设置灵活、计算速度快、输出观测量准确性高的特点；

3、本发明能够满足无源目标定位算法的研究需求。从仿真结果来看，该方法合理、可行，能够保证精度；

4、本发明将为无源目标测向定位算法的研究提供参考和依据；

5、本发明利用初始设置环节，可对卫星轨道、姿态、基线长度、载波频率、目标位置等多项参数进行直接设置，以考察不同参数对定位精度的影响；

6、本发明构造了地球椭球面模型，较之平坦大地和圆球形地球模型更加接近真实场景；

7、本发明中的运动学环境构建与递推环节具有卫星轨道和地球自转递推功能，以时间为线索模拟卫星的实时位置姿态和真实地球的实时姿态，完整模仿了地球、卫星、目标等对象的相互几何位置关系，展现了定位过程的本质，便于在没有卫星实测数据的情况下验证整个过顶过程的定位效果；

8、本发明利用观测数据生成环节，在测向数据生成模块中将测向角度转换成相位差，便于后续定位模块利用相位差还原至目标方位角和高度角，模拟了工程上的实际过程，也便于加入载波频率对定位精度的影响因素。

附图说明

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

图1是本发明定位算法流程框图；

图2是相差测向定位天线布局示意图；

图3为卫星本体坐标系下的目标方向示意图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

本发明提供了一种本发明针对卫星无源测向定位试验成本高的特点，基于无源测向定位理论，设计了一套用于研究无源测向定位算法的数据源生成方法，包括初始设置、地球卫星运动学环境递推、观测数据生成等环节。整个定位仿真系统结构如图1所示，填充有斜线阴影的部分是本发明的涵盖部分，涉及卫星运动模型、地球运动模型、坐标转换等内容，具有参数设置灵活、计算速度快、输出观测量准确性高的特点。本发明解决了星上定位算法设计过程中校验难和校验效率低的工程实际问题，能够为定位算法的验证工作提供大批量可靠数据源，以此得到的定位结果可作为测向定位算法准确性的考察依据。

具体地，本发明是一种脱离第三方软件的自主仿真方法，该方法包括：

初始设置环节，通过输入仿真时间、地球椭球模型、卫星初始轨道、卫星姿态、天线安装位置(至少三个天线组成两组基线方可定位)、测量误差模型(系统误差和随机误差)、目标经纬度、载波信号波长、仿真步长等；

运动学环境构建与递推环节，建立运动学环境，通过引入地球岁差、章动、自转和极移模型，精确计算惯性系和地固系的转换矩阵，继而得到目标随地球转动在惯性空间中的位置；通过二体或J2轨道模型(甚至更加精确的轨道模型)根据时间递推卫星轨道，从而计算卫星及其载荷天线在空间中的位置以及卫星姿态；

观测数据生成环节，将目标相对于各载荷天线的位置转换至卫星本体坐标系，并以此计算目标至各天线的距离，继而得到无线电波从目标至各天线所需要的载波相位和各条基线相位差，并依据误差模型加入误差。

如下对上述三个环节进行更为具体地描述。

1、初始设置环节

计时器时间步进、仿真时长、仿真开始时间；

设置地球参数，包括地球引力常数、赤道半径、极地半径、地球扁率等信息；

设置卫星历元时刻的轨道六要素，基线长度并根据基线长度设置三个天线的安装位置；设置载荷下视角大小；此外，载波信号波长、基线长度、天线安装位置、载荷波束角等信息均作为卫星的属性并予以设置；

设置目标的地理经纬度；

2、运动学环境构建与递推环节

环节中的关键步骤罗列如下：

a)地球的转动

地球的转动由转换矩阵表达如下

HG＝[EP]·[ER]·[NR]·[PR]其中，为地固系下坐标，为J2000惯性系下坐标。HG为所要计算的自转矩阵，先后包括岁差(PR)、章动(NR)、地球自转(ER)和极移(EP)矩阵。它们分别由下列各式表达。

[PR]＝R_z(-z_A)R_y(θ_A)R_z(-ξ_A)

[NR]＝R_x(-Δε)R_y(Δθ)R_z(-Δμ)

[EP]＝R_y(-x_p)R_x(-y_p)

其中，R为绕指定轴的旋转矩阵，Z_A,θ_A,ξ_A由下列公式计算

Δμ,Δθ,Δε等章动量计算公式如下：

Δψ＝(A₀+A₁t)sin(α(t))

Δε＝(B₀+B₁t)cos(α(t))

A₀＝-17.1996″,A₁＝-0.0174.2″,B₀＝9.2025″,B₀＝0.00089″

幅角计算公式为：

α＝125°02′40″.280

其中，Δμ和Δθ为

Δμ＝Δψcos(23°26′21″.448)

Δθ＝Δψsin(23°26′21″.448)

中的是格林威治平恒星时，

S_G＝18^h.697374558+879000^h.051336907T其中，

JD(t)是当前UT1时刻的儒略日形式。对于2007年7月1日的情况，有：

UT1＝UTC-0.1582305s

此外，R_y(-x_p)R_x(-y_p)中的参数已由初始配置给出。

b)卫星的运动

在给定时间下，结合卫星的初始轨道设置可得到卫星的平近点角，其它轨道要素不变。新的平近点角可由下式表示。

M(t)＝M(0)+n·t

式中，n表示卫星的轨道角速度，t表示仿真时间。

从轨道要素到空间位置速度的转换为标准流程，此处略。

3、观测数据生成环节

假设以卫星本体坐标系与轨道坐标系(Z轴指向地心，Y轴沿轨道面负法向，X轴约束沿速方向)重合，以本体系原点建立球坐标系(如图3所示)重合于轨道系，辐射源目标T方向与Z轴的夹角为θ，方向角为(X轴正向为0，绕Z轴旋转为正)，O为地球的球心，G为星下点，S₀，S₁，S₂为三副天线组成的等腰直角三角形天线阵，分别位于原点、X轴和Y轴上，为分析方便，设S₀S₁，S₀S₂为定位基线，基线的长度S₀S₁＝S₀S₂＝d，辐射源目标T的坐标为卫星高度为h，地球半径为R_e。

根据两维相差测向是利用信号传输到三个天线的相位差进行测向，相位差的产生是目标到不同天线的距离差引起的，因此首先研究距离差与信号到达角的关系。如所示，天线S₀和天线S₁的距离差为r-r₁，那么，两者的相位差可由下式求得。

φ₁＝2π(r-r₁)/λ

式中，λ为载波信号波长。同理可得，

φ₂＝2π(r-r₂)/λ

这样便得到了目标辐射源至不同卫星天线的载波相位差，亦可根据后续研究需要去除其整数部分以模仿整周模糊度的产生机理。

本环节还需将其转换至测角信息。对于这一情况，需已知相位差并由其得到距离，有

同理，天线S₀和天线S₂的距离差r-r2为

天线S₀和天线S₁组成定位基线的相差测向方程为

天线S₀和天线S₂组成定位基线的相差测向方程为

式中，λ为载波信号波长，φ_ι(i＝1,2)为相差。

解相差测向方程

由φ和θ可求得目标在卫星本体系下的方向(β_x,β_y,β_z)，β_x,β_y,β_z分别表示该方向矢量与X，Y，Z轴的夹角，此信息和卫星速度、位置、姿态一同作为无源测向定位算法的输入。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (3)

1.一种用于无源测向定位的数据源生成方法，其特征在于，包括如下步骤：

步骤1：初始设置环节；

步骤2：运动学环境构建与递推环节；

建立运动学环境，通过引入地球岁差、章动、自转和极移模型，计算惯性系和地固系的转换矩阵，继而得到目标随地球转动在惯性空间中的位置；通过轨道模型根据时间递推卫星轨道，从而计算卫星及卫星载荷天线在空间中的位置以及卫星姿态；

步骤3：观测数据生成环节；

将目标相对于各载荷天线的位置转换至卫星本体坐标系，并以目标相对于各载荷天线在卫星本体坐标系中的位置计算目标至各天线的距离，继而得到无线电波从目标至各天线所需要的载波相位和各条基线相位差，并依据误差模型加入误差；

在所述步骤3的观测数据生成环节中，假设以卫星本体坐标系与轨道坐标系重合，其中，Z轴指向地心，Y轴沿轨道面负法向，X轴约束沿速方向，以卫星本体系原点建立球坐标系重合于轨道坐标系，辐射源目标T方向与Z轴的夹角为θ，方向角为其中，定义X轴正向为0，绕Z轴旋转为正，轨道坐标系原点O为地球的球心，G为星下点，点S₀，S₁，S₂为三副天线组成的等腰直角三角形天线阵，分别位于原点、X轴和Y轴上，设S₀S₁，S₀S₂为定位基线，基线的长度S₀S₁＝S₀S₂＝d，辐射源目标T的坐标为卫星高度为h，地球半径为R_e；其中，中的r,θ,分别表示目标在球坐标下的距离、方位角和高度角的余角；

设天线S₀和天线S₁同目标的距离的差值为r-r₁，其中，r表示天线S₀到目标的距离，r₁表示天线S₁到目标的距离，那么，天线S₀和天线S₁的相位差由下式求得：

φ₁＝2π(r-r₁)/λ

式中，λ为载波信号波长；

设天线S₀和天线S₂同目标的距离的差值为r-r₂，其中，r表示天线S₀到目标的距离，r₂表示天线S₂到目标的距离，那么，天线S₀和天线S₂的相位差由下式求得：

φ₂＝2π(r-r₂)/λ

其中，φ₁表示天线S₀和天线S₁的相位差，φ₂表示天线S₀和天线S₂的相位差，r₂表示天线S₂到目标的距离；

在所述步骤3的观测数据生成环节中，还需将仿真生成的相位差信息转换至测角信息；

具体地，已知相位差并由相位差得到距离，有

其中，r表示天线S₀到目标的距离，r₁表示天线S₁到目标的距离，d表示基线长度，θ表示目标的高度角的余角，表示目标的方位角；

天线S₀和天线S₂的距离差r-r₂为

天线S₀和天线S₁组成定位基线的相差测向方程为

天线S₀和天线S₂组成定位基线的相差测向方程为

式中，λ为载波信号波长，φ₁为天线S₀和天线S₁之间的相差，φ₂为天线S₀和天线S₂之间的相差；

解相差测向方程

由φ和θ求得目标在卫星本体系下的方向矢量(β_x,β_y,β_z)，β_x,β_y,β_z分 别表示该方向矢量(β_x,β_y,β_z)与X，Y，Z轴的夹角，方向矢量(β_x,β_y,β_z)和卫星速度、位置、姿态一同作为无源测向定位的数据源。

2.根据权利要求1所述的用于无源测向定位的数据源生成方法，其特征在于，所述步骤1包括如下任一个或任多个步骤：

设置计时器时间步进、仿真时长、仿真开始时间；

设置地球参数，地球参数包括地球引力常数、赤道半径、极地半径、地球扁率；

设置卫星历元时刻的轨道六要素、基线长度并根据基线长度设置三个天线的安装位置；

设置载荷下视角大小；

设置卫星的属性，卫星的属性包括载波信号波长、基线长度、天线安装位置、载荷波束角；

设置目标的地理经纬度。

3.根据权利要求1所述的用于无源测向定位的数据源生成方法，其特征在于，在所述步骤2中：

地球的转动由转换矩阵表达如下

HG＝[EP]·[ER]·[NR]·[PR]

其中，为地固系下坐标，为J2000惯性系下坐标，[·]表示坐标转换矩阵，其中，HG为所要计算的自转矩阵，先后包括岁差矩阵PR、章动矩阵NR、地球自转矩阵ER和极移矩阵EP，分别由下列各式表达：

[PR]＝R_z(-z_A)R_y(θ_A)R_z(-ξ_A)

[NR]＝R_x(-Δε)R_y(Δθ)R_z(-Δμ)

[EP]＝R_y(-x_p)R_x(-y_p)

其中，R为绕指定轴的旋转矩阵，R_x表示绕x轴旋转的旋转矩阵，R_y表示绕y轴旋转的旋转矩阵，R_z表示绕z轴旋转的旋转矩阵，x_p表示极移x分量，y_p表示极移y分量；

岁差常数Z_A,θ_A,ξ_A由下列公式计算

章动量Δμ,Δθ,Δε，计算公式如下：

Δψ＝(A₀+A₁T)sin(α(T))

Δε＝(B₀+B₁T)cos(α(T))

A₀＝-17.1996″,A₁＝-0.0174.2″,B₀＝9.2025″,B₁＝0.00089″

α表示幅角，T表示儒略世纪数，Δψ表示黄经章动；

幅角计算公式为：

α＝125°02′40″.280

其中，Δμ和Δθ为

Δμ＝Δψcos(23°26′21″.448)

Δθ＝Δψsin(23°26′21″.448)

中的是格林威治平恒星时，

S_G＝18^h.697374558+879000^h.051336907T

其中，h表示时角，T表示儒略世纪数，由下式计算；

JD(t)是当前UT1时刻的儒略日形式，UT1表示世界时；

平近点角由下式表示；

M(t)＝M(0)+n·t

式中，n表示卫星的轨道角速度，t表示仿真时间，M(t)表示t时刻的平近点角，M(0)表示0时刻的平近点角。