卫星或飞船的全程跟踪测控建模
（海南软件职业技术学院07级计算机网络技术专业 某生）
原题：（来自09年高教社杯数学建模大赛乙组赛题，有所改动。
）
卫星和飞船在国民经济和国防建设中有着重要的作用，对它们的发射和运行过程进行测控是航天系统的一个重要组成部分，理想的状况是对卫星和飞船（特别是载人飞船）进行全程跟踪测控。
而测控设备只能观测到所在点切平面以上的空域，且在与地平面夹角3度的范围内测控效果不
好，实际上每个测控站的测控范围只考虑与地平面夹角3度以上的空域。
请利用模型分析:如果一个卫星或飞船的运行轨道与地球赤道平面有固定的夹角a ，且在离地面高度为H 的球面S 上运行。
那么，至少应该建立多少个测控站才能对该卫星或飞船可能飞行的区域
全部覆盖以达到全程跟踪测控的目的？
问题分析：
由于测控设备只观测与地平面夹角3度以上的空域，显然测控设备观测的空间是个圆锥形空域。
圆锥的顶角是固定不变的，其底面大小随圆锥的高（卫星轨道高H ）变化。
又由于卫星或飞船的运行轨道与地球赤道平面有固定的夹角a ，而且由于地球自转时该卫星或飞船在飞行过程中相继两圈的经度有所变化，这就会导致卫星或飞船的运行轨迹在地球表面的投影将会覆盖北纬a →赤道→南纬a 之间的带状区域如图1。
特别地，当a=0的时候卫星或飞船的轨道和赤道共面。
当a=π/2的时候卫星的投影将会是整个地球表面。
显然夹角a 决定着卫星运行空域在地球表面的总投影区域。
所以最终可把需要求的变量用关于H 和a 的函数表示出来。
图1 卫星正常运行时的示意图 图2 地球和测控区域等在卫星轨道上的截面图
对模中的变量之间的函数关系进行整理。
1.测控点所观测空域在地球表面上的投影区域内的最长测地线L 的计算函数。
如图2，R 是地球半径，外圆是高为H 的卫星运的运行轨道，灰白色部分为测控站J 的盲区，虚线是过测控点的地面切线，L 是测控站J 所观测区域在地球表面上的投影的最长测地线。
卫星轨道与J 的观测临界面交于A 点，在△AJO 中：JO=R,AO=R+H,∠AJO=
π60
61
，解三角函数可得：6060cos arccos ππ-⎪⎭⎫
⎝⎛+=∠H
R R
AOJ 。
弧L 对应的圆心角为2∠AOJ ，于是有：R AOJ L ⋅∠=2
⎪⎪⎭⎫ ⎝
⎛-⎪⎭⎫ ⎝⎛+=⇒6060cos arccos 2ππH R R
R L (1) 2.当卫星或飞船轨道与赤道夹角固定时卫星轨道在地球表面上的投影带的宽度d 关于a 的函数。
d 是2a 对应的弧长，弧半径为地球半径R ，即得：
d=2aR (2)
3.测控模型中要用到多个测控点分布在同一纬线的情况，为了方便表达，定义那些穿过测控点的纬线称为测控纬线。
测控纬线长度用u 表示，控纬线的纬度用W 表示（为了便于数学计算，纬度也用弧度制表示）。
这里算出u 和W 的关系函数。
μ的直径等于2Rcos W ，进而可得出：
μ=2πRcosW (3)
模型建立及求解
1.当测控纬线条数m 已知时同一测控纬线上连接相邻两测控点的测地线长度l 的计算。
如图3，d=AB ，所有圆直径为L ，已知测控纬线数为m ，设AC=x ，DB=y ，从而可得方程关于x,y 的方程组：
L
x
y d
x y m =-=+-2)1(
⇒ 1
12++=
++-=
m L d y m L mL d x (4) 进一步计算l ，l 的计算可以近似地当做平面运
算，即： 图3
442
22x L l -= 2
2
12⎪⎭
⎫ ⎝⎛++--=⇒m L m L d L l (5)
2.测控纬线条数m 、测控纬线纬度w （弧度制单位）、相邻两测控点的测地线长度l 、和该测控
图3纬线上的测控点个数n 的建模。
因为，当测控纬线条数m 为奇数时，由于南北半球对称，赤道也将是一条测控纬线；当测控纬线为偶数时，所有测控纬线和测控点都在赤道两边，与m 为奇数时模型不同，所以要把m 取值的奇数偶数两种情况分开建模。
A.当m 为奇数的时候如图3，虚线圆是单个测控点所测空域的投影区域。
J 1J 2所在的测控纬线到赤道的测地线距离为等式(4)中的y ，结合式(4)得赤道一侧第i 条测控纬线的纬度：
()()
1++=
m R L d i W i ，再由Wi 结合式(2)、式(3)可进一步得出赤道一侧第i 测控纬线长度μi 、和该测控纬线上的测控点个数ni 的值，整理为表2。
图4
表2：测控纬线条数为奇数时N 的求解参照表
i W i μi
n i
0(赤道) —
2πR
()
()[]()
2
2
2112L mL d m L m R +--++π
i
()()
1++m R L d i ()()
1cos
2++m R L d i R π
()
()[]()()()
1cos
2112
2+++--++m R L d i L mL d m L m R π
()
∑-=+=2
/11
02
m i i
n
n N
(6)（n i 的取值规则：当n i 有小数部分时进1取整数。
B.当m 为偶数时如图4。
AE 等于式(4)中的y ，J 1J 2所在的测控纬线到赤道的测地线距离为等于AE/2=y/2，由于赤
道北第二道测控纬线到第一道测控纬线的测地线距离为y ，结合式(4)得可得：
()()()
1212++-=
m R L d i W i ，由W i 结合式(2)、式(3)可进一步得出μi 、n i 的值，整理为表3。
图4
表3：测控纬线条数为偶数时N 的求解参照表
i W i
μi
n i
i ()()()
1212++-m R L d i ()()()
1212cos
2++-⋅m R L d i R π ()
()[]()()()()
1212cos
21122
2++-+--++m R L d i L m L d m L m R π
∑==2/1
2m i i n N (7) （n i 的取值规则：当n i 有小数部分时进1取整数。
）
编程计算最小的N 的值。
1.Matlab 程序设计：
程序输入部分：输入地球半径。
卫星高H ，赤道与卫星轨道夹角a 等。
计算部分：
如果夹角a=0时卫星在赤道上。
N=2πR/L 。
否则用关于m 的for 循环语句进入计算。
for 语句中给m 赋初值，因为测控纬线条数必须大于投影带宽d 和单个测控点测控区域直径L 的比值才能完全覆盖卫星出现的空域。
d/L 苦有小数部分就进1取整数，把所得值赋给m 做初值。
在for 语句主体部分结合前面所建模型，分别计算对应的总测控点个数N 和m+1（即增加一条测控纬线）时需求的总测控点个数N1，当N<N1时N 就是所求的最少的测控点。
这时break 跳出for 循环。
计算单个测控点空域的投影面A 和卫星轨道投影U 只需要依所推出函数试计算。
输出部分：输出A 、U 、N 。
2.程序代码见附码。
附码：程序环境(Matlab6p5)
%输入部分
R=6371;N=0;N1=0;
%算出L 和d
H=input('请输入卫星的平均高度(单位：km)H='); L=2*R*(acos(R/(R+H)*cos(pi/60))-pi/60);
a=input('请输入卫星轨道与赤道夹角(用弧度制)a='); d=2*a*R;
%算出N
if a==0
N=2*pi*R/L;
else for m=ceil(d/L):50 %给m赋初值
if mod(m,2)==1
N=2*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2);
if m>1
for i=1:(m-1)/2
N=ceil(N)+4*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2)*cos(i*(d+L)/(R*m+R));
end
end
for i=1:(m+1)/2
N1=ceil(N1)+4*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2)*cos((2*i-1)*(d+L)/(R*m+R));
end
else for i=1:m/2
N=ceil(N)+4*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2)*cos((2*i-1)*(d+L)/(R*m+R));
end
N1=2*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2);
for i=1:(m+1)/2
N1=ceil(N1)+4*pi*R/sqrt(L^2-((2*d-m*L+L)/(m+1))^2)*cos(i*(d+L)/(R*m+R));
end
end
if N<=N1 break;
else N=0;N1=0;
end
end
end
str=['所需卫星测控站个数N=',num2str(ceil(N))];disp(str);。