用质周期分析法
编算富阳市气象历
实验报告
陕西省气象局退休高工张时钊，20世纪60—70年代在富阳市务农期间，于1977年提出“质周期分析法”，并用手工编算气象历（超过一年的逐日天气预报）。该方法曾于1987年在太原举行的预报方法研讨会上介绍，但除了79年的第三期《富阳科技》上有一篇《用质周期分析法作地震周期性预报的尝试》（署名胥口公社中心学校 张时钊）外，从未在报刊书籍上公开发表。2002年与富阳市气象局施国富、王华卿、尹跃虎合作，用富阳市气象局的气象历史资料进行计算机试验，获得许多新结果。预报准确性也显著提高，今报告如下。
因为“质周期分析法”是未公开发表的数学方法，原理虽不复杂，但数学论证较难看懂。所以本文先对它作一简单明了的介绍，再报告实验结果。
一 质周期分析法简介。
地球大气层，或者天气系统的运动变化，产生各种天气现象和性质的演变。我们定时记录下来，就是时间序列，例如逐日气温序列，逐日雨量序列等。除了气象上的各种时间序列外，无数其他事物，尤其是人类尚未完全掌握其规律的特大系统，都产生各种等待我们去分析研究的时间序列。显然，这种时间序列如果有某种规律性，它必然是周期性。因为客观存在是有规律地运动的。即使是人类难以去影响、不能控制的特大系统也不列外。虽然不能控制，或者在控制之前，人类应该去研究它，认识它，首先做到能够预测它。我们要研究的、可以认识的，就是时间序列的周期性，做预报只能靠这种周期性。
已出现的周期分析方法有千百种，我认为最科学可靠的是谐波分析。谐波分析依据数学定理──福里哀变换，即：任何周期函数都可展开成为一系列正弦函数之和。如果是离散的时间序列，则可展开成为三角函数级数。但是谐波分析用在天气预报上不一定成功，更不比其它方法有效。这里显然有两条不正确、不合理的经验：用作分析的时间序列不能太长，所选取的周期数目不能太多。时间序列愈长，数据资料愈多，应该能提供更多的信息，愈能揭示该事物的性质，为什么反而要短些？按数学定理，时间序列应展成无穷的三角级数，为什么只能取少数几个周期？
因为，气象系统非常复杂，任一种时间序列都将包含数不清的周期，短的可能仅几个小时，长的可以到几年，几万年，所以根本不知道它的总周期有多长。而福里哀变换只适用总周期长度已知的函数，它可以把该函数展开为无穷个倍频的正玄函数之和，对总周期未知的序列则无能为力。
张时钊于1977年提出的“质周期分析”，不仅计算非常简单，还可克服以上两个缺点，而实质上又是根据或等价于谐波分析。可以验证，进行所谓的分相统计得到的周期因子，实际上是以该周期为基频的所有正弦函数之和。另外，如果我们要得到该整数周期的正弦函数，只要对这个周期因子进行多次滑动平均就可以了。其结果与通常的谐波分析是一样的，常常就用这种方法做一般的谐波分析。只有周期非整数时（如在我称之为“精密谐波分析”中所要求的），才需要复杂的解方程运算。
现在举列说明什么叫“分相统计”和“质周期因子”。如果有一个只有18个数据的时间序列：F={3，15，21，6，6，27，18，-33，15，0，-6，-9，6，18，9，3，0，9}，
所谓按周期T=6的“分相统计”是：
3 15 21 6 6 27
18 -33 15 0
-6 -9
相加
6 18 9 3
0 9
──────────────────
求和
27 0 45 9
0 27
平均
9 0 15 3
0 9 此为统计因子S6（F）
这个长度为6的统计因子S6（F）={9，0，15，3，0，9}，它含有周期T=2，周期T=3的因子及平均值（即周期T=1的因子），现在把它们都从统计因子S6（F）中减去，就成为质周期因子S6ο（F）了。
注意；在进行这一步“本质化”运算时一定不要有重复。例如在本列中，如果对S6（F）分别进行T=2和T=3的分相统计得到S2[S6（F）]和S3[S6（F）]，再在S6（F）中减去这两者，那么就减了两次平均值S1[S6（F）]，因为S2和S3中都含有平均值S1。正确的本质化运算可以有不同的途径。本例可取以下的标准过程（有时标准的并不是最优的）：
第一步，先求平均值S1=（9+0+15+3+0+9）/6=6，原统计因子减去这一平均值成为：{3，-6，9，-3，-6，3}
第二步，对第一步结果作T=2的分相统计：
3 -6
9 -3
相加 -6 3
────────────
求和 6 -6
平均 2 -2
这是S2
得S2={2，-2}，从第一步结果中减去，可得到：
第一步结果： 3
-6 9 -3 -6
3
S2 ：
2 -2 2 -2
2 -2
──────────────────────
别相减得差： -1
-4 7 1 -8
5 此为第二步结果：
第三步：对第二步结果作T=3的分相统计：
-1 -4 7
1 -8 5
──────────────
求和： 0 -12
12
平均： 0 -6 6 这是S3
得S3={0，-6，6}，从第二步中减去它，最后即可得到质因子S6：
第二步结果： -1 -4
7 1 -8 5
S3 ：
0 -6 6 0
-6 6
──────────────────────
相减： -1 2 1 1 -2 -1 这就是质因子S6
“统计因子”经过“本质化”就成为“质因子”。当统计因子的周期T是质数时，本质化运算只要减去平均值S1，否则必须除去它所含的所有周期等于真因子的质周期函数。反过来，任何统计因子，或这总周期T已知的时间序列，必等于它的所有质周期因子之和：这里说的质周期因子的周期Ti，是整数T的因子，包括T本身及Ti=1。你可以用实列来验证这个结论，而且不难发现，已知周期的序列不等于普通统计因子之和。
那么，对于像气象数据这样不知道总周期长度的时间序列呢？这里可提出下面的规律可供大家验证。如果序列的长度为L，则只要把周期Ti=2到Ti=2L的平方根的所有质周期因子都求出来，然后把它们都加起来，就是实际序列的良好模拟。如果该实际序列中起作用的周期都在2到根号2L之间，则就可用这些周期因子之和作预报。如果起作用的周期有在这个范围之外的，那么就要找出所有这些周期的质因子。要取多少质周期因子？大概是达到使这些质因子的周期长度之和等L为止。
二 编算富阳市气象历的试验结果。
这次在计算机上的试验,即用富阳市1951年到2000年的16772天逐日气象记录,试用质周期分析法编算气象历。因为最后400天要留作预报验证用，实际参加分析的只有16372天数据。二个多月来的计算取得的新进展，主要有两个:
⑴、最困难的日雨量预报有了一些提高。与气温的波动不同，雨量是突发性的。其序列很难看作是一种振动波。今春就开始用新软件研究，第一次用图示显示2万多日雨量的模拟计算结果，竟同实际拟合的很好（当时是用西安市的资料）。但一进入真正的预报期就完全不同了。过去（在富阳用富阳站3-5年的资料），曾对日雨量进行取对数，取平方根等处理，降低其突发性。但这次证明没有多大作用。只有对其平滑化序列，做预报才好些。近两个月来，提出了一些选取周期的新方法和预报评分的新方法，做了大量试验，才明确一些优化参数的方向，得到稍好一些的结果。后面列出对日均温和日雨量的计算结果。
值得一提的是，对于日雨量序列，使用的资料愈长，预报结果愈好。我是用几十种不同的参数组合试验，采用新的预报指标平均值作比较，才得出此结论的。虽然指标值的变化都不大。无论用评分及标准偏差这两种指标的哪一种，使用资料长度每减少5000日，预报准确性都降低。在温度序列的试验中还看不出这样的结果，而对随机数序列的试验，则表明数据太多没有好处。
(2)、因为由机算机产生的随机数是伪随机数，按同样方式产生的随机数序列，看起来是相同的。为此，我是用这样的方法取得供试验用的随机数序列：取得一个随机数后，不仅下一个随机数不要，还按这个随机数的大小，再丢弃一定数目的连续产生的随机数。对这样的随机数序列，进行质周期分析，并试作预报，竟也有一定的准确性：因为计算值序列与原序列有相关性。按统计规律，由两个不相关的随机序列的对应元素相减，得到的差序列的标准差应是原序列标准差的sqr(2)=1.4倍，但我们对随机序列的预报误差的标准差却是1.2倍，比预期的小。
附图1:日均温.红色为计算值,上部为1999.11.26日前400天的模拟值,下部为后400天预报值,黑色为实际值.
附图1:日雨量.红色为计算值,上部为1999.11.26日前400天的模拟值,下部为后400天预报值,黑色为实际值.
注:用富阳市1955-2000共46年16772天资料,最后400天留作预报验证.