1，引言

“天人合一”，是东方文化中最重要的思想，也是中医理论的基础。在这里，“天”不仅包括可以观察到、触摸到的世界，也包括我们目前尚不能直接观察到的世界，如图1所示。如果把这个“一”比作一棵大树，那么 “天”就是这棵大树的树干，世界上每个人都是这棵大树上的一片树叶。我也不例外，也是这棵大树上的一片树叶。

图1 “天人合一”中“天”的图解

“天人合一”思想的具体应用，最基本的一条就是“天”与“人”具有统一的律动、谐振。中国的整个易经易学文化体系都在阐述这个问题。其中每个人的个体特征则表现在其律动幅度与相位（时延）的不同上。用“大树的树叶”来形容这种状况是很恰当的。基于“一叶知秋”的道理，我从1998年2月开始，对自己的体温进行了长期、频繁、规范、精密的测量、记录，每日不辍，至今已超过22个年头，希望藉此叩开“中医科技化”大门。

最近中国武汉爆发的“新型冠状病毒肺炎”疫情，距2003年“非典”爆发时间跨度是17年。我惊奇地发现，我在我的体温数据处理过程中，也呈现出一个“17年周期”。是巧合？还是必然？我现在将我的数据处理过程，及关键图形放在下面，也请各方面的专家、学者研究、思考。

2，数据来源与研究方法

2.1，数据来源及其预处理

图2-1是所测得的体温原始数据趋势图。

数据截止时间段为1998年2月20日～2019年11月9日，共计7933天。采样间隔2小时（1时辰），每天12个采样点，故数据总长度为7933*12=95196点。

测量工具为市场所购“电子体温计”，小数点后面保留两位数字的那种。数据记录时略去十位数“3”及小数点，例如“36.73°C”即记为“673”。测量部位为舌下。

图2.1 体温原始数据趋势图

下面是数据处理时的matlab程序代码：

--------------------以下代码--------------------

x=TW793312(:)';% %确保x为行向量

figure(5);plot(x)

title('原始体温body temperature数据TW793312趋势图)')

%下面利用“经验模态分解”的方法去掉原始数据中的趋势项，

%并使数据均值为0。

imf=emd(x);% 经验模态分解

s=size(imf,1);

x0=sum(imf(1:s-1,:));%去掉最后一项“残余项”

--------------------以上代码--------------------

2.2，求体温数据自相关函数：

--------------------以下代码--------------------

Rx0=xcorr(x0,'unbiased');%求自相关函数

%下面在自相关函数上加个两端收敛的三角窗函数ww，

%以消除数据两端过度发散现象

Lx=length(x);

LRx=length(Rx0);%=2*Lx-1

ww=ones(1,LRx);%单位矩形窗的长度等于LRx

%window is equal to LRx

n=100;%两端各取100个采样点

nn=2*n;

w=triang(nn);%三角窗

ww(1:n)=w(1:n); %窗函数左边

ww(LRx-n+1:end)=w(n+1:end);%窗函数右边

Rx=Rx0.*ww;%加窗

figure(20);plot(Rx)

title('TW793312自相关Autocorrelation函数(加窗)')

--------------------以上代码--------------------

图2.2.1是体温数据的自相关函数图。

可以看到，图形上下包络线由很多“锯齿”构成。由于从图形中心线到两边端头，时间跨度（时延）约为22年，因此很容易得知，每个“锯齿”的时间跨度（时延）为1年。

同时很容易看出，从图形的上包络线的正中间，到两边“大波”的最高处，时间跨度约为17个“锯齿”，也就是约为17年。

将图2.2.1横轴放大，看看图形的短时间内结构。（见图2.2.2）：

图2.2.1 体温数据的自相关函数

图2.2.2 自相关函数图横轴放大

图2.2.3 自相关函数图横轴二次放大

继续放大横轴（见图2.2.3）：图中，一个周期对应1天时间。

2,3，提取自相关函数包络线

前面说了，自相关函数上包络线的“大波”周期约为17年，但这只是根据极粗糙的目测。要知道其准确数值，必须经过一套更严谨的算法。常用的方法有“最小二乘曲线拟合法”。

为了采用这一方法，先将自相关函数上下包络线提取出来。从上面“横轴二次放大”图中可知，其上下包络线分别是由每天体温周期当中的最高值与最低值组成。

--------------------以下代码--------------------

t=1:LRx;

[env,envmoy] = cenvelope(t,Rx,1,'spline');%提取函数图形上下包络线

envb=env(1,:);%下包络线

enva=env(2,:);%上包络线

figure(30)

plot(enva,'r')%上包络线

hold on

plot(envb,'b')%下包络线

title('体温数据自相关函数上、下包络线Envelope')

--------------------以上代码--------------------

图2.3：很容易看出，上包络线是由两条幅值、频率各异的余弦函数叠加合成。

下包络线好像稍微复杂些。由于我们此文的主要目的是提取上包络线的“17年周期”的准确数值，故下包络线暂时不研究。

图2.3 自相关函数图上下包络线

2.4，对上包络线进行最小二乘拟合

单独绘制上包络线图：

--------------------以下代码--------------------

figure(35)

plot(enva,'r')% 上包络线

title('体温数据自相关函数上Upper包络线')

--------------------以上代码--------------------

构造拟合上包络线的原型函数（函数m文件）：

--------------------以下代码--------------------

function csa = Ptt20200123( a0,tt )

%为最小二乘法拟合构造原函数m文件

A10=a0(1);%280;

T1a0=a0(2);%75830;

A20=a0(3);%50;

T2a0=a0(4);%4390;

Ha0=a0(5);%480

cs1a=A10*cos(tt*2*pi/T1a0);

cs2a=A20*cos(tt*2*pi/T2a0);

ha=Ha0*ones(1,length(tt));

csa=cs1a+cs2a+ha;

--------------------以上代码--------------------

在命令窗口输入并运行：

--------------------以下代码--------------------

tt=-Lx+1:Lx-1;

--------------------以上代码--------------------

不断修改“a0=[A10,T1a0,A20,T2a0,Ha0]”的各项系数“A10”,“T1a0”,“A20”,“T2a0”,“Ha0”，关闭前一轮的figure(35)，并运行下述指令：

--------------------以下代码--------------------

figure(35)

plot(enva,'r')%上包络线

A10=280;

T1a0=75830;

A20=50;

T2a0=4390;

Ha0=480;

a0=[A10,T1a0,A20,T2a0,Ha0];

hold on

csa = Ptt20200123( a0,tt )

plot(csa)

title('体温数据自相关函数上Upper包络线')

--------------------以上代码--------------------

目测csa与enva的拟合程度见下图：觉得拟合得很好了，停止，得到一组初始值“a0”（向量）。

在命令窗口调用最小二乘函数，输入变量“a0”,“tt”,“enva”，运行，得到“a”的一组值。

--------------------以下代码--------------------

a=lsqcurvefit(@Ptt20200123,a0,tt,enva);

--------------------以上代码--------------------

.

它就是最终拟合得到的系数（向量）。绘制与“a”对应的csa图形：

--------------------以下代码--------------------

plot(Ptt20200123( a,tt ),'--k')

--------------------以上代码--------------------

见图2.4中黑色虚线：

图2.4 最小二乘法拟合上包络线

2.5，“17年周期”的呈现

看看拟合的结果----工作空间里系数向量“a”的具体数值，见图2.5.1：

将上面两个周期改写成日常习惯的“年”与“天”，见图2.5.2：

看见没有！T1a=16.9798年，几乎是严格等于“17年”。而T2a=365.1918天，也几乎是严格等于“1年（365.2422天）”。

图2.5.1 在工作空间里查看“a”值

图2.5.2 两个周期时间单位换算

3，结果

2003年的“非典”，与现在正在爆发的“新冠”肺炎疫情，两次极为相似的重大事件，相隔17年。这是一个事实。从1998年2月，到2019年11月，在这近22年中，我测量了自己的高精度体温数据，从中提取出了“17年周期”。它也是一个事实。两个事实，两个“17年”，正好相等。你相信这是偶然的巧合，两个事实之间毫无关系吗？我不相信。

略具信号处理知识的人都知道，周期信号的自相关函数周期，与原信号的周期是完全一致的。上述“17年周期”虽然是从自相关函数图的上包络线中提取的，但从图2.2.3可知，上、下包络线其实就是由“1天”周期分量的极大值、极小值组成的。它们的周期信息无疑反映了原信号的周期信息。

在原信号自相关函数上包络线中，有3个大能量周期分量。其中第2大能量的“17年周期”分量，它的周期取值域是时间轴上除“1天”、“1年”两点外的整个正实数。文中数据处理结果“16.9798年”，距“17年整”，相对误差只有0.12%，绝对误差只差了7天零9个小时，从概率统计意义上讲属于“极小概率事件”。“极小概率事件”在一次试验中就发生了，则可认为是“必然事件”。两个“17年”是同一件事情在不同方面的反映。

4，结论分析

4.1，“17年周期”说明此次疫情爆发是一个必然事件。

我就是“天”这棵大树上的一片树叶，我身体中的“17年周期”，也是“天”这棵大树的的周期。这里使用“周期”一词仅就此次数据处理过程中呈现的“余弦波”形状而言，并不代表说数据无限延长后，此余弦波必定仍然会无限延长。从一个余弦波周期“发育”完成的必然性，说明了从2003年“非典”疫情到2020年“新冠”疫情爆发的必然性。

4.2，此次疫情爆发是“天道”逐步“量变”的结果。

同时说明，此次疫情爆发，并不是一个突发的事件。它像我身体中的（高精度）体温一样，是经过了日积月累的、一步一步的“量变”过程，最后才达到疫情爆发的“质变”临界点。由于这种特性，人们要对下次疫情爆发的时间点进行预测，在理论上便具备了可行性。从由密集的数据点形成的波形曲线研究它的未来走向，比由零散数据点预测未来趋势要容易得多。

此文不是专门研究疫情爆发时间规律的，不过从1918年疫情大流行距2003年是5个17年，距2020年是6个17年的历史来看，“17年”的确是研究疫情爆发时间规律的一个很值得重视的时间间隔。

4.3，“天道”节律可以用科技方法获取或验证。

我们人体中的节律，也就是“天（天道）”的节律，除了像“一天24小时（地球自转）”、“一年365.2422天（地球绕日公转）”这样的可以从天文观察中得知的周期，还存在着无法从天文观察中得知的周期。此次研究结果表明，中国人自古以来关于“天”的概念是对的，是可以用现代科技方法验证的。本次数据处理中呈现的“17年周期”可以说是“意料之外，情理之中”，我相信以后从既有经典理论中必可找到其具体、针对性的解释。

5，展望

我对自己身体的测量表明，对个人生理参数长期、高精度的测量，所得数据中含有大量极其重要的（个体及全体）医学生理学信息，绝不都是杂乱无章的噪声。我坚定地相信，如果参加测量的人数足够多，时间足够长，数据精度足够高，最后将所得数据进行处理，结合《黄帝内经》等医易经典所述规则，结合医、易、道学界高人的经验对数据结果意义的解读，可以建立人体系统的数学模型。这既是完成“中医科技化”的必经之途，也是促使已经陷入发展困境中的“现代（循证）医学”产生医学革命的必然之举。那时候人们对疾病的预测、预防、诊断、治疗将变得极为简单而准确。

（本文完）