【卫星】用于实时预测和跟踪地球静止轨道卫星的轨道附Matlab代码

最新推荐文章于 2025-10-04 09:12:06 发布

原创最新推荐文章于 2025-10-04 09:12:06 发布 · 574 阅读

23 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #数据结构 #开发语言

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

🍊个人信条：格物致知,完整Matlab代码及仿真咨询内容私信。

🔥 内容介绍

地球静止轨道卫星：轨道预测与跟踪的重要性

地球静止轨道位于赤道平面上空约 36000 公里处，运行于此轨道的卫星绕地球公转周期与地球自转周期严格同步，从地面上看，卫星就像 “悬停” 在天空某一固定位置。通信领域中，全球各地的电视信号、互联网数据传输，很大程度依赖这些卫星搭建起的 “太空桥梁”，实现跨洋、跨洲的信息快速传递。气象监测时，卫星能持续捕捉地球云图、温度、湿度等气象数据，为天气预报、灾害预警提供一手资料，帮助人类提前防范台风、暴雨等极端天气。军事侦察中，其凭借居高临下的视野优势，对地面目标进行全天候监控，守护国家安全。

但地球静止轨道并非 “平静地带”，卫星会受到地球非球形引力、日月引力摄动、太阳光压等多种复杂因素干扰，这些干扰日积月累，会使卫星逐渐偏离预定轨道。一旦轨道偏差超出允许范围，卫星通信信号可能中断、气象监测数据出现偏差、军事侦察精度大打折扣，后果不堪设想。所以，实时、精准预测和跟踪卫星轨道，及时调整卫星姿态与轨道参数，维持卫星稳定运行，成为航天领域亟待攻克的关键课题。

传统轨道预测与跟踪方法：困境与局限

在早期航天探索阶段，科学家主要依靠基于天体力学的经典轨道理论预测和跟踪卫星轨道。这些方法以开普勒定律为基础，结合牛顿万有引力定律，建立卫星运动方程，通过求解方程预测卫星位置。但实际应用中，传统方法暴露出诸多弊端。

复杂环境下，地球非球形引力会使卫星轨道产生复杂摄动，经典理论难以精确描述。太阳光压对卫星轨道的影响，因卫星形状、表面材料特性等因素变得复杂，传统模型无法准确考虑。这导致预测结果与卫星实际轨道偏差逐渐增大，长期预测精度严重不足。并且，传统方法在处理大量观测数据时效率低下。早期地面跟踪站收集的卫星观测数据，需人工整理、分析，计算过程繁琐，数据处理速度远远赶不上卫星轨道变化速度，难以满足实时性要求。

⛳️ 运行结果

📣 部分代码

%from celestrak.com automatically every programmable period of time and

%propagate it with SGP4 code. With this program we can plot multiple satellites

%simultaneously, also we can change many configurations and save it.

%Author: (Tony) Chao Huang

%Major: Electronic Engineer

%From: Santa Cruz - Bolivia

%Email: tony.c.huang@outlook.com (personal public email), chao.huang@abe.bo (only for official contact)

%LinkedIn: https://www.linkedin.com/in/tony-chao-huang-l-6a3a4929

%Current Work: Space Segment Engineer, Bolivian Space Agency - ABE (www.abe.bo)

%Current Study: Master Candidate in Space Technology Application - Remote Sensing and GIS

%BeiHang University, Regional Centre for Space Science and

%Technology Education in Asia and the Pacific (China)(Affiliated to the United Nations) - RCSSTEAP

%Current City: Beijing - China

%January - 2016

function propagation_result=propagate_TLE(sat_cells,min_before_epoch,min_after_epoch,timestep,propagation_mode,wgs_const,leap_second)

%input

%sat_cells: struct with TLE

%min_before_epoch: propagate from minutes before TLE epoch (can be 0 or negative number)

%min_after_epoch: propagate to minutes after TLE epoch

%output

%propagation_result: struct with the propagation results

%base on:

% script testmat.m

% This script tests the SGP4 propagator.

% Author:

% Jeff Beck

% beckja@alumni.lehigh.edu

% Version Info:

% 1.0 (051019) - Initial version from Vallado C++ version.

% 1.0 (aug 14, 2006) - update for paper

% 2.0 (apr 2, 2007) - update for manual operations

% 3.0 (3 jul, 2008) - update for opsmode operation afspc or improved

% 3.1 (2 dec, 2008) - fix tsince/1440.0 in jd update

% these are set in sgp4init

global tumin mu radiusearthkm xke j2 j3 j4 j3oj2

global opsmode

propagation_result = sat_cells;

% // ------------------------ implementation --------------------------

% add operation smode for afspc (a) or improved (i)

% afspc='a', improved= 'i' initl.m

opsmode= propagation_mode; %'i'; %improved mode - Config

whichconst = wgs_const; %84; %WGS-84 getgravc.m - config

% //typerun = 'c' compare 1 year of full satcat data

% //typerun = 'v' verification run, requires modified elm file with

% //typerun = 'm' maunual operation- either mfe, epoch, or dayof yr

% // start stop and delta times

%typerun = 'v'; %v=verification run without input().... for twoline2rv .. sugestion: modiy this function

%we use this mode because we don't need the startmfe, stopmfe, deltamin

%output information from twoline2rv()

typerun = 'c'; %use typerun = 'c' catalog mode (propagates at 20 min timesteps from one day before epoch to one day after)

typeinput='e'; %any value; check twoline2rv.m

rad = 180.0 / pi;

% // ---------------- setup files for operation ------------------

% // input 2-line element set file

% infilename = 'tk.txt';

% infile = fopen(infilename, 'r');

% if (infile == -1)

% fprintf(1,'Failed to open file: %s\n', infilename);

% return;

% end

% outfile = fopen('tmatver.out', 'wt');

error_log=fopen('TLE_Data\progation_error.txt', 'a');

% global idebug dbgfile

% // ----------------- test simple propagation -------------------

total_sat=size(sat_cells,2);

for (sat_n=1:total_sat)

%while (~feof(infile))

% longstr1 = fgets(infile, 130);

% while ( (longstr1(1) == '#') && (feof(infile) == 0) )

% longstr1 = fgets(infile, 130);

% end

% if (feof(infile) == 0)

% longstr2 = fgets(infile, 130);

% if idebug

% catno = strtrim(longstr1(3:7));

% dbgfile = fopen(strcat('sgp4test.dbg.',catno), 'wt');

% fprintf(dbgfile,'this is the debug output\n\n' );

% end

% // convert the char string to sgp4 elements

% // includes initialization of sgp4

%we use typerun = 'c' because we don't need the startmfe, stopmfe, deltamin

%output information from twoline2rv()

[satrec, startmfe, stopmfe, deltamin] = twoline2rv( whichconst, ...

sat_cells{sat_n}.line1, sat_cells{sat_n}.line2, typerun, typeinput);

startmfe = min_before_epoch;

stopmfe = min_after_epoch;

deltamin = timestep;

% fprintf(outfile, '%d xx\n', satrec.satnum);

% fprintf(1,' %d\n', satrec.satnum);

% // call the propagator to get the initial state vector value

[satrec, ro ,vo] = sgp4 (satrec, 0.0);

% fprintf(outfile, ' %16.8f %16.8f %16.8f %16.8f %12.9f %12.9f %12.9f\n',...

% satrec.t,ro(1),ro(2),ro(3),vo(1),vo(2),vo(3));

% fprintf(1, ' %16.8f %16.8f %16.8f %16.8f %12.9f %12.9f %12.9f\n',...

% satrec.t,ro(1),ro(2),ro(3),vo(1),vo(2),vo(3));

%tsince = startmfe;

tsince = startmfe - deltamin;

% // check so the first value isn't written twice

% if ( abs(tsince) > 1.0e-8 )

% tsince = tsince - deltamin;

% end

% // loop to perform the propagation

n=1; % output data line 1 first epoch

output_data=zeros(1,23);

while ((tsince < stopmfe) && (satrec.error == 0))

tsince = tsince + deltamin;

if(tsince > stopmfe)

tsince = stopmfe;

end

[satrec, ro, vo] = sgp4 (satrec, tsince); %ECI\TEME coordinate

if (satrec.error > 0)

fprintf(error_log,'# *** error: t:= %f *** code = %3i\n', tsince, satrec.error);

end

if (satrec.error == 0)

jd = satrec.jdsatepoch + tsince/1440.0;

[year,mon,day,hr,minute,sec] = invjday ( jd );

output_data(n,1:7) = [tsince,ro(1),ro(2),ro(3),vo(1),vo(2),vo(3)]; %change code to improve velocity - preallocate

% fprintf(outfile, ' %16.8f %16.8f %16.8f %16.8f %12.9f %12.9f %12.9f',...

% tsince,ro(1),ro(2),ro(3),vo(1),vo(2),vo(3));

% fprintf(1, ' %16.8f %16.8f %16.8f %16.8f %12.9f %12.9f %12.9f \n',...

% tsince,ro(1),ro(2),ro(3),vo(1),vo(2),vo(3));

[p,a,ecc,incl,node,argp,nu,m,arglat,truelon,lonper] = rv2coe (ro,vo,mu); %ECI\TEME coordinate

output_data(n,8:21) = [ a, ecc, incl*rad, node*rad, argp*rad, nu*rad, m*rad,year,mon,day,hr,minute,sec,jd]; %change code to improve velocity - preallocate

% fprintf(outfile, ' %14.6f %8.6f %10.5f %10.5f %10.5f %10.5f %10.5f %5i%3i%3i %2i:%2i:%9.6f \n',...

% a, ecc, incl*rad, node*rad, argp*rad, nu*rad, m*rad,year,mon,day,hr,minute,sec );

%improve this using iers.org data

%Polar Motion Data December-18 2015

%http://datacenter.iers.org/eop/-/somos/5Rgv/getTX/6/bulletina-xxviii-051.txt

timezone = 0;

dut1 = 0.10536; % ut1 - utc sec

dat = leap_second; %leap second tai - utc sec

lod = 0.0015; %lod - excess length of day sec

xp = 0.0818/3600*pi/180; %xp - polar motion coefficient rad

yp = 0.2496/3600*pi/180; %yp - polar motion coefficient rad

[ut1, tut1, jdut1, utc, tai, tt, ttt, jdtt, tdb, ttdb, jdtdb ] ...