【SRC-Matlab】GNSS 时间系统的转换代码实现

Part.I Introduction

本文分享了常用的时间转换的函数，主要包括判断闰年、年月日 / 年积日 / GPS 周 / 儒列日 互转、时分秒 / 天内秒 / 周内秒 互转… 编程语言为：Matlab。

Chap.I 基础知识

我们常用的表示时间的方式是北京时『时分秒』，即hh:mm:ss；而且北京市并不是UTC，它和UTC存在这样的转换关系：UTC+8小时=北京时

• 判断闰年：年份是 4 的倍数但不是 100 的倍数或者年份是 400 的倍数。
• GPS 周：GPS Week，GPSW GPS系统内部描述时间的一种方式，它的起点为 1980 年 1 月 5 日夜晚与 1980 年 1 月 6 日凌晨之间 0 点（这天是周日，老外经常周日表示一周的第一天，以0表示）。
• 儒列日：Julian Day，JD 是指由公元前 4713 年 1 月 1 日，协调世界时中午 12 时开始所经过的天数。
• 简化儒列日：Modified Julian Day，MJD 后来经国际天文学联合会（1973 年）商讨，采用简化儒列日，其起点是 1858 年 11 月 17 日世界时 0 时；MJD = JD - 2400000.5
• 年积日：Day Of Year，DOY 一年当中的第几天，其取值范围为[1,365]
• 天内秒：Second Of Day，SOD 一天中的第几秒，其取值范围为[1,86400]
• 周内秒：Second Of Week，SOW 一周中的第几秒，其取值范围为[1,604800]
• 周内分：Hour Of Week，HOW 一周中的第几小时，其取值范围为[1,168]
• 协调世界时（Universal Time of Coordination，UTC）：以原子秒长为计量单位，在时刻上与平太阳时之差小于 0.9 秒的时间系统。UT0 是完全按照天体运行计算出来的时间，UT1 是在 UT0 的基础上做了一些调整，UT2 是在 UT0 和 UT1 的基础上又进行了一些调整。天体运动并不是十分稳定的，比如极移和章动的存在。
• 国际原子时（International Atomic Time，TAI）：取 1958 年 1 月 1 日 0 时 0 分 0 秒世界时(UT)的瞬间作为同年同月同日 0 时 0 分 0 秒TAI。
• GPS 时（GPS Time，GPST）：由 GPS 星载原子钟和地面监控站原子钟组成的一种原子时基准，与国际原子时保持有 19s 的常数差，并在 GPS 标准历元 1980 年 1 月 6 日 0 时与 UTC 保持一致。
• 北斗时（BDS Time，BDT） ：同 GPST 一样由原子钟保持基准，在 2006 年 1 月 1 日 0 时与 UTC 保持一致。因为从 1980 年到 2006 年共有 14 次跳秒发生，所以 BDT 和 GPST 相差 14 秒且基本恒定不变。
• GLONASS 时（GLONASS Time，GLST）：以莫斯科本地协调时 UCTsu 定义，其值与 UTC 存在 3 小时时差。
• Galileo 时（Galileo Time，GST）：同 GPST 保持一致。
• 跳秒：leap second，在 1980 年 1 月 6 日 0 时 GPS 时与 UTC 时对齐，GPS 时是依靠稳定的原子钟来维持的，也就是说它的单位时间长度是很稳定的；而协调世界时是根据天文确定的，和地球自转有关，但是地球自转速度在不断变慢，也就是说协调世界时的单位时间长度并不是恒定的。但是总不能他们稍微不一致就调吧，这样也太麻烦了，所以就规定，当两者相差接近1秒时，就让UTC跳一秒。
• 跳秒对是 UCT 的操作，原子时并不会跳秒。并且跳秒通常是把 UTC 『拨快』一秒（因为地球自转速度越来越慢），比如23:59:58下一秒是下一天的00:00:00，它不经过23:59:59；但不排除将其『拨慢』一秒的可能。
• 为什么TAI - GPST = 19？这是因为国际原子时是在 1958 年与 UTC 对齐的，而 GPS 时是在 1980 年与 UTC 对齐的，从 1958 年到 1980 年已经跳了19秒。它们两个的单位时间长度都是靠原子钟维持的，只不过TAI靠的是全球分布的约 240 台原子钟维持的，而 GPS 时是靠数十台原子钟来维持的；但是原子钟的精度本身就很高了，所以它们两个的差值可以说是恒定的。

更多的关于时间系统的介绍可参看博文。

Chap.II 函数总览

下面是涉及的函数总览图。

需要注意的有一下几点：

• 目前关于UTC的跳秒，下面编写的函数并没有考虑到

Part.II Matlab 版代码实现

编写不易，积分多的可给点积分。

Chap.I 完整源码

% Determine if a year is a leap year
function whe=leapyear(year)
whe = 0;
if rem(year,4) == 0 && rem(year,100) ~= 0
    whe = 1;
end
if rem(year,400) == 0
    whe = 1;
end
end

% convert an illegal doy to a normal doy
function [year,doy]=norm_doy(year,doy)
while(1)
    if doy <= 0
        year = year - 1;
        doy = doy + 365 + leapyear(year);
    elseif doy >= 365 + leapyear(year)
        doy = doy - (365 - leapyear(year));
        year = year + 1;
    else
        break;
    end
end
end

% Change year, day of year to year, month, day
function [year,mon,day]=doy2ymd(year, doy)
day = doy;
mon = 0;
monthdays = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
for i=1:13
    monthday = monthdays(i);
    if i==3 && leapyear(year) == 1
        monthday = monthday + 1;
    end
    if day > monthday
        day = day - monthday;
    else
        mon = i - 1;
        break
    end
end
end

% Change year, month, day to year, day of year
function [year,doy]=ymd2doy(year, mon, day)
doy=day;
monthdays = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
for i=1:mon
    doy=doy+monthdays(i);
end
if mon>2
    doy=doy+leapyear(year);
end
end

% Change year, month, day to Modified Julian Day
function mjd=ymd2mjd(year, mon, day)
if mon <= 2
    mon = mon + 12;
    year = year - 1;
end
mjd = 365.25 * year - rem(365.25 * year , 1.0) - 679006.0;
mjd = mjd + floor(30.6001 * (mon + 1)) + 2.0 - ...
    floor(year / 100.0) + floor(year / 400) + day;
end

% Change year, day of year to Modified Julian Day
function mjd=doy2mjd(year, doy)
[year,doy] = norm_doy(year, doy);
[year, mon, day] = doy2ymd(year, doy);
mjd=ymd2mjd(year, mon, day);
end

% Change year,month,day to GPS weeks
function [week,day]=ymd2gpsweek(year, mon, day)
mjd = ymd2mjd(year, mon, day);
week = floor((mjd - 44244.0) / 7.0);
day = floor(mjd - 44244.0 - week * 7.0);
end

% Change Modified Julian Day to year, day of year
% The input date must after 1952.0
function [year,doy]=mjd2ydoy(mjd)
year = 1952;
dd = mjd + 1 - 34012;
yday = 366;
while(dd > yday)
    dd = dd - yday;
    year = year + 1;
    yday = 365 + leapyear(year);
end
year = floor(year);
doy = floor(dd);
end

% Change seconds of day to hours, minutes and seconds
function [hh,mm,ss]=sod2hms(sod)
hh = floor(sod/3600);
mm = floor((sod-hh*3600)/60);
ss = floor(sod - hh * 3600 - mm * 60);
end

% Change hours:minutes:seconds to seconds of day
function sod = hms2sod(hh, mm, ss)
sod = floor(hh)*3600 + floor(mm)*60 + ss;
end

% Change (Modified Julian Day, seconds fo day) to (GPS week, hours of week)
function [week, how]=sod2how(mjd,sod)
[year,doy]=mjd2ydoy(mjd);
[year,mon,day]=doy2ymd(year, doy);
[week,day]=ymd2gpsweek(year, mon, day);
how=(sod+day*86400)/3600.0;
end

Chap.II 调用示例

将上面的源码贴到文件中，然后在文件头部追加下面的内容：

year = 2024; mon = 4; day = 28;
hor=7; min=0; sec=0;
mjd=ymd2mjd(year, mon, day)
sod=hms2sod(hor,min,sec)
doy=ymd2doy(year, mon, day)
[w,d]=ymd2gpsweek(year, mon, day)
[y,m,d]=doy2ymd(2023, 310)

运行得到如下结果

mjd =

       60428


sod =

       25200


doy =

        2024


w =

        2312


d =

     0


y =

        2023


m =

    11


d =

     6