网站链接: element-ui dtcms
当前位置: 首页 > 技术博文  > 技术博文

WGS84 坐标系转到J2000坐标系

2021/5/16 16:58:51 人评论

在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH。 具体转换茹下 step1 WGS 84 转换到协议地球坐标系。 function XYZ wgs842xyz(BLH) %UNTITLED2 大地坐标系转北天东坐标系 % 此处显示详细说明 thita 为当地恒星时…

在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH。 具体转换茹下

step1  WGS 84 转换到协议地球坐标系。 

function XYZ = wgs842xyz(BLH)
%UNTITLED2 大地坐标系转北天东坐标系
%   此处显示详细说明   thita 为当地恒星时
ae  = 6378137; %wgs84坐标系场半轴
ee = 0.00669437999013 ;  %第一偏心率
%ee=0.00669437999; %84系

thitaG = pi/2; %格林尼治面与春分面夹角 
%归一化量纲(地球)

% 地平系坐标
thita = thitaG+BLH(2);
temp =ae/sqrt(1-ee^2*sin(BLH(1)));
temp2 = (temp+BLH(3))*cos(BLH(1));
x = temp2*cos(thita)  ; %东
y = temp2*sin(thita) ;  %北
z = (temp*(1-ee^2)+BLH(3))*sin(BLH(1)) ;%天

CL2I =[-sin(thita) -cos(thita)*sin(BLH(1)) cos(thita)*cos(BLH(1));
       cos(thita) -sin(thita)*sin(BLH(1))  sin(thita)*cos(BLH(1));
       0               cos(BLH(1))              sin(BLH(1))       ];  %地平系转地惯系转换矩阵
XYZ=CL2I * [ x y z]' ;% X 为地心指向春分点,Y为 地心指向赤道面顺向90度,Z为地心指向北


end

step 2  协议地球坐标系 转换为瞬时坐标系 

这主要涉及查询 给定时刻的  xp,yp ,  通过查询IERS公报获得,如果 对计算要求不高,可以忽略 此类型 的 坐标 转换。

earthFixedXYZ=ordinateSingleRotate('x',yp)*ordinateSingleRotate('y',xp)*earthFixedXYZ;

其中  ordinateSingleRotate 为 坐标转换 函数 ,该函数第一个参数为给定的坐标轴,第二个参数为 旋转的角度 ,第三个参数为 旋转的角度的度量单位。 后续要多次调用该函数。

function [R] = ordinateSingleRotate( axis,angle_deg,angleType)
%坐标轴旋转
%   axis 表示围绕旋转的坐标轴 
% '1'或者'x' 表示围绕 x轴逆时针旋转
% '2'或者'y' 表示围绕 y轴逆时针旋转
% '3'或者'z' 表示围绕 x轴逆时针旋转
% angle_deg 表示旋转的角度
% angleType 表示角度类型,可以为 'rad'和’deg‘,默认deg

if nargin()<3
    angleType='deg';
end
switch lower(angleType)
    case 'deg'
    case 'rad'
       angle_deg=angle_deg*180/pi; 
     otherwise
        error(message('ordinateRotate:unknownRotation:angleType=', angleType));  
end
switch axis
    case {'x','X',1,'1'}
        R =[1         0             0  ;...
            0  cosd(angle_deg) sind(angle_deg);...
            0  -sind(angle_deg) cosd(angle_deg) ];
    case {'y','Y',2,'2'}
         R =[cosd(angle_deg)   0  -sind(angle_deg);...
                  0            1     0;...
            sind(angle_deg)   0   cosd(angle_deg) ];
     case {'z','Z',3,'3'}
         R =[cosd(angle_deg)  sind(angle_deg)  0  ;...
            -sind(angle_deg)   cosd(angle_deg)  0;...
                   0              0             1 ];     
    otherwise
        error(message('ordinateRotate:unknownRotation:axis=',axis));
end
end

step 3  瞬时地球坐标系 转换为 瞬时真天球坐标系 

xyz= ordinateSingleRotate('z',-gst_deg)*earthFixedXYZ;

该步骤主要涉及到格林尼治恒星时角的计算,关于格林尼治恒星时角 计算方法 很多,下面给出 一个较为精确的计算方法。

function [gst_deg,JDTDB ] = utc2gst(UTC,dUT1,dAT)

%将utc时间转换为格林尼治恒星时
%参数
%utc 格式为y m d 其中d的数值为小数,将h m s 按(sec/3600+min/60+h)/24转换成小数,并累加到day 上
%dUT1 为ut1-utc 的差,数值不超过正负1秒,查iers可获得数值
%UNTITLED 计算当地恒星时,返回值以时秒为单位
%  UTC世界协调时  协调地球时与原子时的时间差,采用闰秒的办法进行协调
%dAT 润秒数
%TAI国际原子时间 该时间最准  TAI=UTC+dAT;  %国际原子时间
%TT 地球时  TT = TAI+32.184 至2014年的时候;
%TDT 地球动力学时间
% ET 历书时间
%地球时=地球动力学时间=历书时间
J2000=2451545;
if nargin()<2
    dUT1=0;
end
if nargin()<3
    dAT=37.0;
end

JDutc=YMD2JD(UTC(1),UTC(2),UTC(3));
JDUT1  = JDutc+dUT1/86400 ;  %UT1 为世界时,世界时由于自转的不均匀,因此与UTC时间有dUT1的差别,dUT1在各国卫星授时信号中会以0.1秒的精度给出IERS经过处理后,会以1e-5的精度给出。
dT=32.184+dAT-dUT1;
JDTT=JDUT1+dT/86400;


% JDTT= YMD2JD(UTC(1),UTC(2),UTC(3))+dUT1/86400;
%首先计算TDB,TDB是质心动力学时,是太阳月球行星等天体星历表中的时间尺度
T=(JDTT-J2000)/36525;
temp= 0.001657*sin(628.3076*T+6.2401) ...
    +0.000022*sin(575.3385*T+4.2970) ...
    +0.000014*sin(1256.6152*T+6.1969) ...
    +0.000005*sin(606.9777*T+4.0212)...
    +0.000005*sin(52.9691*T+0.4444)...
    +0.000002*sin(21.3299*T+5.5431)...
    +0.00001*T*sin(628.3076*T+4.2490);
JDTDB=JDTT+temp/86400;

%下面计算UT2,UT2是在UT1的基础上修正地球周期性季节变化后得到的世界时间
%       %根据UT1计算Tb,Tb为以贝塞尔年为单位,从历元B2000.0起算ff
%       B2000=2451544.033;
%       Tb=(YMD2JD(UT1(1),UT1(2),UT1(3))-B2000)/365.2422;
%       dTs=0.022*sin(2*pi*Tb)-0.012*cos(2*pi*Tb)-0.006*sin(4*pi*Tb)+0.007*cos(4*pi*Tb);
%       dTs=dTs/86400;
%       UT2=UT1+[0 0 dTs];

%下面计算地球转动角thita,thita 单位为圈,
%下面计算格林尼治平恒星时GMST; 单位为时秒,采用张洪波教材方法
T=(JDTDB-J2000)/36525;
T2=T*T;T3=T2*T;T4=T3*T;T5=T4*T;
Du=JDUT1-J2000;
thita = 0.7790572732640+1.00273781191135448*Du; %自J2000起地球转过的圈数;
GMST_deg=(-0.0000000368*T5-0.000029956*T4-0.00000044*T3+1.3915817*T2+4612.156534*T+0.014506) /3600+(thita-floor(thita))*360  ;

[epthilongA_deg,dertaPthi_deg] = nutationInLongitudeCaculate(JDTDB);

%计算赤经章动引起的误差eclipticObliquitygama,单位为角秒
%计算月球平近点角   太阳平近点角 月球纬度辐角 日月平角距(月球平黄经-太阳平黄经) 月球升交点平黄经 单位为角秒
F_deg =(0.00000417*T4-0.001037*T3-12.7512*T2+1739527262.8478*T+335779.526232)/3600;
F_deg=mod(F_deg,360);
D_deg = (-0.00003169*T4+0.006593*T3-6.3706*T2+1602961601.2090*T+1072260.70369)/3600 ;
D_deg=mod(D_deg,360);
Omiga_deg =(-0.00005939*T4+0.007702*T3+7.4722*T2-6962890.5431*T+450160.398036 )/3600   ;
Omiga_deg=mod(Omiga_deg,360);

epthilongGama_deg=(dertaPthi_deg*cosd(epthilongA_deg)...
    +0.00264096*sind( Omiga_deg )...
    +0.00006352*sind(2*Omiga_deg)...
    +0.00001175*sind(2*F_deg-2*D_deg+3*Omiga_deg)...
    +0.00001121*sind(2*F_deg-2*D_deg+Omiga_deg)...
    -0.00000455*sind(2*F_deg-2*D_deg+2*Omiga_deg)...
    +0.00000202*sind(2*F_deg+3*Omiga_deg)...
    +0.00000198*sind(2*F_deg+Omiga_deg)...
    -0.00000172*sind(3*Omiga_deg)...
    -0.00000087*T*sind(Omiga_deg))/3600;
gst_deg=epthilongGama_deg+GMST_deg;
end

该函数要调用nutationInLongitudeCaculate函数来计算章动交、黄赤交角等 ,该函数在 步骤4中给出;

同时 YMD2JD函数为 年月日转换为儒略日,具体说明 见https://blog.csdn.net/qq_24172609/article/details/112244135?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522162115452716780357279082%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=162115452716780357279082&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~first_rank_v2~rank_v29-3-112244135.first_rank_v2_pc_rank_v29&utm_term=%E5%85%AC%E5%85%83%E7%BA%AA%E5%B9%B4%E6%B3%95&spm=1018.2226.3001.4187

代码如下

    function [jd] =YMD2JD(y,m,d)  %函数公历年月日转儒略日
    if size(y,2)==3  %如果输入为年月日数组
        d=y(3);
        m=y(2);
        y=y(1);        
    end 
if m<3
    m=m+12;
    y=y-1;
end
B=0;
if y>1582||(y==1582&&m>10)||(y==1582&&m==10&&d>=15)
    B=2-floor(y/100)+floor(y/400);
end
jd=floor(365.25*(y+4712))+floor(30.6*(m+1))+d-63.5+B;

    
    
    end

step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系

xyz=ordinateSingleRotate('x',-epthilongA_deg)*ordinateSingleRotate('z',dertaPthi_deg)*ordinateSingleRotate('x',dertaEpthilong_deg)*xyz;
其中  epthilongA_deg 为平黄赤交角, dertaPthi_deg为黄经章动, dertaEpthilong_deg为交角章动, epthilong_deg为瞬时黄赤交角

这些角度的计算方法为

function [epthilongA_deg,dertaPthi_deg,dertaEpthilong_deg,epthilong_deg] = nutationInLongitudeCaculate(JD,accuracy)
%计算平黄赤交角,黄经章动、和交角章动、瞬时黄赤交角。
%   参数
% epthilongA_deg平黄赤交角
% dertaPthi_deg黄经章动
% dertaEpthilong_deg和交角章动
% epthilong_deg瞬时黄赤交角
%JD 儒略日
% accuracy 表示计算的精度要求 数值为’normal‘ 和’high‘ 或者用'N'和’H'表示一般精度和高精度 。默认为高精度计算
if nargin()<2
    accuracy='h';
end
T=(JD-2451545)/36525;
T2=T*T;
T3=T2*T;
T4=T3*T;
T5=T4*T;
%计算月球平近点角lunarMeanAnomaly_deg (l_deg) 太阳平近点角SolarMeanAnomaly_deg(solarl_deg)
%月球纬度辐角lunarLatitudeAngle_deg(F_deg)  日月平角距diffLunarSolarElestialLongitude_deg(D_deg月球平黄经-太阳平黄经)
%月球升交点平黄经SolarAscendingNodeElestialLongitude_deg ( Omiga_deg)
l_deg = (-0.00024470*T4+0.051635*T3+31.8792*T2+1717915923.2178*T+485868.249036)/3600;
l_deg=mod(l_deg,360);
solarl_deg =(-0.00001149*T4+0.000136*T3-0.5532*T2+129596581.0481*T+1287104.79305)/3600;
solarl_deg=mod(solarl_deg,360);
F_deg =(0.00000417*T4-0.001037*T3-12.7512*T2+1739527262.8478*T+335779.526232)/3600;
F_deg=mod(F_deg,360);
D_deg = (-0.00003169*T4+0.006593*T3-6.3706*T2+1602961601.2090*T+1072260.70369)/3600;
D_deg=mod(D_deg,360);
Omiga_deg =(-0.00005939*T4+0.007702*T3+7.4722*T2-6962890.5431*T+450160.398036 )/3600;
Omiga_deg=mod(Omiga_deg,360);
basicAngle_deg=[l_deg solarl_deg F_deg D_deg Omiga_deg];

epthilongA_deg=(-0.0000000434*T5-0.000000576*T4+0.00200340*T3-0.0001831*T2-46.836769*T+84381.406)/3600;
epthilongA_deg=epthilongA_deg-floor(epthilongA_deg/360)*360;
switch lower(accuracy)
    case {'h','high'}    %IAU2000模型有77项
        elestialLonNutation= elestialLonNutationCaculate();
        dertaPthi_deg =-3.75e-08;
        dertaEpthilong_deg =0.388e-3/3600;
         for i = 1:77  %计算日月周期项
             sumAngle_deg=0;
             for j=1:5
                 sumAngle_deg=sumAngle_deg+elestialLonNutation(i,j)*basicAngle_deg(j);
             end
             sumAngle_deg=sumAngle_deg-floor(sumAngle_deg/360)*360;
             
             dertaPthi_deg=dertaPthi_deg+((elestialLonNutation(i,6)+elestialLonNutation(i,7)*T)...
             *sind(sumAngle_deg)+elestialLonNutation(i,8)*cosd(sumAngle_deg))*1e-7/3600;
             dertaEpthilong_deg=dertaEpthilong_deg+((elestialLonNutation(i,9)+elestialLonNutation(i,10)*T)...
             *cosd(sumAngle_deg)+elestialLonNutation(i,11)*sind(sumAngle_deg))*1e-7/3600;
         end
    case{'n','l','normal','low'}
     
                Omiga_deg=basicAngle_deg(5);
        F_deg=basicAngle_deg(3);
        D_deg=basicAngle_deg(4);
        solarl_deg=basicAngle_deg(2);
        dertaPthi_deg=((-17.1996+0.01742*T)*sind(Omiga_deg)...
                       +(0.2062+0.00002*T)*sind(2*Omiga_deg)...
                       -(1.3187+0.00016*T)*sind(2*F_deg-2*D_deg+2*Omiga_deg)...
                       +(0.1426-0.00034*T)*sind(solarl_deg)...
                       -(0.2274+0.00002*T)*sind(2*F_deg+2*Omiga_deg))/3600;
        dertaEpthilong_deg=((9.2025+0.00089*T)*cosd(Omiga_deg)...
                          -(0.0895-0.00005*T)*cosd(2*Omiga_deg)...
                          +(0.5736-0.00031*T)*cosd(2*F_deg-2*D_deg+2*Omiga_deg)...
                          +(0.0054-0.00001*T)*cosd(solarl_deg)...
                          +(0.0977-0.00005*T)*cosd(2*F_deg+2*Omiga_deg) )/3600;
    otherwise
        error(message(':unknownaccuracy:accuracy=',accuracy));
end
epthilong_deg=epthilongA_deg+dertaEpthilong_deg;


end

其中 elestialLonNutationCaculate()  函数 返回一个 77行11列 的 二维数组,函数 如下

function elestialLonNutation=elestialLonNutationCaculate()                              

elestialLonNutation = ...
  [0 0 0 0 1 -1.72064161E+8 -174666 33386 9.2052331E+7 9086 15377;
   0 0 2 -2 2 -1.3170906E+7 -1675 -13696 5.730336E+6 -3015 -4587;
   0 0 2 0 2 -2.276413E+6 -234 2796 978459 -485 1374;
   0 0 0 0 2 2.074554E+6 207 -698 -897492 470 -291;
   0 1 0 0 0 1.475877E+6 -3633 11817 73871 -184 -1924;
   0 1 2 -2 2 -516821 1226 -524 224386 -677 -174;
   1 0 0 0 0 711159 73 -872 -6750 0 358;
   0 0 2 0 1 -387298 -367 380 200728 18 318;
   1 0 2 0 2 -301461 -36 816 129025 -63 367;
   0 -1 2 -2 2 215829 -494 111 -95929 299 132;
   0 0 2 -2 1 128227 137 181 -68982 -9 39;
   -1 0 2 0 2 123457 11 19 -53311 32 -4;
   -1 0 0 2 0 156994 10 -168 -1235 0 82;
   1 0 0 0 1 63110 63 27 -33228 0 -9;
   -1 0 0 0 1 -57976 -63 -189 31429 0 -75;
   -1 0 2 2 2 -59641 -11 149 25543 -11 66;
   1 0 2 0 1 -51613 -42 129 26366 0 78;
   -2 0 2 0 1 45893 50 31 -24236 -10 20;
   0 0 0 2 0 63384 11 -150 -1220 0 29;
   0 0 2 2 2 -38571 -1 158 16452 -11 68;
   0 -2 2 -2 2 32481 0 0 -13870 0 0;
   -2 0 0 2 0 -47722 0 -18 477 0 -25;
   2 0 2 0 2 -31046 -1 131 13238 -11 59;
   1 0 2 -2 2 28593 0 -1 -12338 10 -3;
   -1 0 2 0 1 20441 21 10 -10758 0 -3;
   2 0 0 0 0 29243 0 -74 -609 0 13;
   0 0 2 0 0 25887 0 -66 -550 0 11;
   0 1 0 0 1 -14053 -25 79 8551 -2 -45;
   -1 0 0 2 1 15164 10 11 -8001 0 -1;
   0 2 2 -2 2 -15794 72 -16 6850 -42 -5;
   0 0 -2 2 0 21783 0 13 -167 0 13;
   1 0 0 -2 1 -12873 -10 -37 6953 0 -14;
   0 -1 0 0 1 -12654 11 63 6415 0 26;
   -1 0 2 2 1 -10204 0 25 5222 0 15;
   0 2 0 0 0 16707 -85 -10 168 -1 10;
   1 0 2 2 2 -7691 0 44 3268 0 19;
   -2 0 2 0 0 -11024 0 -14 104 0 2;
   0 1 2 0 2 7566 -21 -11 -3250 0 -5;
   0 0 2 2 1 -6637 -11 25 3353 0 14;
   0 -1 2 0 2 -7141 21 8 3070 0 4;
   0 0 0 2 1 -6302 -11 2 3272 0 4;
   1 0 2 -2 1 5800 10 2 -3045 0 -1;
   2 0 2 -2 2 6443 0 -7 -2768 0 -4;
   -2 0 0 2 1 -5774 -11 -15 3041 0 -5;
   2 0 2 0 1 -5350 0 21 2695 0 12;
   0 -1 2 -2 1 -4752 -11 -3 2719 0 -3;
   0 0 0 -2 1 -4940 -11 -21 2720 0 -9;
   -1 -1 0 2 0 7350 0 -8 -51 0 4;
   2 0 0 -2 1 4065 0 6 -2206 0 1;
   1 0 0 2 0 6579 0 -24 -199 0 2;
   0 1 2 -2 1 3579 0 5 -1900 0 1;
   1 -1 0 0 0 4725 0 -6 -41 0 3;
   -2 0 2 0 2 -3075 0 -2 1313 0 -1;
   3 0 2 0 2 -2904 0 15 1233 0 7;
   0 -1 0 2 0 4348 0 -10 -81 0 2;
   1 -1 2 0 2 -2878 0 8 1232 0 4;
   0 0 0 1 0 -4230 0 5 -20 0 -2;
   -1 -1 2 2 2 -2819 0 7 1207 0 3;
   -1 0 2 0 0 -4056 0 5 40 0 -2;
   0 -1 2 2 2 -2647 0 11 1129 0 5;
   -2 0 0 0 1 -2294 0 -10 1266 0 -4;
   1 1 2 0 2 2481 0 -7 -1062 0 -3;
   2 0 0 0 1 2179 0 -2 -1129 0 -2;
   -1 1 0 1 0 3276 0 1 -9 0 0;
   1 1 0 0 0 -3389 0 5 35 0 -2;
   1 0 2 0 0 3339 0 -13 -107 0 1;
   -1 0 2 -2 1 -1987 0 -6 1073 0 -2;
   1 0 0 0 2 -1981 0 0 854 0 0;
   -1 0 0 1 0 4026 0 -353 -553 0 -139;
   0 0 2 1 2 1660 0 -5 -710 0 -2;
   -1 0 2 4 2 -1521 0 9 647 0 4;
   -1 1 0 1 1 1314 0 0 -700 0 0;
   0 -2 2 -2 1 -1283 0 0 672 0 0;
   1 0 2 2 1 -1331 0 8 663 0 4;
   -2 0 2 2 2 1383 0 -2 -594 0 -2;
   -1 0 0 0 2 1405 0 4 -610 0 2;
   1 1 2 -2 2 1290 0 0 -556 0 0];
end

step5  瞬时平天球坐标系转换为协议天球坐标系(J2000)

xyz=ordinateSingleRotate('Z',zetaA)*ordinateSingleRotate('y',-thitaA)*ordinateSingleRotate('z',zA)*xyz;

该步主要涉及到 岁差角的计算。下面给出计算方法。

function [zetaA,thitaA,zA] = precessionAngle(JDTDB)
%zA、thitaA、zetaA为赤道岁差角,计算赤道岁差角
T=(JD-2451545)/36525;
T2=T*T;
T3=T2*T;
%张玉祥 人造卫星测轨技术方法
% zetaA=(2306.218*T+0.30188*T2+0.017998*T3)/3600;
% zA=(2306.218*T+1.09468*T2+0.018203)/3600;
% thitaA=(2004.3109*T-0.42665*T2-0.041833*T3)/3600;
T4=T3*T;
T5=T4*T;
%张洪波  
zetaA=(-0.0000003173*T5-0.000005971*T4+0.01801828*T3+0.2988499*T2+2306.083227*T+2.650545)/3600;
thitaA=(-0.0000001274*T5-0.000007089*T4-0.04182264*T3-0.4294934*T2+2004.191903*T)/3600;
zA=(-0.0000002904*T5-0.000028596*T4+0.01826837*T3+1.0927348*T2+2306.077181*T-2.650545)/3600;

end

函数中 JDTDB 为给定时刻 的地球动力学时对应的儒略日,其计算方法由步骤三种的函数给出。

相关资讯

    暂无相关的数据...

共有条评论 网友评论

验证码: 看不清楚?