GNSS 时间系统的转换代码实现（Matlab/Python）

2023-03-10 04:09| 来源: 网络整理| 查看: 265

文章目录 Part.I IntroductionChap.I 基础知识Chap.II 函数总览 Part.II Python 版代码实现Part.III Matlab 版代码实现

Part.I Introduction 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 Python 版代码实现

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

import math import time monthdays = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] def leapyear(year): """ Determine if a year is a leap year """ if year % 4 == 0 and year % 100 != 0: return 1 if year % 400 == 0: return 1 return 0 def norm_doy(year, doy): """ Determine if a doy is legal """ while True: if doy 365 + leapyear(year): doy -= (365 + leapyear(year)) year += 1 else: break return year, doy def doy2ymd(year, doy): """ Change year,day of year to year,month,day """ day = doy mon = 0 for i in range(13): monthday = monthdays[i] if i == 2 and leapyear(year) == 1: monthday += 1 if day > monthday: day -= monthday else: mon = i break return mon, day def ymd2doy(year, mon, day): """ Change year,month,day to year,day of year """ doy = day for i in range(1, mon): doy += monthdays[i] if mon > 2: doy += leapyear(year) return doy def ymd2mjd(year, mon, day): """ Change year,month,day to Modified Julian Day """ mjd = 0.0 if mon yday: dd -= yday year += 1 if leapyear(year) == 0: yday = 365 else: yday = 366 return int(dd), int(year) def mjd2ymd(mjd): """ Change year-mon-day to Modified Julian Day. """ doy, year = mjd2ydoy(mjd) mon, day = doy2ymd(year, doy) return year, mon, day def sod2hms(sod): """ Change seconds of day to hours, minutes and seconds """ sod = float(sod) hh = math.floor(sod / 3600) mm = math.floor((sod - hh * 3600) / 60.0) ss = int(sod - hh * 3600 - mm * 60) return hh, mm, ss def hms2sod(hh, mm=0, ss=0.0): """ Change hours:minutes:seconds to seconds of day """ return int(hh) * 3600 + int(mm) * 60 + float(ss) def sod2how(mjd, sod): """ Change (Modified Julian Day, seconds fo day) to (GPS week, hours of week) """ """ not consider second slip now """ doy, year = mjd2ydoy(mjd) mon, day = doy2ymd(year, doy) week, day = ymd2gpsweek(year, mon, day) return week, (sod + day * 86400) / 3600.0 def sod2sow(mjd, sod): """ Change (Modified Julian Day, seconds fo day) to (GPS week, seconds of week) """ """ not consider second slip now """ doy, year = mjd2ydoy(mjd) mon, day = doy2ymd(year, doy) week, day = ymd2gpsweek(year, mon, day) return week, sod + day * 86400 def sow2hms(sow): """ Change seconds of week to (day of week, hours, minutes and seconds) """ """ not consider second slip now """ sod = sow % (24*3600) dow = math.floor(sow/(24*3600)) + 1 return dow, sod2hms(sod) def sow2sod(sow): """ convert Second of Week to Second of Day """ return sow % (24*3600) Part.III Matlab 版代码实现

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

% 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 = 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 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

【本文地址】

公司简介

联系我们