CN115426027B - 一种基于干扰分布的受扰系统地球站位置计算方法 - Google Patents

Abstract

本发明公开了一种基于干扰分布的受扰系统地球站位置计算方法，该方法包括：步骤1)获取干扰系统卫星S_I、干扰系统地球站E_I和受扰系统卫星S的位置数据，以及干扰地球站E_I外设定的区域范围[0,length]，令H＝length；步骤2)计算在干扰系统地球站E_I外设定距离为H时，干扰最小的受扰地球站E_x位置、最小干扰值I/N，以及受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角α_H；步骤3)判断步骤2)得到的最小干扰值I/N是否满足限值要求，判断为是，令α_max＝α_H，采用二分法计算I/N满足限值要求的夹角最小值α_min，令α＝α_min，转至步骤4)；判断为否，无满足的受扰地球站位置，转至步骤5)；步骤4)根据α计算受扰系统地球站E_x的位置坐标；步骤5)结束。

Description

一种基于干扰分布的受扰系统地球站位置计算方法

技术领域

本发明属于空间干扰规避技术领域，尤其涉及一种基于干扰分布的受扰系统地球站位置计算方法。

背景技术

随着商业航天的蓬勃发展，Ka、Q/V频段资源日趋紧张，各国卫星操作者纷纷将目光投向更高频段。长期在轨运行的对地静止轨道(Geostationary Orbit，GSO)卫星，信息接收、管理、分发须依赖地球站，地球站之间的干扰评估是选址的重要因素之一，开展GSO卫星系统地球站址间干扰分析对W频段共享具有重要意义。

目前针对GSO系统卫星和地球站布设中的同频干扰问题，主要面向X、Ka等较低频段，采用栅格化方法遍历计算干扰值以评估。董苏惠等人提出基于干扰函数极值的评估方法，通过建立随机分布的地球站集合，分析不同轨位间隔的干扰、受扰卫星系统间的干扰分布；李建欣等人基于遍历方位角的方法，计算工作在8450-8500MHz频段的卫星固定业务和近地空间研究业务地球站间干扰保护距离范围。针对W频段的频谱兼容性分析主要面向通信链路传播衰减影响，Riva分析了W频段卫星通信中降雨、云雾等天气条件引起的链路衰减影响；Cuervo分析了75GHz低轨卫星地-空链路的大气信道衰减模型；Nessel初步估算71-86GHz的气体和云衰减影响。

现有技术主要采用遍历的方法，如将距离已知干扰地球站E_I一定距离范围内(球面圆内)进行栅格划分，或在干扰地球站E_I一定距离上(球面圆上)采样取点的方式，遍历计算每一栅格处或采样点对应的干扰值，基于数据拟合等方法找到满足要求的地球站位置，计算复杂度为O(n²)或O(n)，存在计算过程复杂的问题。

发明内容

本发明的目的在于克服现有技术缺陷，本发明提出了一种基于干扰分布的受扰系统地球站位置计算方法，用于在干扰系统地球站E_I外设定的区域范围内，确定满足干扰限值要求的最近受扰系统地球站E_x的位置，所述方法包括：

步骤1)获取干扰系统卫星S_I、干扰系统地球站E_I和受扰系统卫星S的位置数据，以及干扰系统地球站E_I外设定的区域范围[0,length]，令H＝length；

步骤2)计算在干扰系统地球站E_I外设定距离为H时，干扰最小的受扰系统地球站E_x位置、最小干扰值I/N，以及受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角α_H；

步骤3)判断步骤2)得到的最小干扰值I/N是否满足限值要求，判断为是，令α_max＝α_H，采用二分法计算I/N满足限值要求的夹角最小值α_min，令α＝α_min，转至步骤4)；判断为否，无满足的受扰系统地球站位置，转至步骤5)；

步骤4)根据α计算受扰系统地球站E_x的位置坐标；

步骤5)结束。

作为上述方法的一种改进，所述步骤1)具体包括：

获取干扰系统卫星S_I的位置坐标为(x_SI,y_SI,0)，干扰系统地球站E_I的坐标为(x_EI,y_EI,z_EI)，经纬度坐标为受扰系统卫星S的位置坐标为(x_S,y_S,0)，受扰系统卫星S的接收天线离轴角为θ，地球半径为R，S_I和S均为GSO卫星，轨道高度均为h，干扰系统地球站E_I外设定的距离为H，干扰系统地球站E_I外设定的区域范围为指定的球面区间[0,length]，令H＝length。

作为上述方法的一种改进，所述步骤2)具体包括：

步骤s2-1)计算干扰系统地球站E_I与E_I外设定距离H处到地心的夹角α_H；

步骤s2-2)基于几何关系得到干扰分布特征，根据干扰分布特征计算夹角为α_H时对应的最小干扰值I/N；

步骤s2-3)根据α_H计算受扰系统地球站E_x的位置坐标。

作为上述方法的一种改进，所述步骤s2-1)的干扰分布特征为：对于一个给定的夹角α，I/N在η＝90°处有最小值，η为受扰系统地球站E_x与受扰系统卫星s以及干扰系统地球站E_I构成的夹角。

作为上述方法的一种改进，所述步骤s2-2)最小干扰值I/N满足下式：

其中，I表示干扰信号功率；N表示受扰系统接收机的等效噪声功率；α为受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角，β为受扰系统地球站E_x、受扰系统卫星S和地心的夹角；p_I表示干扰系统地球站E_I的发射功率，g_I(γ)表示干扰系统地球站天线的发射增益；g(θ)表示受扰系统卫星天线S的接收增益；t_S表示受扰系统卫星S接收端的等效噪声温度；b表示受扰系统上行链路的通信带宽；k表示波尔兹曼常数；l_W表示干扰系统链路的传输损耗。

作为上述方法的一种改进，所述步骤s2-3)具体包括：

根据下式得到受扰系统地球站E_x的位置坐标(x,y,z)：

其中，p、q、m、n均为中间变量，分别满足下式：

(d₄)²＝2R²-2R²cosα_H

(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

作为上述方法的一种改进，所述步骤3)的限值要求为I/N≤-12.2dB。

作为上述方法的一种改进，所述步骤4)具体包括：

根据下式得到受扰系统地球站E_x的位置坐标(x,y,z)：

其中，p、q、m、n均为中间变量，分别满足下式：

(d₄)²＝2R²-2R²cosα

(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

其中，d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

与现有技术相比，本发明的优势在于：

1、本发明提出了的基于干扰分布的受扰系统地球站位置计算方法，分别针对在干扰系统地球站E_I外设定的区域范围内，确定满足干扰限值要求的最近受扰地球站E_x的位置，以及在干扰系统地球站E_I外设定的距离上，确定干扰最小的受扰地球站E_x位置两种情况；

2、通过本发明所提方法，可将上述两者情况的计算复杂度分别降至O(n)和O(1)，有效降低计算复杂度。

附图说明

图1是GSO卫星系统间上行链路干扰分析场景；

图2是本发明实施例1的区域内满足干扰限值的最近受扰系统地球站计算流程；

图3是基于干扰分布的最小干扰地球站位置计算流程；

图4是基于二分法计算满足干扰限制要求的最小α；

图5(a)是受扰系统地球站位置二元函数模型中参数α角度值分布；

图5(b)是受扰系统地球站位置二元函数模型中参数β角度值分布；

图5(c)是受扰系统地球站位置二元函数模型中参数θ角度值分布；

图5(d)是受扰系统地球站位置二元函数模型中参数η角度值分布。

具体实施方式

1.GSO系统间上行链路干扰分析模型构建

构建两GSO(对地静止轨道)卫星系统上行链路干扰分析场景如图1所示，其中，卫星S与地球站E_x作为受扰通信系统(系统1)，卫星S_I与地球站E_I作为潜在干扰通信系统(系统2)。

采用干扰噪声比I/N表示受扰卫星S接收到来自地球站E_I的潜在干扰如公式(1)所示。根据ITU-R S.1432-1建议书，采用I/N作为评估干扰的门限值为-12.2dB。

其中：I表示干扰信号功率；N表示受扰系统接收机的等效噪声功率；p_I表示干扰地球站的发射功率；g_I(γ)表示干扰地球站天线的发射增益；g(θ)表示受扰卫星天线的接收增益；t_S表示受扰卫星接收端的等效噪声温度；b表示受扰系统上行链路的通信带宽；k表示波尔兹曼常数；l_W表示干扰链路的传输损耗。

设受扰系统地球站E_x的经纬度坐标为如图1所示为GSO卫星系统间上行链路干扰分析场景。若固定干扰系统地球站E_I、干扰系统卫星S_I、受扰系统卫星S三点，变换受扰系统地球站的位置E_x，可改变卫星S的接收天线离轴角θ，从而影响干扰指标I/N。设已知E_I坐标为/>令Δλ＝|λ_x-λ|，/>则图1中的站-站地心角α(即，受扰地球站E_x与干扰地球站E_I间的地心角)如公式(2)所示。

同理可计算站-星地心角β(即，受扰地球站E_x与受扰卫星S间的地心角)。由此，将E_x的位置坐标转换为(α,β)二元组，且α、β∈[0,π]。

根据余弦定理可知，图1中卫星S的接收天线离轴角θ计算如公式(3)所示。

式中：d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

基于受扰地球站E_x的相对位置二元组(α,β)对公式(3)进行变换，得如公式(4)所示。

式中：R为地球半径，R_g为GSO卫星轨道半径。

结合公式(1)(4)，可得受扰系统地球站位置二元组(α,β)与干扰指标I/N的数学模型如公式(5)所示。

2干扰分布特征分析

晴空条件下，干扰链路的传输损耗l_W近似为自由空间损耗，若考虑天气因素，则传输损耗l_W为自由空间损耗叠加干扰链路的降雨、云雾等衰减损耗。因此，在仅改变受扰系统地球站E_x的场景下，干扰链路的传输损耗l_W为常量，则公式(5)中部分可视为常系数。由此，在受扰系统地球站位置二元函数模型I/N(θ(α,β))中，随着(α,β)的变化产生变化的主要是θ(α,β)部分。

1)分析θ关于α的变化特性，基于公式(4)进行变换，如公式(6)所示，根据复合函数求偏导法则可知，即在β一定的条件下，随着α的增加，θ单调递增。

2)类似地，分析θ关于β的变化特性，将公式(4)进行变换，如公式(7)所示。

令η＝∠SE_xE_I，由公式(7)可知，当即η＜90°时，/>表明在α一定的条件下，随着β的增加，θ单调递减；当/>即η＞90°时，/>表明在α一定的条件下，随着β的增加，θ单调递增。

综上，如图1所示，若β固定不变，则受扰地球站E_x的位置轨迹将构成以O为顶点，OS为轴，β为半顶角的球顶锥体的球面圆C₂，θ在圆C₂上随着α角的增大而增大。同理，若α固定不变，则受扰地球站E_x的位置轨迹将构成以O为顶点，OE_I为轴，β为半顶角的球顶锥体的球面圆C₁，θ在圆C₁上一部分(η＞90°时)随着β角的增大而增大，一部分(η＜90°时)随着β角的增大而减小。而I/N模型中g(θ)为受扰卫星天线的接收增益函数，与卫星采用的天线方向图有关，若g(θ)随θ的增大而减小，则根据复合函数单调性可知，I/N关于α、β变化的趋势与θ相反。

因此，当α固定，即E_x与E_I的距离L固定时，I/N在η＝90°处有最小值。

下面结合附图和实施例对本发明的技术方案进行详细的说明。

实施例1

如图2所示，本发明的实施例1提供了一种基于干扰分布的受扰系统地球站位置计算方法，用于在干扰地球站E_I外一定区域内，确定满足干扰限值的最近受扰地球站位置E_x，所述方法包括：

步骤1)获取干扰系统卫星S_I、干扰系统地球站E_I和受扰系统卫星S的位置数据，以及干扰地球站E_I外设定的区域范围[0,length]，令H＝length；

具体包括：

获取干扰系统卫星S_I的位置坐标为(x_SI,y_SI,0)，干扰系统地球站E_I的坐标为(x_EI,y_EI,z_EI)，经纬度坐标为受扰系统卫星S的位置坐标为(x_S,y_S,0)，受扰系统卫星S的接收天线离轴角为θ，地球半径为R，S_I和S均为GSO卫星，轨道高度均为h，干扰系统地球站E_I外设定的距离为H，干扰系统地球站E_I外设定的区域范围为指定的球面区间[0,length]，令H＝length。

步骤2)计算在干扰系统地球站E_I外设定距离为H时，干扰最小的受扰地球站E_x位置、最小干扰值I/N，以及受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角α_H；

如图3所示，具体包括：

步骤s2-1)计算干扰系统地球站E_I与E_I外设定距离H处到地心的夹角α_H；

步骤s2-2)基于几何关系得到干扰分布特征，根据干扰分布特征计算夹角为α_H时对应的最小干扰值I/N；

步骤s2-3)根据α_H计算受扰系统地球站E_x的位置坐标。

具体包括：

根据下式得到受扰系统地球站E_x的位置坐标(x,y,z)：

其中，p、q、m、n均为中间变量，分别满足下式：

(d₄)²＝2R²-2R²cosα_H

(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

其中，d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

步骤3)判断步骤2)得到的最小干扰值I/N是否满足限值要求，判断为是，令α_max＝α_H，采用二分法计算I/N满足限值要求的夹角最小值α_min，令α＝α_min，转至步骤4)；判断为否，无满足的受扰地球站位置，转至步骤5)；二分法具体过程见图4；

步骤4)根据α计算受扰系统地球站E_x的位置坐标；

步骤5)结束。

对于步骤4)参见图3，基于二分法计算满足干扰限制要求的最小α_min。

如图4为基于干扰分布的最小干扰地球站位置计算流程。

其中，

(d₄)²＝2R²-2R²cosα

(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

以下是上述解对应的方程组的求解过程：

设受扰系统地球站E_x的坐标为(x,y,z)，已知干扰系统卫星S_I的位置坐标为(x_SI,y_SI,0)、干扰系统地球站E_I的坐标为(x_EI,y_EI,z_EI)、受扰系统卫星S的位置坐标为(x_S,y_S,0)，地球半径为R，GSO卫星S_I、S的轨道高度为h。

根据干扰分布特征，地球站E_x满足如下方程组

其中，(d₄)²＝2R²-2R²cosα，(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²。

令③-①，得

令③-②，得

令

则，y＝px+q，z＝mx+n，将其带入③得：

(1+p²+m²)x²+(2pq+2mn)x+q²+n²-R²＝0

令：a＝1+p²+m²，b＝2pq+2mn，c＝q²+n²-R²，则：

求得方程组的解：

实施例2

如图3所示，本发明的实施例2提供了一种基于干扰分布的受扰系统地球站位置计算方法，该方法也是实施例1的一个步骤；在本实施例中用于基于干扰地球站E_I外设定的距离H上，确定干扰最小的受扰地球站E_x位置，所述方法包括：

步骤1)获取干扰系统卫星S_I、干扰系统地球站E_I和受扰系统卫星S的位置数据，以及干扰地球站E_I外设定的距离H；

步骤2计算干扰系统地球站E_I与E_I外设定距离H处到地心的夹角α_H；

步骤3)基于几何关系得到干扰分布特征，根据干扰分布特征计算夹角为α_H时对应的最小干扰值I/N；

步骤4)根据α_H计算受扰系统地球站E_x的位置坐标。

根据下式得到受扰系统地球站E_x的位置坐标(x,y,z)：

其中，p、q、m、n均为中间变量，分别满足下式：

(d₄)²＝2R²-2R²cosα_H

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

其中，d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

实验验证与分析

基于图1中建立的GSO卫星系统间上行链路干扰分析场景，以ASIASAT-105.5P(轨道经度105.5°E)作为受扰系统的GSO卫星S，选取干扰系统的GSO卫星S_I与受扰卫星S轨位间隔0°和2°(S_I轨道经度分别为105.5°E和107.5°E)，干扰系统地球站E_I以密云(40.3N,116.8E)站为例，分析受扰系统地球站位置变化对干扰指标I/N的影响。其中，受扰系统卫星参数、干扰系统卫星参数如表1所示，卫星天线方向图参考建议书ITU-R S.672-4，地球站天线方向图参考建议书ITU-R S.580-6。

表1干扰卫星上行链路的波束及空口参数

基于构建的干扰仿真场景，计算受扰地球站E_x在60°E～160°E经度范围，20°S～70°N纬度范围内满足I/N≤-12.2dB的最近距离位置。

干扰分布特征验证：

基于传统栅格化方法，一方面作为本发明方法的比对验证，另一方面验证本发明的干扰分布特征的正确性。通过以1°间隔为单位进行栅格化计算α、β、θ、η的分布如图5(a)、图5(b)、图5(c)和图5(d)所示。

从图5中可见，随着受扰地球站E_x的位置变换，α角以干扰地球站E_I为中心延展增大，β角以受扰系统卫星S为中心延展增大。由第2节中对I/N的函数分析得知，θ关于β的变化趋势，与η值相关，计算η的角度分布，并将η＝90°等值线标出如图5(d)所示。选择图5(a)中的某一α等值线，以α＝40°为例，在该等值线上取A、B、C三点，分析β从A到B、从B到C变化对θ的影响，A、B、C三点相关数据如表2所示。结合图5及表2数据可知，AB段位于η＞90°的区域内，从A到B，β递增，且θ递增；BC段位于η＜90°区域，从B到C，β递增，但θ递减。同理，结合图5(a)、5(b)、5(c)容易看出，在β一定的条件下，随着α的增加，θ单调递增。以上仿真分析结果与函数特性分析一致。

表2 A、B、C三点的坐标、角度及干扰指标数据

实验结果对比分析：

实验(1)：区域上满足干扰限值的最近距离点

基于构建的干扰仿真场景，计算受扰地球站E_x在60°E～160°E经度范围，20°S～70°N纬度范围内满足I/N≤-12.2dB的最近距离位置。

基于传统栅格化方法，栅格化方法以1°间隔进行栅格，通过数据拟合的方式，计算满足I/N≤-12.2dB(精度要求ε≤0.05dB)的最近地球站位置坐标为(73.80E，22.99N)，I/N＝-12.21dB，α＝40.00°。相邻栅格数据如表3所示。

表3传统栅格化方法的部分栅格数据结果

基于本发明所提方法计算得出两个坐标值，分别为(1633.4,5636.7,2479.8)、(5053.1,3782.4,865.8)，对应经纬度为(73.84E,23.05N)、(143.18E,7.86N)，I/N＝-12.21dB，α＝39.99°。

综上，两种方法计算结果相近，差别在于栅格化方法在数据拟合时会存在误差。栅格化方法计算次数为100*90＝9000次，本发明所提方法，计算次数为9次。

实验(2)：指定距离上的最小干扰值及其位置计算

以α＝25°的距离为例，计算该距离上的最小干扰值。

基于传统取点方式，取236个点，通过插值方法，计算最小干扰为-7.27dB，对应的坐标位置是(23.75E,139.07N)。

基于本发明所提方法计算最小干扰值I/N＝-7.21dB，两处坐标值分别是(408.3,5157.0,3718.7)、(-4412.3,3820.1,2555.1)，经纬度为(35.89E,85.47N)、(23.79E,139.11N)。

综上，两种方法计算结果相近，传统取点方法计算次数为236次，本发明所提方法，计算次数为1次。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (2)

1.一种基于干扰分布的受扰系统地球站位置计算方法，用于在干扰系统地球站E_I外设定的区域范围内，确定满足干扰限值要求的最近受扰系统地球站E_x的位置，所述方法包括：

步骤1)获取干扰系统卫星S_I、干扰系统地球站E_I和受扰系统卫星S的位置数据，以及干扰系统地球站E_I外设定的区域范围[0,length]，令H＝length；

步骤2)计算在干扰系统地球站E_I外设定距离为H时，干扰最小的受扰系统地球站E_x位置、最小干扰值I/N，以及受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角α_H；

步骤3)判断步骤2)得到的最小干扰值I/N是否满足限值要求，判断为是，令α_max＝α_H，采用二分法计算I/N满足限值要求的夹角最小值α_min，令α＝α_min，转至步骤4)；判断为否，无满足的受扰系统地球站位置，转至步骤5)；

步骤4)根据α计算E_I到E_x的距离d₄，采用步骤s2-3)计算受扰系统地球站E_x的位置坐标；

步骤5)结束；

所述步骤1)具体包括：

获取干扰系统卫星S_I的位置坐标为(x_SI,y_SI,0)，干扰系统地球站E_I的坐标为(x_EI,y_EI,z_EI)，经纬度坐标为受扰系统卫星S的位置坐标为(x_S,y_S,0)，受扰系统卫星S的接收天线离轴角为θ，地球半径为R，S_I和S均为GSO卫星，轨道高度均为h，干扰系统地球站E_I外设定的距离为H，干扰系统地球站E_I外设定的区域范围为指定的球面区间[0,length]，令H＝length；

所述步骤2)具体包括：

步骤s2-1)计算干扰系统地球站E_I与E_I外设定距离H处到地心的夹角α_H；

步骤s2-2)基于几何关系得到干扰分布特征，根据干扰分布特征计算夹角为α_H时对应的最小干扰值I/N；

步骤s2-3)根据α_H计算受扰系统地球站E_x的位置坐标；

所述步骤s2-2)的干扰分布特征为：对于一个给定的夹角α，I/N在η＝90°处有最小值，η为受扰系统地球站E_x与受扰系统卫星s以及干扰系统地球站E_I构成的夹角；

所述步骤s2-2)最小干扰值I/N满足下式：

其中，I/N(α,β)表示与受扰系统地球站位置(α,β)二元组相关的最小干扰值，I表示干扰信号功率；N表示受扰系统接收机的等效噪声功率；α为受扰系统地球站E_x、干扰系统地球站E_I与地心的夹角，β为受扰系统地球站E_x、受扰系统卫星S和地心的夹角；p_I表示干扰系统地球站E_I的发射功率，g_I(γ)表示干扰系统地球站天线的发射增益；g(θ(α,β))表示与(α,β)二元组相关的受扰系统卫星S天线的接收增益；t_S表示受扰系统卫星S接收端的等效噪声温度；b表示受扰系统上行链路的通信带宽；k表示波尔兹曼常数；l_W表示干扰系统链路的传输损耗；

所述步骤s2-3)具体包括：

根据下式得到受扰系统地球站E_x的位置坐标(x,y,z)：

其中，p、q、m、n均为中间变量，分别满足下式：

(d₄)²＝2R²-2R²cosα_H

(d₃)²＝(d₂)²-(d₄)²

(d₂)²＝(x_S-x_EI)²+(y_S-y_EI)²+(z_S-z_EI)²

其中，d₂为E_I到S的距离，d₃为E_x到S的距离，d₄为E_I到E_x的距离。

2.根据权利要求1所述的基于干扰分布的受扰系统地球站位置计算方法，其特征在于，所述步骤3)的限值要求为I/N≤-12.2dB。