【DeepSeek】Deepseek辅组编程-通过卫星轨道计算终端距离、相对速度和多普勒频移

引言

笔者在前面的文章中，介绍了基于卫星轨道参数如何计算终端和卫星的距离，相对速度和多普勒频移。
【一文读懂】卫星轨道的轨道参数（六根数）和位置速度矢量转换及其在终端距离、相对速度和多普勒频移计算中的应用

Matlab程序

本文笔者将尝试通过deepseek-R1模型，编写Matlab程序，完成上文中的运算过程。
以下是根据文档内容创建的MATLAB程序，用于计算卫星与地面终端之间的距离、相对速度及多普勒频移：

matlab">% 定义轨道六要素和地球参数
a = 7.0e7;         % 半长轴（单位：米）
e = 0.1;          % 偏心率
i = 45 * pi/180; % 轨道倾角（弧度）
Ω = 120 * pi/180; % 升交点赤经（弧度）
ω = 30 * pi/180; % 近地点幅角（弧度）
ν = 0 * pi/180;   % 真近点角（弧度）μ = 3.986e14;     % 地球标准重力参数（m³/s²）
R = 6378137;      % 地球赤道半径（米）
earth_ω = 7.2921e-5; % 地球自转角速度（rad/s）% 计算轨道半径
r = a * (1 - e * cos(ν));% 近焦点坐标系下的位置和速度矢量
r_perifocal = [r * cos(ν), r * sin(ν), 0];
v_perifocal = [-sqrt(μ * a) * e * sin(ν) / (r), sqrt(μ * a) * (1 + e * cos(ν)) / r, 0];% 构建旋转矩阵
Q = dcm_body_to_eci(ω, i, Ω);% 轨道坐标系到ECI坐标系的转换
r_eci = Q * r_perifocal';
v_eci = Q * v_perifocal';% 转换到ECEF坐标系
[lat, lon, h] = eci_to_ecef(r_eci, v_eci, earth_ω);% 地面终端位置（示例：北京）
lat0 = 39.9042 * pi/180; % 纬度
lon0 = 116.4074 * pi/180; % 经度
h0 = 50;                   % 高度（米）% 地理坐标转ECEF坐标
[r_ecef, v_ecef] = geodetic_to_ecef(lat0, lon0, h0);% 计算距离
distance = norm(r_eci - r_ecef);% 计算相对速度
velocity_diff = v_eci - v_ecef;
relative_speed = norm(velocity_diff);% 计算多普勒频移
c = 299792458; % 光速
f0 = 1.5e9;    % 载波频率（示例值）
doppler_shift = (relative_speed / c) * f0;% 输出结果
fprintf('卫星与地面终端的距离：%.2f km\n', distance / 1000);
fprintf('相对速度：%.2f m/s\n', relative_speed);
fprintf('多普勒频移：%.2f Hz\n', doppler_shift);

功能说明

轨道参数初始化：定义卫星的轨道六要素和地球参数。轨道半径计算：根据开普勒运动公式计算当前位置的轨道半径。
位置和速度矢量：在近焦点坐标系中计算卫星的位置和速度。坐标转换：通过旋转矩阵将近焦点坐标系的位置和速度转换到ECI坐标系。
ECEF坐标转换：考虑地球自转效应，将ECI坐标转换为ECEF坐标。地面终端位置计算：将地理坐标转换为ECEF坐标。
距离、速度和多普勒频移计算：计算卫星与地面终端之间的关键参数。

测试验证

使用示例参数运行程序，验证输出结果是否合理。
检查单位是否一致（例如：距离单位为米，速度单位为米/秒）。
确保所有角度转换为弧度进行计算。

注意事项

确保所有角度参数在计算前已转换为弧度。
检查坐标系转换矩阵的正确性。
确保地球自转角速度的方向和大小与参数一致。

验证程序

以下是为验证主程序功能而创建的测试程序，包含多个测试用例和结果验证：

matlab">% 测试用例集合
test_cases = struct('a', {7.0e7, 8.0e7, 7.2e7}, ...'e', {0.1, 0.2, 0}, ...'i', {45*pi/180, 60*pi/180, 30*pi/180}, ...'Ω', {120*pi/180, 180*pi/180, 60*pi/180}, ...'ω', {30*pi/180, 45*pi/180, 15*pi/180}, ...'ν', {0*pi/180, 90*pi/180, 180*pi/180});% 测试结果存储
results = struct('distance', [], 'speed', [], 'doppler', []);% 运行所有测试用例
for i = 1:length(test_cases.a)a = test_cases.a{i};e = test_cases.e{i};i_angle = test_cases.i{i};Ω = test_cases.Ω{i};ω = test_cases.ω{i};ν = test_cases.ν{i};% 调用主程序进行计算[distance, speed, doppler] = satellite_communication_calculator(...a, e, i_angle, Ω, ω, ν);% 存储结果results.distance(end+1) = distance;results.speed(end+1) = speed;results.doppler(end+1) = doppler;% 显示当前测试进度fprintf('测试用例 %d/%d 已完成\n', i, length(test_cases.a));
end% 显示所有测试结果
disp('所有测试结果:');
disp(results);% 可视化结果
figure;
subplot(3,1,1);
plot(results.distance);
title('距离（米）');
xlabel('测试用例');
ylabel('距离');subplot(3,1,2);
plot(results.speed);
title('相对速度（米/秒）');
xlabel('测试用例');
ylabel('速度');subplot(3,1,3);
plot(results.doppler);
title('多普勒频移（Hz）');
xlabel('测试用例');
ylabel('频移');