基于MATLAB的ATK二次开发指南：卫星对地目标可见性分析优化
一、ATK与卫星可见性任务
ATK（航天任务设计工具箱）是我校自主研发的航天任务仿真平台，对标国际主流工具STK，支持轨道设计、可见性分析、多域任务建模等功能。在卫星对地观测任务中，可见性分析是核心环节，直接影响卫星与地面站的通信效率、数据回传窗口规划等关键指标。通过MATLAB二次开发，可实现自动化参数配置、批量化仿真计算、数据智能分析，大幅提升任务设计效率。本篇文章基于如下任务，对基于MATLAB的ATK二次开发进行详解：
任务内容：针对一颗24小时内会经过文昌站(19.62, 110.75)上空的圆轨道卫星，卫星高度为500km，现有三个地面测站目标：Seoul(37.59850， 126.9780)；Taipei(25.04780， 121.5320)；Ukrainsk(48.10， 37.38330); 请设计卫星轨道，使得该卫星对三个地面站在24小时内观测时间最长。指定起始时间为2025-03-21 00:00:00。
二、任务过程分析
该任务本质为一个目标优化任务，优化目标为卫星对三个地面站的可见时间最长，同时24小时内对三个地面站均可见。任务对目标的约束为：卫星高度500km；圆轨道；24小时内过文昌上空。
从0时刻开始计，假定时间后，卫星经过文昌站正上方；由于圆轨道卫星高度已定，所以此时的卫星速度大小可以唯一确定。如果在此时定义好速度的方向，则可以获得此刻卫星的位置速度，进而可以计算出卫星初始的轨道根数。该问题为转化为“二自由度优化问题+由初始轨道根数预报卫星可见性”的分步问题。
由于约束的离散型，采用粒子群算法particleswarm，定义自变量为经过文昌上空的时间与卫星的速度夹角。定义输出变量为三个站点可见时间总和；若有站点不可见则总可见时间为0（即当前轨道无效）。
至此，问题的输入输出已确认完毕。针对轨道预报与可见时间计算，采用matlab对ATK用connect模式进行连接，输入初始轨道根数进行计算。
下面介绍MATLAB对ATK的Connect操作代码。
三、MATLAB&ATK连接
-
路径配置
在MATLAB中添加ATK的二次开发接口路径，确保调用接口函数时无冲突：1
addpath('D:\ATK-3.2.0\IntegratingWithATK\connect\win\Matlab\Win_2015b');
-
进程管理
启动ATK前需关闭残留进程，避免端口占用。由于ATK在多次开关下存在僵尸后台进程，故需要用Windows的taskkill与wmic命令，关闭其父进程与后台僵尸子进程。进而用命令行调用ATK启动：1
2
3
4[~,~] = system('taskkill -f -im ATK.exe');
[~, ~] = system(strcat('wmic process where "name', "='ATK.exe'", '" delete'));
system('call start /min "n" "D:\ATK-3.2.0\ATK.exe" &');
pause(5) % 等待ATK初始化完成 -
建立连接
通过IP和端口实现MATLAB与ATK的通信：1
2conID = atkOpen('127.0.0.1', 6655);
atkConnect(conID, 'New', '/ Scenario InclinationChange1'); -
计算起止时间与创建新场景
首先用MATLAB的
datetime函数定义起止时间，进而用datestr函数转化为需要的字符串格式1
2
3
4starttime = input_time();
stoptime = starttime + seconds(86400);
starttime_str = datestr(starttime, 'dd mmm yyyy HH:MM:SS');
stoptime_str = datestr(stoptime, 'dd mmm yyyy HH:MM:SS');创建名为
InclinationChange1的新场景，并定义场景时间：1
2
3
4
5
6% 创建新场景
atkConnect(conID, 'New', '/ Scenario InclinationChange1');
% 设置场景起止时间
atk_starttime_str = ['* "', starttime_str, '" "', stoptime_str, '"'];
atkConnect(conID, 'SetAnalysisTimePeriod', atk_starttime_str);
atkConnect(conID, 'Animate', '* Reset');
四、卫星轨道与地面站参数设置
-
轨道参数
• 传入卫星的初始轨道六根数Struct结构体，其结构格式为OrbEmt.a; OrbEmt.e; OrbEmt.i; OrbEmt.Omega; OrbEmt.omega; OrbEmt.f; 定义轨道六根数• 新建卫星，并使用
TwoBody模型初始化卫星轨道：1
2
3
4
5
6
7
8
9
10
11% 字符串化轨道六根数
atk_a = sprintf('%.10f ', OrbEmt.a * 1000);
atk_e = sprintf('%.10f ', OrbEmt.e);
atk_i = sprintf('%.10f ', OrbEmt.i/pi*180);
atk_Omega = sprintf('%.10f ', OrbEmt.Omega/pi*180);
atk_omega = sprintf('%.10f ', OrbEmt.omega/pi*180);
atk_f = sprintf('%.10f', OrbEmt.f/pi*180);
%% 新建卫星
atkConnect(conID, 'New', '/ Satellite Satellite1');
% 设置参数
atkConnect(conID, 'SetState', [' */Satellite/Satellite1 Classical TwoBody "', starttime_str,'" "', stoptime_str, '" 1 J2000 "', starttime_str, '" ', atk_a, atk_e, atk_i, atk_omega, atk_Omega, atk_f]); -
地面站部署
• 定义多个地面站，设置地理坐标：1
2
3
4
5
6
7
8
9
10%% 新建地面站
atkConnect(conID, 'New', '/ Facility Seoul');
atkConnect(conID, 'New', '/ Facility Taipei');
atkConnect(conID, 'New', '/ Facility Ukrainsk');
atkConnect(conID, 'New', '/ Facility WenChang');
% 设置地面站位置
atkConnect(conID, 'SetPosition', [' */Facility/Seoul Geodetic ', sprintf('%.10f %.10f %.10f', Seoul0(1), Seoul0(2), Seoul0(3))]);
atkConnect(conID, 'SetPosition', [' */Facility/Taipei Geodetic ', sprintf('%.10f %.10f %.10f', Taipei0(1), Taipei0(2), Taipei0(3))]);
atkConnect(conID, 'SetPosition', [' */Facility/Ukrainsk Geodetic ', sprintf('%.10f %.10f %.10f', Ukrainsk0(1), Ukrainsk0(2), Ukrainsk0(3))]);
atkConnect(conID, 'SetPosition', [' */Facility/WenChang Geodetic ', sprintf('%.10f %.10f %.10f', WenChang0(1), WenChang0(2), WenChang0(3))]);
五、传感器视场约束与可见性分析
-
传感器配置
• 添加卫星传感器并设置视场角（如20°半锥角）：1
2
3
4
5
6
7
8% 定义半锥角
half_theta = 20;
%% 添加敏感器
atkConnect(conID, 'New', ' / */Satellite/Satellite1/Sensor Sensor1');
% 勾选视场约束
atkConnect(conID, 'SetConstraint', ' */Satellite/Satellite1/Sensor/Sensor1 fieldofview On');
% 设置视角
atkConnect(conID, 'Define', [' */Satellite/Satellite1/Sensor/Sensor1 SimpleCone ', sprintf('%.2f', half_theta)]); -
可见性计算
• 调用Report_RM生成各地面站的访问时间窗口报告。注意，设置TimeStep较小，计算较为准确，但是同样耗时较高。若较大，则可能丢失数据。：1
2
3
4
5
6%% 计算可见性并导出报告
return1 = atkConnect(conID, 'Report_RM', [' */Satellite/Satellite1/Sensor/Sensor1 Style "Access" AccessObject */Facility/Seoul TimePeriod "', starttime_str, '" "', stoptime_str, '" TimeStep 1']);
return2 = atkConnect(conID, 'Report_RM', [' */Satellite/Satellite1/Sensor/Sensor1 Style "Access" AccessObject */Facility/Taipei TimePeriod "', starttime_str, '" "', stoptime_str, '" TimeStep 1']);
return3 = atkConnect(conID, 'Report_RM', [' */Satellite/Satellite1/Sensor/Sensor1 Style "Access" AccessObject */Facility/Ukrainsk TimePeriod "', starttime_str, '" "', stoptime_str, '" TimeStep 1']);
return4 = atkConnect(conID, 'Report_RM', [' */Satellite/Satellite1/Sensor/Sensor1 Style "Access" AccessObject */Facility/WenChang TimePeriod "', starttime_str, '" "', stoptime_str, '" TimeStep 1']); -
数据解析
• 由于Report_RM得到的报告并非完全格式化，这里要用正则提取需要的数据：自定义函数parseVisibilityReport解析返回的可见性数据，提取时间窗口和持续时间：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73accessData1 = parseVisibilityReport(return1);
all_time1 = readAccess(accessData1); % 统计总可见时长
function accessData = parseVisibilityReport(reportStr)
% 解析卫星可见性报告字符串
% 输入参数：
% reportStr - 包含可见性报告的字符串
% 输出参数：
% accessData - 结构体数组，包含以下字段：
% target: 目标名称
% startTime: 开始时间(datetime)
% endTime: 结束时间(datetime)
% duration: 持续时间(秒)
% 初始化结构体数组
accessData = struct('target', {}, 'startTime', {}, 'endTime', {}, 'duration', {});
% 正则表达式模式（关键解析逻辑）
% pattern_locate = [
% 'Satellite\d+-Sensor\d+对(\w+)的可见性报告[^\n]*'... % 匹配目标名称
% '[\s\S]*?'... % 非贪婪匹配任意字符
% '(\d+)\s+(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\.\d{3})'... % 匹配序号和开始时间
% '\s+(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\.\d{3})'... % 匹配结束时间
% '\s+(\d+\.\d{3})' % 匹配持续时间
% ];
pattern_locate = 'Satellite\d+-Sensor\d+对(\w+)的可见性报告[^\n]*';
pattern = [
'(\d+)\s+(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\.\d{3})'... % 匹配序号和开始时间
'\s+(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\.\d{3})'... % 匹配结束时间
'\s+(\d+\.\d{3})' % 匹配持续时间
];
% 执行正则匹配
[matches, tokens_locate] = regexp(reportStr, pattern_locate, 'match', 'tokens');
[matches, tokens] = regexp(reportStr, pattern, 'match', 'tokens');
% 遍历所有匹配项
indexs = [];
for i = 1:length(tokens)
target = tokens_locate{1, 1};
target = target(1);
index = str2num(tokens{1, i}{1});
if ismember(index, indexs)
continue
else
indexs = [indexs, index];
end
startTime = datetime(tokens{1, i}{2}, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSS');
endTime = datetime(tokens{1, i}{3}, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSS');
duration = str2double(tokens{1, i}{4});
% 添加到结构体数组
accessData(end+1) = struct(...
'target', target,...
'startTime', startTime,...
'endTime', endTime,...
'duration', duration);
end
% 数据有效性验证
if isempty(accessData)
accessData(1).duration = 0;
end
end
function all_time = readAccess(accessData)
N = length(accessData);
timelisti = zeros(N, 1);
for i = 1:N
timei = accessData(i).duration;
timelisti(i) = timei;
end
all_time = sum(timelisti);
end
六、定义优化函数与优化计算
-
定义由速度倾角与时间计算初始轨道根数函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48% 轨道根数
function [OrbEmt, check] = getOrbEmt(eopdata, const, v_theta, time_pass)
% v_theta: 轨道倾角，rad
% time_pass: 经过时间后相遇，单位s
% clear;clc;
% v_theta = 10 / 180 * pi;
% time_pass = 0 * 3600;
% 文昌
aim_station = [19.620833000000005, 110.751389000000003];
h = 500;
R_e = 6378.140;
mu = 3.986004415e5;
% 计算半长轴以及周期
a = h + R_e;
T = sqrt(a^3*4*pi^2/mu);
time = input_time();
v_norm = sqrt(mu/a);
% 经过时间
time1 = add_time(time, time_pass);
% MJDi = JDD(time1);
MJDi = juliandate(time1) - 2400000.5;
Me = ECI2ECF(MJDi);
[X, Y, Z] = Geodetic2Cartesian(aim_station(1)/180*pi, aim_station(2)/180*pi, h);
rei = [X, Y, Z];
ri = Me \ (rei');
[Y1] = ECEF2ECI(eopdata, const, MJDi, [rei, [1,1,1]]);
ri = Y1(1:3);
% 这里用z轴叉乘去定义一个x和y的单位矢量，去仿真
zi = [0, 0, 1];
rzy = cross(ri, zi);
rzy_ = rzy/norm(rzy);
rzx = cross(ri, rzy);
rzx_ = rzx/norm(rzx);
dvi = rzy_ * cos(v_theta) + rzx_ * sin(v_theta);
v = v_norm * dvi;
[OrbEmt] = PosVel2OrbEmt(ri', v);
% 反推0点的f
df = time_pass / T * 2*pi;
f1 = OrbEmt.f - df;
f1 = Angle0_2pi(f1);
OrbEmt.f = f1;
check = 1;
% fprintf('%.9f\n', [OrbEmt.a*1000, OrbEmt.e, OrbEmt.i/pi*180, OrbEmt.Omega/pi*180, OrbEmt.omega/pi*180, OrbEmt.f/pi*180]) -
定义被优化函数
1
2
3
4
5
6M = sum(station_time == 0);
if M > 0
values = 86400;
else
values = -sum(station_time);
end -
定义优化算法与相关细节数据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21clear;clc;
tic;
ObjectiveFunction = @(x) objective_function(x);
options = optimoptions('particleswarm', 'UseParallel', true, 'SwarmSize',100,'MaxStallIterations',100, 'HybridFcn',@fmincon);
nvars = 2;%变量个数
LB = [0,0];%定义域下限
UB = [2*pi,24];%定义域上限
[x,fval,exitflag,output] =particleswarm(ObjectiveFunction,nvars,LB,UB,options);%调用particleswarm函数
toc;
% 60-32s 0.882514534157354 0.135176571902176
% 20-85s 3.97034896303301 24
% 1-1405.7 0.875097451387032 0.292774115756876
%%
fprintf('Final Conclusion:\n')
fprintf('-----------------------------------\n')
fprintf('x = [')
fprintf('%.16f\t', x)
fprintf(']\n')
fprintf('Time = %.2f s\n', -fval)
save('out.mat')
七、优化结果展示与分析
-
最终优化后的初始轨道根数为：
Orbital Elements:
--------------------------------------
a : 6873.3468182288 km
e : 0.0003485573
i : 50.4178267682 deg
Omega : -86.7600402771 deg
omega : 205.6135266226 deg
f : 101.0355340694 deg
-
可见性时间为：
目标 覆盖次数 时间(s) Seoul 2 97.324 Taipei 1 50.330 Ukrainsk 2 89.738 文昌 1 53.313
七、结语
通过ATK的MATLAB二次开发，工程师可快速构建卫星任务闭环仿真系统，实现从参数配置到结果分析的全流程自动化。未来可进一步集成AI算法，优化轨道规划和资源调度策略，推动航天任务设计的智能化转型。
八、单次ATK操控代码附录
1 | clear all |
All articles on this blog are licensed under CC BY-NC-SA 4.0 unless otherwise stated.
