罗兰关于MATLAB的艺术

将想法转化为MATLAB

到底什么是闰秒呢?

发布的罗兰美国舒尔，

29次(过去30天 喜欢| 9日评论
今天的帖子来自Peter Perkins，他是MathWorks开发团队的成员。
在做过MATLAB的一些时间和日期函数后，MathWorks的人有时会问我一些关于日历和计时的问题。通常情况是这样的:1月初(或7月初)，我在午餐时盯着手机一会儿，脱口而出:“ 公报C 发表!明年6月(或12月)没有新的闰秒了! 再一次 ，有些人会说:“闰秒到底是什么东西?”我不是这方面的专家，但我的标准答案是这样开始的:“所以，事实证明地球正在减速。”接下来是故事的其余部分。
在此过程中，我将展示一些对处理时间序列数据有用的时间表和日期时间的用法。

所以，事实证明地球正在减速……但不是最近

在日常生活中，我们认为一天的长度是太阳回到天顶所需要的时间，换句话说，地球绕它的轴旋转一次。我们认为这是24小时，但事实上它是不同的(作为结果 开普勒第二定律 )，由于地球轨道的偏心(开普勒第一定律)，季节相差约30秒/-21秒。那将是很不方便的长度随“天”在过去的一年,因为如果我们认为每天24 * 60 * 60 = 86400秒长,然后第二个需要不同的长度比6月12月,和测量时间恰恰是不可能的。这些对太阳日长度的季节性影响可以在一年中被消除，导致“平均太阳日”，它的长度在一年内不会变化。以平均太阳日为基础的时间系统被称为“世界时”(UT)。
UT的问题是，除了偏心轨道造成的“几何效应”之外，地球的自转速度并不随时间而恒定。事实上，由于潮汐与月球的相互作用，自转有一个长期的放缓，据估计，在过去6亿年里，天长了几个小时。还有一个短期因素，规模较小，但方向相反，据信是由于近几千年来大陆冰盖的融化。从历史上看，太阳自转的减慢意味着太阳日的长度平均每世纪增加约1.4-1.7毫秒。还有季节性循环，以及地质事件和其他原因造成的不可预测的波动。在任何给定的几十年的时间里，即使是一个 的意思是 太阳日可能增加，减少，或两者兼有。
出于科学目的，第二个在19世纪被定义为 平均太阳日的1/86400 ，但如果一天的长度随时间而波动，那么一秒实际上等于多少?这个定义是基于18世纪和19世纪收集的天文数据，有效地利用了1820年左右的平均太阳日。但即使到了19世纪末，人们也认识到当时的平均太阳日略长86400秒。到20世纪40年代和50年代，石英和原子钟使地球自转速度的变化更加明显，并导致 更精确的定义 一个国际单位制(SI)秒的长度。不足为奇的是，地球的自转没有注意到这一点，在20世纪，平均太阳日的长度继续上下波动。最近几十年，它大多比86400秒长1-3毫秒。这一差异被称为“超长日(LOD)” 国际地球自转和参考系统服务中心 ．
文件=“https://datacenter.iers.org/data/latestVersion/EOP_14_C04_IAU2000A_one_file_1962-now.txt”；
IERSdata = readtable(文件,“NumHeaderLines”14);
IERSdata.Properties。变量名([4 8])= [“MJD”“ExcessLOD”]；
IERSdata。日期= datetime(IERSdata.MJD,“ConvertFrom”，“mjd”)；%修改儒略日期->日期时间
IERSdata。ExcessLOD= seconds(IERSdata.ExcessLOD);
IERSdata = table2timetable (IERSdata (:,ExcessLOD“日期”)));
情节(IERSdata.Date IERSdata.ExcessLOD,“b -”)；
ylabel (“超额LOD”)；

好吧，当然，但是是是一闰秒吗?

几个额外的毫秒似乎不重要，但过多的LOD会累积。一个将“一天”定义为86400国际单位秒的时钟实际上相对于地球自转运行得很快，相对于太阳每天获得时间。例如，太阳在6月21日出现在天顶的时钟上的时间每年都要晚零点几秒。不过，把20世纪60年代以来的过量LOD加起来，就会发现这种变化还不到一分钟，甚至在日常生活中都不明显，所以谁在乎呢?
IERSdata。CumulativeELOD = [0;cumsum (IERSdata.ExcessLOD (1: end-1)];
stackedplot (IERSdata [“ExcessLOD”“CumulativeELOD”]）;
但是经过几个世纪的积累，LOD累积到中午太阳明显不再出现在天顶附近。换句话说，基于每天86400 SI秒的计时将偏离物理经验，即UT。这被认为是一个问题,尽管显而易见的事实,中午的太阳位置的变化城市之间在不同时区的“结束”,甚至在标准时间和夏令时在同一个城市,远远大于在上个世纪的转变。1972年之前, 各种特殊调整 保持民用时间与UT大致同步。到1972年, 时间系统的当前版本 当时制定了协调世界时(UTC)，它使用原子钟定义的86400 SI秒，但偶尔会插入额外的“闰秒”调整，以保持UTC与UT ( UT1 更准确地说)。想法是,因为86400年代有点太短对地球的自转,UTC时钟运行有点太快了对太阳的开销,所以你只需要慢下来,使它列举一个额外的第二个偶尔。也有可能地球的自转在短期内会加速，使86400秒也增加一点 长 在这种情况下，闰秒可能最终需要 删除 ．到目前为止，这种情况从未发生过;下面有更多相关内容。
地球的自转是如此不可预测，以至于 而闰秒只会提前6个月公布 在一月和七月。这种特殊的时间安排使得在每个MATLAB版本中获取最新信息变得非常困难。我们试(检查第二个输出 leapseconds 功能)。

等一下!

UTC的当前定义是从1972年1月1日午夜开始采用的，所以我将只选择1971年之后的数据。
IERSdata1972 = IERSdata (timerange (“1 - 1月- 1972”正)“ExcessLOD”)；
20世纪60年代对民用计时进行的特别调整意味着在1972年1月1日UTC与UT1相差45毫秒。但是，如果1972年之后没有闰秒调整，累积的过量LOD将导致UTC在多年后远离UT1。到2017年，差距将扩大到约26秒。
DiffUT1_1972 =秒(0.0454859);
IERSdata1972。UnadjDiff= DiffUT1_1972 + [0; cumsum(IERSdata1972.ExcessLOD(1:end-1))];
IERSdata1972 (timerange (”29日- 12月- 2016，“4 - 1月- 2017”):)
ans = 6×2时间表
日期 ExcessLOD UnadjDiff
1 29日- 12月- 2016 0.0008055秒 26.402秒
2 2016年- 12月30日 0.0008525秒 26.402秒
3. 2016年- 12月31日 0.0009173秒 26.403秒
4 01 - 1月- 2017 0.001016秒 26.404秒
5 02 - 1月- 2017 0.0011845秒 26.405秒
6 03 - 1月- 2017 0.0013554秒 26.406秒
这种不断增长的差异正是闰秒的目的所在。的 leapseconds 函数枚举自1972年以来发生的每个“闰秒事件”。由 leapseconds 包含时间戳和事件的描述(到目前为止，所有的闰秒都是“正的”插入，而不是“负的”删除)，以及截至事件和包括该事件的累计闰秒数。
lsEvents = leapseconds ()
lsEvents = 27×2时间表
日期 类型 CumulativeAdjustment
1 1972年- 6月30日 + 1秒
2 1972年- 12月31日 + 2秒
3. 1973年- 12月31日 + 3秒
4 1974年- 12月31日 + 4秒
5 1975年- 12月31日 + 5秒
6 1976年- 12月31日 + 6秒
7 1977年- 12月31日 + 7秒
8 1978年- 12月31日 + 8秒
9 1979年- 12月31日 + 9秒
10 1981年- 6月30日 + 10秒
11 1982年- 6月30日 + 11秒
12 1983年- 6月30日 + 12秒
13 1985年- 6月30日 + 13秒
14 1987年- 12月31日 + 14秒
使用闰秒时间戳索引包含LOD数据的时间表，我将添加如果没有插入闰秒，在每个瞬间UTC和UT1之间未调整的差值。
lsEvents。UnadjDiff = IERSdata1972.UnadjDiff (lsEvents.Date)
lsEvents = 27×3的时间表
日期 类型 CumulativeAdjustment UnadjDiff
1 1972年- 6月30日 + 1秒 0.63494秒
2 1972年- 12月31日 + 2秒 1.1864秒
3. 1973年- 12月31日 + 3秒 2.2976秒
4 1974年- 12月31日 + 4秒 3.289秒
5 1975年- 12月31日 + 5秒 4.2718秒
6 1976年- 12月31日 + 6秒 5.3334秒
7 1977年- 12月31日 + 7秒 6.3469秒
8 1978年- 12月31日 + 8秒 7.398秒
9 1979年- 12月31日 + 9秒 8.3526秒
10 1981年- 6月30日 + 10秒 9.6281秒
11 1982年- 6月30日 + 11秒 10.389秒
12 1983年- 6月30日 + 12秒 11.249秒
13 1985年- 6月30日 + 13秒 12.451秒
14 1987年- 12月31日 + 14秒 13.634秒
查看每个闰秒事件的未调整差异，很明显，只要额外的LOD积累了大约1秒，就会插入闰秒。将闰秒绘制为LOD数据上覆盖的点，从视觉上证实了这一点。
情节(IERSdata1972.Date IERSdata1972.UnadjDiff);
ylabel (“累积过多的LOD”)；
持有在
情节(lsEvents.Date IERSdata1972.UnadjDiff (lsEvents.Date),“r”。)；
持有从

闰秒为何有效

每一闰秒都是在特定时间发生的瞬间事件。另一种表示它们的方法是将变量附加到LOD数据中，表示任意给定日期的累计闰秒之和。闰秒发生在一天的结束，所以首先给每个事件的日期之后的一天分配1秒的值…
IERSdata1972.Adjustment (lsEvents.Date + caldays(1)) =秒(1);
IERSdata1972 (timerange (”29日- 12月- 2016，“4 - 1月- 2017”):)
ans = 6×3的时间表
日期 ExcessLOD UnadjDiff 调整
1 29日- 12月- 2016 0.0008055秒 26.402秒 0秒
2 2016年- 12月30日 0.0008525秒 26.402秒 0秒
3. 2016年- 12月31日 0.0009173秒 26.403秒 0秒
4 01 - 1月- 2017 0.001016秒 26.404秒 1秒
5 02 - 1月- 2017 0.0011845秒 26.405秒 0秒
6 03 - 1月- 2017 0.0013554秒 26.406秒 0秒
．．．然后用累积的和来填充 调整 在数据的所有其他时间。
IERSdata1972。调整= cumsum(IERSdata1972.Adjustment);
IERSdata1972 (timerange (”29日- 12月- 2016，“4 - 1月- 2017”):)
ans = 6×3的时间表
日期 ExcessLOD UnadjDiff 调整
1 29日- 12月- 2016 0.0008055秒 26.402秒 26秒
2 2016年- 12月30日 0.0008525秒 26.402秒 26秒
3. 2016年- 12月31日 0.0009173秒 26.403秒 26秒
4 01 - 1月- 2017 0.001016秒 26.404秒 27秒
5 02 - 1月- 2017 0.0011845秒 26.405秒 27秒
6 03 - 1月- 2017 0.0013554秒 26.406秒 27秒
这个新变量每次都在数据中定义，因此可以使用它来计算UTC和UT1之间的实际差值(给定所有的闰秒调整)，方法是将它从上面计算的未调整差值中减去。按照惯例，它被计算为UT1-UTC，与之相反的符号是我们所称的未调整的差异，表示为DUT1。
IERSdata1972。DUT1= -(IERSdata1972.UnadjDiff - IERSdata1972.Adjustment);
IERSdata1972 (timerange (”29日- 12月- 2016，“4 - 1月- 2017”):)
ans = 6×4时间表
日期 ExcessLOD UnadjDiff 调整 DUT1
1 29日- 12月- 2016 0.0008055秒 26.402秒 26秒 -0.40159秒
2 2016年- 12月30日 0.0008525秒 26.402秒 26秒 -0.40239秒
3. 2016年- 12月31日 0.0009173秒 26.403秒 26秒 -0.40324秒
4 01 - 1月- 2017 0.001016秒 26.404秒 27秒 0.59584秒
5 02 - 1月- 2017 0.0011845秒 26.405秒 27秒 0.59482秒
6 03 - 1月- 2017 0.0013554秒 26.406秒 27秒 0.59364秒
闰秒调整的目的是保持UTC与UT1大致同步。策划 调整 随着时间的推移，分段常数步长函数说明了调整是如何工作的:步长函数大致逼近累积的多余LOD，因此减去它会保持DUT1较小。沿着同一时间轴绘制DUT1，确认调整保持差值小于1秒。
stackedplot (IERSdata1972 {[“UnadjDiff”“调整”］“DUT1”}）;

上升时间和下降时间

只有在需要时才会插入闰秒。看看上面的一些数据，跳跃秒之间的平均时间似乎自1972年以来一直在增加。计算每十年发生的事件的数量就可以证明这一点。
十年=离散化(lsEvents。目前为止,“十年”，“分类”)；
柱状图(十年);
为了了解原因，我们将回到多余的LOD数据。
自20世纪60年代以来，曾出现过几个短期内过量LOD下降的时期。平滑原始LOD数据可以更容易地看到局部行为。
IERSdata。SmoothedELOD = smoothdata(秒(IERSdata.ExcessLOD),“黄土”，“SmoothingFactor”、。4);
情节(IERSdata.Date IERSdata.ExcessLOD,“b -”)；
持有在
情节(IERSdata.Date IERSdata.SmoothedELOD,“g -”，“线宽”2);
持有从
ylabel (“超额LOD”)；
从平滑的数据中，识别短期趋势改变方向的波峰和波谷，并得到每个波谷的日期和平滑后的超额LOD值。
山峰=找到(islocalmax (IERSdata.SmoothedELOD));
波谷=找到(islocalmin (IERSdata.SmoothedELOD));
日期= IERSdata.Date((波峰、波谷));
值= IERSdata.SmoothedELOD((波峰、波谷));
接下来，我将创建一个包含高峰/低谷事件的时间表，并在情节中标记每个高峰/低谷事件。时间表的每一行都包含事件日期、事件类型和该事件的超额LOD值。
类型=分类([0(大小(山峰));(大小(波谷))],[0 1]、[“峰”“槽”]）;
extremaEvents = sortrows(时间表(日期、类型、值))
extremaEvents = 9×2时间表
日期 类型 价值
1 24 - 7 - 1972 峰 0.0030
2 1975年- 7月31日 槽 0.0028
3. 25 - 6 - 1977 峰 0.0029
4 22日- 1月- 1987 槽 0.0013
5 10 - 11月- 1993 峰 0.0023
6 11月20 - - 2003 槽 0.0003
7 10 - 2月- 2008 峰 0.0009
8 25 - 7 - 2010 槽 0.0007
9 28日- 12月- 2015 峰 0.0012
持有在
isPeak = (extremaEvents。类型= =“高峰”)；
情节(extremaEvents.Date (isPeak) extremaEvents.Value (isPeak),“y ^”，“MarkerFaceColor”，“y”)；
情节(extremaEvents.Date (~ isPeak) extremaEvents.Value (~ isPeak),“青年志愿”，“MarkerFaceColor”，“y”)；
持有从
上面是每一天的超额LOD的图，用一个平滑的版本粗略地跟踪了每年的平均值。它说的是，一年中白天的平均长度已经变得差不多等于，甚至 略低于 这意味着UTC-UT1之间的差异实际上是持平或缩小的，而不是增长的。前面所示的累积过量LOD曲线表明它已经趋于平稳。这就解释了为什么自2016年12月以来没有增加任何正闰秒。一个缩小的UTC-UT1增加了A 负 “闰秒”最终还是需要的。还没有，耐心等待，但如果这种趋势继续下去，可以想象，多余的LOD会在几年内累积接近-1和负闰秒 将 然后根据UTC的定义是必要的。这在以前从来没有发生过，并且可能会造成很多软件很多的混淆(MATLAB不会被混淆)。一个消极的闰秒可能就是动机 终于放弃闰秒 ．然而，许多专家认为，这种趋势不会持续下去，而只是周期的一部分。
这个数字里有很多东西。这里有多余的LOD(本质上是太阳日长)，它的斜率(日长变化的速度)和它的积分(UTC和UT1之间的差)。

与时间赛跑

从峰谷到谷谷，过量LOD逐渐减小，说明地球自转速度加快。我将在一个表中存储这些递减周期的开始和结束。最后一个峰值没有匹配的波谷;我将使用数据中的最后一次时间作为最后一个间隔的结束时间。
StartDate可以= IERSdata.Date(峰);
EndDate = [IERSdata.Date(槽);IERSdata.Date(结束)];
decreasingEvents =表(StartDate可以EndDate)
decreasingEvents = 5×2表
StartDate可以 EndDate
1 24 - 7 - 1972 1975年- 7月31日
2 25 - 6 - 1977 22日- 1月- 1987
3. 10 - 11月- 1993 11月20 - - 2003
4 10 - 2月- 2008 25 - 7 - 2010
5 28日- 12月- 2015 06 - 4月- 2021
持有在
(decreasingEvents startEnd =。StartDate可以decreasingEvents。EndDate]；
h = fill(startEnd(:，[1 2 2 1])，[-。002 -。002．005．005］，“红色”，“FaceAlpha”2,“线型”，“没有”)；
持有从
最后，我将计算每个间隔期间超额LOD的平均减少量(以每年“每日超额LOD的秒数”为单位)，并将该信息添加到表中的间隔事件中。
dTime = decreasingEvents。EndDate- decreasingEvents.StartDate;
dExcess = IERSdata.SmoothedELOD(supplingevents . enddate) - IERSdata.SmoothedELOD(supplingevents . startdate);
decreasingEvents。annualavgreduce = seconds(dExcess ./ years(dTime))
decreasingEvents = 5×3表
StartDate可以 EndDate AnnualAvgDecrease
1 24 - 7 - 1972 1975年- 7月31日 -8.6974 e-05秒
2 25 - 6 - 1977 22日- 1月- 1987 -0.00016705秒
3. 10 - 11月- 1993 11月20 - - 2003 -0.0001993秒
4 10 - 2月- 2008 25 - 7 - 2010 -5.6298 e-05秒
5 28日- 12月- 2015 06 - 4月- 2021 -0.00030701秒
换句话说，在过去几年里，平均太阳日，也就是一整年的平均太阳日，以每年0.3毫秒的速度减少。正如我们上面看到的，它目前非常接近86400。如果最近的减少持续下去，一年后平均太阳日比86400秒小0.3毫秒，这意味着负的过量LOD将在每年大约十分之一秒的时间内积累。但按照这个速度，需要近十年的时间才能实现负闰秒。

刻不容缓

地球会停止加速吗?从长远来看，是的。你不可能战胜潮汐的相互作用。
在近年来的许多时期，有一个很短的时期，原始过剩LOD是显著负的，相对而言。这些只是短暂的波动，但在这些波动期间，地球的自转速度高达1.5毫秒 快 超过86400 SI秒。
情节(IERSdata.Date IERSdata.ExcessLOD,“b -”)；
ylabel (“超额LOD”)；
持有在
线(IERSdata。日期([1 end]),[0 0],“颜色”，“k”，“线型”，”:“）
持有从
ylabel (“超额LOD”)；
我将确定哪一年的额外LOD在任何一天都小于-0.5ms，找出每一年的最小额外LOD，并创建这些最小值的表。
negLOD = IERSdata (IERSdata。ExcessLOD< milliseconds(-0.5),“ExcessLOD”)；
groupsummary (negLOD“日期”，“年”，“最小值”）
ans = 11×3表
year_Date GroupCount min_ExcessLOD
1 2001 8 -0.0007064秒
2 2002 13 -0.0007436秒
3. 2003 31 -0.0009769秒
4 2004 29 -0.0010672秒
5 2005 25 -0.0010809秒
6 2007 4 -0.0006192秒
7 2010 12 -0.000784秒
8 2018 6 -0.0006457秒
9 2019 19 -0.0009571秒
10 2020 88 -0.0014663秒
11 2021 11 -0.0007026秒
那么明年会发生什么呢?闰秒之所以只提前6个月公布，是因为很难准确预测地球自转速度的变化。但是,尽管 日前发表的 并不能预测未来会有过多的LOD 发表DUT1预测 出去一年。
t = readtable (“https://datacenter.iers.org/data/latestVersion/finals.all.iau2000.txt”)；
T = T (:，[3 10 12]);
t.Properties.VariableNames = [“日期”“DUT1”“LOD”]；
t.Date = datetime (t。目前为止,“ConvertFrom”，“mjd”)；
t.DUT1 =秒(t.DUT1);
t.LOD =毫秒(t.LOD);
t = t (t。日期> =“1 - 1月- 2020”:)
t = 908×3表
日期 DUT1 LOD
1 01 - 1月- 2020 -0.17716秒 0.0004379秒
2 02 - 1月- 2020 -0.17763秒 0.0004915秒
3. 03 - 1月- 2020 -0.17811秒 0.0004743秒
4 04 -简- 2020 -0.17859秒 0.0004946秒
5 05 - 1月- 2020 -0.17908秒 0.0004482秒
6 06 - 1月- 2020 -0.17943秒 0.0002298秒
7 07 - 1月- 2020 -0.17955秒 4.14 e-05秒
8 08 - 1月- 2020年 -0.17953秒 -8.77 e-05秒
9 09 - 1月- 2020 -0.17938秒 -0.0002198秒
10 10 - 1月- 2020 -0.17912秒 -0.0002715秒
11 11 - 1月- 2020 -0.17888秒 -0.0001832秒
12 12 - 1月- 2020 -0.17879秒 2.29 e-05秒
13 13 - 1月- 2020 -0.17894秒 0.0002682秒
14 14 - 1月- 2020 -0.17933秒 0.0005166秒
预测= (t.Date >=“1 - 4月- 2021)；
情节(t.Date(~)预测,t.DUT1预测(~),“b”，．．．
t.Date(预测),t.DUT1(预测),“r”)；
ylabel (“DUT1”）
传奇([“实际”“预测”]，“位置”，“西北”）
请记住，DUT1为UT1-UTC，当DUT1小于0.5s时将增加正闰秒。从过去一年的斜率可以看出，太阳日长度先是大于86400秒，然后小于86400秒，最后大致等于86400秒。下一年的斜率预测太阳日长度将小于86400秒，然后等于86400秒，然后大于86400秒。如果正确，这意味着……一段时间内没有闰秒，消极或积极!在接下来的几年里，关注这个问题将会很有趣。如果自1970年以来的模式继续，则为负闰秒 将 是需要的，但问题是什么时候:过量的LOD会很快再次出现，作为最近十年周期的一部分，还是会像2015年以来那样继续减少?

阅读时间

我对时间系统和闰秒的大部分了解都来自于此 这些非常全面的页面 史蒂夫·艾伦 UCO /利克天文台 ．其他关于UTC和闰秒的历史的非常可读的叙述也在后面 莱文 ， 麦卡锡 ， 奎因 , 纳尔逊等 ．

该结束了

当你计算跨越一个或多个闰秒的两个瞬间之间的差值时，闰秒就会显得很丑。在MATLAB中，你可以选择使用闰秒 datetime 计算通过使用 UTCLeapSeconds 创建日期时间:
datetime(2016, 12日31,36岁,0,0,“时区”，“UTCLeapSeconds”31) - datetime(2016年,12日,12日,0,0,“时区”，“UTCLeapSeconds”）
ans =持续时间24:00:01
但我们发现，大多数人只会对随机出现的额外秒感到恼火，所以无需明确要求闰秒计算，你就会看到
datetime(2016, 12日31,36岁,0,0,“时区”，“UTC”31) - datetime(2016年,12日,12日,0,0,“时区”，“UTC”）
ans =持续时间24:00:00
其他时区也一样，比如 美国/ New_York ．(我们还发现，许多人觉得一年两次多出一个小时或少出一个小时很烦人，所以你也必须选择使用时区来使用日光节约时间。)
你的工作在计算中需要闰秒吗?你用过闰秒吗 UTCLeapSeconds 吗?或者你只是觉得他们很烦人?应该放弃闰秒吗?与UTC不同的是，GPS时间系统没有插入闰秒，因此它是一个“统一”的时间系统，但是在UTC中每插入一个闰秒，它就与UTC偏离1。你的工作是否使用GPS时间?让我们知道你的MATLAB时间体验(悲哀?) 在这里 ．
|
类别:
数据类型,
时间
  • 打印
  • 发送电子邮件

另请参阅

评论

要留下评论，请点击在这里登录到您的MathWorks帐户或创建一个新帐户。