CN114925330A

CN114925330A - 一种不依赖重构参数的自适应跨尺度因果分析方法

Info

Publication number: CN114925330A
Application number: CN202210543472.6A
Authority: CN
Inventors: 邓佳; 王关锁; 蒋暑民; 乔方利
Original assignee: First Institute of Oceanography MNR
Current assignee: First Institute of Oceanography MNR
Priority date: 2022-05-19
Filing date: 2022-05-19
Publication date: 2022-08-19
Anticipated expiration: 2042-05-19
Also published as: CN114925330B

Abstract

本发明涉及一种不依赖重构参数的自适应跨尺度因果分析方法，属于数据处理领域，所述方法包括数据搜集、对搜集的数据进行尺度分解、对关注的模态重构状态空间、在重构的状态空间计算目标统计量和设置置信检验；所述方法在置信检验中提出的判别指标为“有效面积S_ΔMCR(τ,m)”，所述指标是关于重构参数τ和m的函数，是判别两变量是否存在因果关系的最终依据，该指标能够克服重构参数以及替代数据算法偏差对结论的影响，使得结论具有鲁棒性。本发明技术方案具备良好的并行特征，可通过并行计算极大地缩短计算时间。

Description

一种不依赖重构参数的自适应跨尺度因果分析方法

技术领域

本发明属于数据处理领域，具体地涉及一种不依赖重构参数的自适应跨尺度因果分析方法。

背景技术

因果关系即某一事件的发生会导致另一事件的发生，并且前者的改变会导致后者的改变，但反之则不然，比如地形可以决定局地水汽输送的强弱，但反之则不然。因果关系不同于相关关系，它呈现不对称特征，具有方向性。比如地形可以决定局地水汽输送的强弱，但反之不然。两事件之间即使不存在因果关系，但仍然可能表现出虚假的相关关系；相反，它们之间存在因果关系，也不一定表现出相关关系，比如闪电和雷声高度相关，但它们并不存在因果关系；又如在大尺度上，气候变化才是人类历史进程的主要原因，甚至是决定性因素。科学研究的最终目标是阐明事物之间的因果关系，解释所观察到的现象是在哪些因素的影响下产生和形成的。一旦真正理解了因果思维背后的逻辑，就可以在现代计算机上进行模拟，预测事件将来在什么条件下可能发生并得出干预措施中的控制手段，从而为制定对策提供依据。

从观测数据甄别两系统之间的因果关系是一项极具挑战的任务。数据之间因果关系的定义最初由Wiener(Wiener N.The theory of prediction.Modern Mathematics forthe Engineer，1956.)与Granger(Granger C W J.1969.Investigating CausalRelations by Econometric Models and Cross-spectral Methods.Econometrica：J.Econom.Soc，37：424.)先后提出，后被广泛运用于分析经济变量之间的因果关系。因果分析主要分为线性因果分析和非线性因果分析两大类。

(一)线性因果分析方法

Wiener-Granger因果不仅是线性因果分析的基础，也是整个因果分析体系的基石，其核心概念在于“原因”要发生在“结果”之前，也就是说，引入作为“原因”变量的历史数据能够显著减小对作为“结果”变量的预测误差。这种测量因果关系的方法涉及两个关键点：预测与预测误差。Norbert Wiener最早给出了这种因果的数学形式，但其定义中存在的一些困难使其难以应用到数据分析中。之后，Granger采用线性自回归模型进行拟合，以拟合误差作为预测误差，实现了将Wiener-Granger因果应用到数据分析中。但Granger默认变量之间的关系是线性的，其应用也局限在线性因果分析领域。

(二)非线性因果方法

非线性因果分析方法主要包括如下几类：基于信息论的非线性因果分析方法，基于状态空间(phase space，or state space)的非线性因果分析方法，以及基于递归分析的非线性因果分析方法。

大多数非线性因果分析方法以信息论中的“熵”为基础，并由此衍生出很多不同的因果关系检测方法。这些方法的基本思想是一致的：依据一定的标准，将采样的待分析数据粗粒化至不同的网格中，此时，每一个采样数据只能取值于各种情况中某一种。通过统计每种情况出现的频率以及使用频率代替概率，便可得到每种情况出现的概率值，也就获知了变量的统计学特性，进而可以利用熵值的公式进行计算，并最终得到期望的估测结果。基于信息论的因果关系检测方法对大样本问题较为准确，但在样本数量有限且含噪声的情况下，这类方法存在较大偏差。

基于状态空间的非线性因果分析方法需要使用时间-延迟法重构状态空间，再根据“原因”变量流形的临近点是否可用来确定“结果”变量流形的临近点，来判断系统之间是否存在因果关系。这类方法适用于非同步区域，但无法判断广义同步系统的因果关系。此外，重构状态空间面临预先设置重构参数(时间延迟和嵌入维数)的困难，不恰当的参数重构出来的状态空间存在失真的问题。对于随机动力系统，时间延迟参数的选取对重构状态空间至关重要：在实际应用中，比较理想的τ既要保证状态空间中的相点相互独立，但又不能丢失相点之间的相互作用信息；太小的τ会使相邻的相点强相关，从而导致重构的状态空间产生冗余信息；太大的τ又会使相邻相点之间毫无关联，从而导致重构的状态空间丢失信息。常用的选取时间延迟参数τ的方法包括自相关函数法和互信息法。但自相关函数方法并不适用于对非线性系统选取时间延迟参数，互信息法对数据样本要求较为严格，适用于大样本数据。常用的选取嵌入维数的方法为伪临近点法，但该方法对伪临近点的判别存在较强的主观性。

递归或重现(Recurrence)是许多动力系统的基本性质，它由亨利·庞加莱在1890年首次提出。近些年，随着计算机技术的迅猛发展，混沌理论的相关研究也得到了飞速发展，与递归过程相关的一些冗长计算也得以实现。1987年，Eckmann等(Eckmann J.P，Kamphorst S.O，Ruelle D.1987.Recurrence plots of dynamical systems，Europhys.Lett，5：973-977)介绍了递归图(recurrence plot)方法，实现了动力系统递归过程的可视化分析。递归图方法对有限样本、含噪声数据同样具有普适性，因此基于递归条件概率的因果分析方法比基于信息论的因果分析方法更为稳健。但递归图方法同样也需要重构系统的状态空间，因此也面临重构参数选取的困难。

无论是基于信息论的因果分析方法还是基于递归条件概率的因果分析方法，其定义的统计量均为非对称量，在估算时无法避免数值偏差，因此给出置信水平十分必要。替代数据法是分析非线性系统的一种有效方法，由J.Theiler(Theiler J，Eubank S，LongtinA，etal.1992.Testing for nonlinearity in time series：The method of surrogatedata.Physica D，58：77-94.)首次提出用于检测单变量数据是否包含非线性成分，后又被发展用于检验多变量之间的相关性，主要通过传统实现法和受限实现法构造。传统实现法需要从数据中提取完整的模型方程，然后通过这些方程得到替代数据，其困难在于模型方程的阶数和类别难以选定。受限实现法是对原始数据的重排得到的替代数据，把所需的结构强加到随机时间信号上以避免挑选模型方程，可以说是对原始数据的一种“复制”。受限实现法已有的具体算法包括：以傅里叶变换为基础的生成算法，比如FT算法，调节幅度的傅里叶(Amplitude Adjusted Fourier Transform)算法等；伪周期替代数据算法，比如循环相位扰动(cyciic phase permutation)方法使用希尔伯特变换提取数据相位，再对相位进行扰动生成相应的替代数据；孪生替代数据(Twin surrogate)法，该算法通过对数据重构相空间得到系统所有可能的轨迹，再通过这些轨迹和不同的初始条件生成相应的替代数据。该算法适用于混沌系统，但仍然面临选取适合重构参数的困难。此外，任何替代数据算法都存在系统偏差。

发明内容

本发明要解决的技术问题在于提供一种不依赖重构参数的自适应跨尺度因果分析方法，所述方法是在递归图(recurrence plot，RP)方法的基础上提出，适用于小样本、含噪声数据。本发明通过构造具有自适应特征的替代数据，以及创造性地提出了“有效面积”作为判别系统之间因果关系的最终指标，从而克服了重构参数(时间延迟参数和嵌入维数)以及替代数据算法偏差对结论的影响，使得结论具有鲁棒性。

本发明方法可定量评估各种可观测变量之间的因果关系(即两过程是否存在某种非线性耦合作用)与因果方向，即哪个过程是对另一过程的响应，比如通过冰芯资料与地球轨道资料评估末次冰期气候变化与地球轨道参数之间的因果关系与因果方向，通过心率、血压数据与呼吸频率数据研究心血管疾病与呼吸系统之间的因果关系与因果方向等。该方法是从数据分析的角度甄别驱动因子与响应因子，可为建立或改进相关的数学模式(比如地球气候系统模式)提供理论基础，也可为专业人员制定相关方案提供理论依据。

本发明是通过如下技术方案来实现的：

一种不依赖重构参数的自适应跨尺度因果分析方法，所述方法具体步骤如下：

(一)数据搜集

对关注的两个系统

和系统

分别选取可观测变量，并对观测变量进行数据收集工作，记系统

的观测变量为x(t)＝{x(t₁)，x(t₂)，…，x(t_N)}，系统

的观测变量为y(t)＝{y(t₁)，y(t₂)，…，y(t_N)}，其中t＝t_k(k＝1，2，...，N)表示观测时刻，若数据为非均匀时间格点上的数据，也就是Δt_k＝t_k-t_k-1为变量，则将数据插值到均匀时间格点上，使得Δt_k为常数；

(二)对搜集的数据进行尺度分解

对x(t)和y(t)分别做集合经验模态(EEMD)分解，得到两组如公式(1.1)-(1.2)所示的本征模态函数(Intrinsic mode function，简写IMF)，每个IMF都表征了数据中的一种内蕴的固有尺度：

其中，c_j(t)＝{c_j(t₁)，c_j(t₂)，…，c_j(t_N)}和d_j(t)＝{d_j(t₁)，d_j(t₂)，…，d_j(t_N)}表示第j个IMF；r_x(t)＝{r_x(t₁)，r_x(t₂)，…，r_x(t_N)}和r_y(t)＝{r_y(t₁)，r_y(t₂)，…，r_y(t_N)}表示长期趋势项。在公式(1.1)-(1.2)中，x(t)分解得到的模态个数不必等于y(t)分解得到的模态个数；

(三)对关注的模态重构状态空间

通过变量x(t)和变量y(t)的本征模态去重构系统

和系统

的状态空间。对选定的模态c_p(t)(1≤p≤J)和模态

使用时间-延迟法(time-delay embeddingtechnique)重构状态空间；给定时间延迟参数τ和嵌入维数m，其中，τ和m为正整数，选取2≤τ≤40以及2≤m≤40，并在此范围内对每一个τ和m计算目标统计量，得到关于目标统计量的一张曲面；

进一步，所述步骤(三)中c_p(t)重构状态空间的具体步骤：给定时间延迟参数τ和嵌入维数m，其中，τ和m为正整数，按照公式(2)构造向量

由这组向量

构成的集合称之为模态c_p(t)的状态空间，记为

其中向量

称为状态空间

中的相点(phase point)，表征模态c_p(t)的一个可能状态；M表示

中相点的总数；按照同样方法构造模态d_q(t)的状态空间

其中的相点记作

(四)在重构的状态空间计算目标统计量：平均递归条件概率之差(thedifference of mean conditional probabilities of recurrence，ΔMCR)。

基于状态空间

和

中的相点计算递归图(Recurrence Plot，RP)，递归图的定义如公式(3)所示：

在公式(3)中，对c_p(t)和d_q(t)选取相同的重构参数，其各自重构状态空间中的相点数

ε是预先设置的误差容限，取值满足重现率(recurrence rate，RR)不超过1％；

然后，再对模态c_p(t)和模态d_q(t)计算如公式(5)所示的联合递归矩阵

根据

计算如公式(6)所示的递归条件概率(MCR)：

在公式(6)中，

表示条件概率，即在t＝t_i时刻，已知模态c_p(t)处于状态

时，模态d_q(t)处于状态

的概率；根据公式(6)计算得到与模态c_p(t)与模态d_q(t)的因果关系为：

关系式(7)等价地表述如下：定义

ΔMCR＝MCR(d_q|c_p)-MCR(c_p|d_q)， (8)，

于是：

(五)置信检验

在置信检验中提出的判别指标：有效面积S_ΔMCR(τ，m)，所述指标是关于重构参数τ和m的函数，是判别两变量是否存在因果关系的最终依据，置信检验步骤如下：

(5.1)同时使用3种不同的算法对模态c_p(t)和模态d_q(t)分别生成K_s组替代数据(surrogate data)，所述3种算法为改善的时间平移替代数据(improved-time-shiftedsurrogate，ITS)、保持幅度的块相位平移替代数据(preserve-amplitude-block-phase-shifted surrogate，PABPS)，以及调节幅度的傅里叶变换(Amplitude adjusted Fouriertransform，AAFT)替代数据，用

和

(其中k＝1，2，...，K_s)表示原数据的ITS替代数据，

和

表示原数据的PABPS替代数据，

和

表示原数据的AAFT替代数据；

(5.2)给定τ和m，τ＝a，m＝b，a和b为大于2的两个正整数，对K_s组ITS替代数据

和

计算平均递归条件概率之差，记为ΔMCR^ITS，k(a，b)，k＝1，2，...，K_s；通过ΔMCR^ITS，k(a，b)，k＝1，2，...，K_s计算得到它们的0.05-分位数，记为

对K_s组PABPS替代数据

和

计算ΔMCR^PABPS，k(a，b)，k＝1，2，...，K_s，从而得到

对K_s组AAFT替代数据

和

计算ΔMCR^AAFT，k(a，b)，k＝1，2，...，K_s，从而得到

用|ΔMCR_0.05(a，b)|表示|ΔMCR(a，b)|的95％-置信水平，于是|ΔMCR_0.05(a，b)|定义为：

(5.2.2)但若对某些τ和m，比如

存在|ΔMCR(a_i，b_j)|＞|ΔMCR_0.05(a_i，b_j)|，其中i＝1，2，...，N_τ，j＝1，2，...，N_m，存在以下几种情况：若

在(τ，m)坐标系下为离散的点或为线，此时有效面积等于零，即S_ΔMCR(τ，m)＝0，c_p(t)和d_q(t)的因果关系在统计上不显著；若

在(τ，m)坐标系下为连续的面，此时有效面积为大于零的某个正整数，即S_ΔMCR(τ，m)＞0，c_p(t)和d_q(t)存在因果关系，因果方向由ΔMCR(τ，m)在该连续曲面

内的正负符号决定。上述连续均是指自然数意义下的连续性。

进一步，步骤(五)中的ITS生成数据算法包含如下两个步骤：

第一步：对于快速变化的模态，采用时间平移替代数据生成ITS替代数据；

第二步：对于慢变模态，找出其极大和极小值点，将极大和极小值点随机打乱位置，再重新通过三次样条插值等插值方法插值得到ITS替代数据。

进一步，步骤(五)中对时间序列x(t)构造PABPS替代数据，包含如下三个步骤：

第一步：使用希尔伯特变换(Hilbert-transform，HT)提取x(t)的瞬时相位θ(t)和瞬时振幅a(t)；

第二步：保持a(t)不变，对θ(t)按照如下规则生成块相位平移的替代数据(block-phase-shifted surrogaet，用θ^BPS(t)表示)：计算θ(t)的极小值或极大值点个数，记为N_min或N_max个极小或极大值点；将θ^BPS(t)分为大约[M_min/2](或[N_max/2])个块(block)；将这[N_min/2](或[N_max/2])打乱生成θ^BPS(t)替代数据；这里[N_min/2]表示对N_min/2取整数；

第三步：将a(t)和θ^BTS(t)按照公式(10)合成为x(t)的PABPS数据，用x^PABPS(t)表示如下，

x^PABPS(t)＝a(t)cosθ^BPS(t)， (10)。

本发明与现有技术相比的有益效果：

本发明提出了一种不依赖参数的自适应跨尺度因果分析方法(关键创新点包括提出了具有自适应特征的置信水平计算方法和评估因果关系的“有效面积”指标)。本发明属于非线性因果分析方法，适合用于分析复杂的非线性系统之间的因果关系，真实的系统通常都是复杂的非线性系统。现有的非线性因果分析方法主要以信息论为基础，这类方法对数据的样本数量要求较高(样本数量通常以万计)，对有限样本且含噪声数据偏差大(比如地质数据的样本量常常只有几百)。本发明是以对系统的递归分析为基础，适用于样本数量有限且含噪声数据。

递归类因果分析方法需要重构系统的状态空间，本发明采用简单的、常用的时间延迟方法重构系统的状态空间。但时间-延迟方法面临选取合适重构参数(时间延迟和嵌入维数)的困难，不恰当的重构参数会给出失真的状态空间，从而导致截然相反的结论。若数据的样本数量足够大且不含噪声，那么只需嵌入维数大于某一临界值，对任意的时间延迟，通过时间-延迟法重构的状态空间都能够如实地反应系统的拓扑性质，但嵌入维数的这一临界值通常未知。若数据的样本数量有限且含噪声，则时间延迟的选取至关重要，太大或太小的时间延迟都会导致重构的状态空间失真；另外，重构真实状态空间所需的时间延迟与嵌入维数可能满足某种未知的函数关系。时间延迟通常采用自相关函数法和互信息法估计。自相关函数法不适合非线性系统，互信息法不适合有限样本且含噪声数据。嵌入维数通常采用伪邻近点法估计，但该方法在判断临近点时存在主观性。本发明的贡献之一在于对重构参数采用搜索式方案，即让重构参数遍历某一连续的正整数范围。该方案的优势在于无需预先知道重构参数的临界值以及可能存在的函数关系，但该方案会引入不恰当的重构参数。而本发明的关键贡献则在于提出了能够排除错误重构参数的自适应置信水平计算方法和“有效面积”指标。具体而言，本发明对时间延迟在区间[2，40]连续取值，对嵌入维数也在区间[2，40]连续取值；对每一对重构参数(时间延迟和嵌入维数)计算目标统计量；再根据本发明提出的自适应置信水平计算方法和有效面积指标判别目标统计量是否统计显著。

本发明使用替代数据方法计算目标统计量的置信水平。替代数据方法是估计置信水平的一种常用方法，该方法由一个零假设和一套特定的数学算法组成。任何一种替代数据算法都存在系统偏差，单独使用某一种替代数据去研究来自复杂系统的数据，通常存在低估或高估目标统计量的问题，从而导致错误的结论。本发明借鉴集合平均的思想，提出同时使用三套替代数据进行置信检验，以减小不同算法系统误差带来的影响，具体包括：使用已发表的调节幅度的傅里叶变换(AAFT)替代数据检验观测变量的非线性特征，以及使用本发明提出的两套替代数据检验观测变量之间是否存在非线性或线性相互作用关系。在检验两变量是否存在相互作用时，常用的替代数据有孪生(Twin)替代数据、时间平移(time-shifted)和循环相位置换(Cyclic phase permutation，CPP)替代数据。孪生替代数据算法也需要重构状态空间，同样存在选取重构参数的困难，并且孪生替代数据算法还包含其他待定参数。此外，孪生替代数据计算量较大，对存储量要求高，不适用于样本有限的缓变波形数据。时间平移替代数据需要对数据做截断边界的预处理，以避免生成的替代数据产生不连续性。该方法对样本有限的缓变波形数据不适用，会因截断预处理使得数据丢失过多的信息。CPP替代数据可能会过高地估计目标统计量从而得到相反的结论，这是因为当数据出现快波骑在慢波上的现象时，CPP算法对相位的随机打乱不够充分，导致数据中的部分关联性未被打乱并高估目标统计量。在本发明提出的两套替代数据中，其中一套是对时间平移替代数据的改进，主要提出了如何对样本有限的缓变波形数据生成时间平移替代数据；另一套则是对数据相位的极值点进行随机置换，以保证变量间的关联性被充分打乱，克服了CPP替代数据中高估目标统计量的困难。

在理想情况下，对错误重构参数计算的目标统计量，本发明提出的自适应置信水平计算方法能够准确地判别出其统计不显著。但在实际操作中，由于观测数据受到各种噪声干扰(比如观测噪声和系统的动力噪声)、计算过程存在不可避免的误差(主要表现为构造替代数据算法的系统误差)，结果是可能出现偏差的。也就是说，对某些错误重构参数计算的目标统计量，任何一种置信水平计算方法都可能给出其统计显著的错误结论。但本发明提出的自适应置信水平计算方法具有如下特征：尽管对某些错误的重构参数计算的目标统计量，本发明的置信水平方法也无法识别出其统计不显著，但这些错误的参数呈现出或为孤立的点或为独立线条的几何特征。因此，本发明创造性地提出“有效面积”这一指标作为判别系统之间因果关系的最终依据。这是因为，在数学上，点或线的面积测度均为零。事实上，本发明提出的自适应置信水平计算方法和“有效面积”指标都具备相应的理论依据：在自适应置信水平计算中，替代数据的设计利用了系统的相位更易受到外界影响这一事实；“有效面积”指标则是利用了重构的状态空间对重构参数存在某种连续性。

综上所述，本发明提出自适应置信水平计算方法和评估因果关系的“有效面积”指标，使得结论不再依赖重构参数的选取，保证了结论的鲁棒性。本发明技术方案具备良好的并行特征，可通过并行计算极大地缩短计算时间，因此也具备实际的应用潜力。本发明为不依赖重构参数的因果分析提供了一个可行的框架，其关键之处在于对重构参数进行搜索式的计算思想，以及由此对置信水平引申出的“有效面积”这一概念。在本发明的技术方案中，对数据进行尺度分解的集合经验模态分解(EEMD)方法可替换为变分模态分解(VMD)方法等其他分解方法；计算置信水平所使用的替代数据构造算法也可进一步发展，但仍然遵从本发明提出的技术框架。

附图说明

图1为来自两相互独立系统的变量(因此x(t)与y(t)不存在因果关系)的EEMD分解图。

图2为图1中模态c₂(t)和d₅(t)在重构参数(τ，m)坐标系下的ΔMCR(τ，m)曲面(可见ΔMCR(τ，m)曲面上的值可正可负)图。

图5为来自格陵兰岛NGRIP冰芯的氧同位素(δ¹⁸O，单位为‰)数据以及黄赤交角数据图，其单位为度(°)：黄赤交角在110-0千年至今的变化范围为约22°至24°之间。这里为方便对比，我们将黄赤交角减去一个常数值，平移地画在了(时间，δ¹⁸O)坐标系下。

图6为氧同位素(δ¹⁸O，单位为‰)数据的集合经验模态(EEMD)分解以及黄赤交角模态图：由于黄赤交角已经为规则波形，因此我们仅将其减去均值就可得其模态d₁(t)。

图7为模态c₃(t)和模态d₁(t)在67-0千年至今这一时期的替代数据图。

具体实施方式

下面通过实施例来对本发明的技术方案做进一步解释，但本发明的保护范围不受实施例任何形式上的限制。

实施例1

(一)数据搜集

对关注的两个系统

和系统

分别选取恰当的可观测变量，并对观测变量进行数据收集工作。记系统

的观测变量为x(t)＝{x(t₁)，x(t₂)，…，x(t_N)}，系统

的观测变量为y(t)＝{y(t₁)，y(t₂)，…，y(t_N)}，其中t＝t_k(k＝1，2，...，N)表示观测时刻。若数据为非均匀时间格点上的数据，也就是Δt_k＝t_k-t_k-1为变量，则将数据插值到均匀时间格点上，使得Δt_k为常数。

(二)对搜集的数据进行尺度分解

对x(t)和y(t)分别做集合经验模态(EEMD)分解，得到两组如公式(1.1)-(1.2)所示的本征模态函数(Intrinsic mode function，IMF)，每个IMF都表征了数据中的一种内蕴的固有尺度：

c_j(t)＝{c_j(t₁)，c_j(t₂)，…，c_j(t_N)}和d_j(t)＝{d_j(t₁)，d_j(t₂)，…，d_j(t_N)}表示第j个IMF；r_x(t)＝{r_x(t₁)，r_x(t₂)，…，r_x(t_N)}和r_y(t)＝{r_y(t₁)，r_y(t₂)，…，r_y(t_N)}表示长期趋势项，为单调函数或者常数，在分析中不予考虑。此外，在公式(1.1)-(1.2)中，x(t)分解得到的模态个数不必等于y(t)分解得到的模态个数，也就是说，公式(1.1)中的J与公式(1.2)中的

可以不相等。作为示例，对两个相互独立的系统(因此来这两个系统的任何观测变量都不存在因果关系)分别取观测变量x(t)与y(t)(因此，这两个变量不存在因果关系)，图1展示了x(t)和y(t)的EEMD分解。

(三)对关注的模态重构状态空间。

因果关系是在系统的状态空间中计算得到，因此需要重构系统

和系统

的状态空间。研究证明：整个系统状态空间的动力学可以通过测量系统某个变量的一个时间序列进行重构。因此，通过变量x(t)和变量y(t)的本征模态去重构系统

和系统

的状态空间。对选定的模态c_p(t)(1≤p≤J)(比如图1中的c₂(t))和模态

使用时间-延迟法(time-delay embedding technique)重构状态空间。这里以c_p(t)为例阐明其具体步骤：给定时间延迟参数τ和嵌入维数m(这里，τ和m为正整数)，按照公式(2)构造向量

由这组向量

构成的集合称之为模态c_p(t)的状态空间(记为

)；其中向量

称为状态空间

中的相点(phase point)，表征模态c_p(t)的一个可能状态；M表示

中相点的总数。类似地，可构造模态d_q(t)的状态空间

其中的相点记作

前人研究证明：对从某一确定性的非线性动力系统所收集的观测变量，当数据不包含噪声时，只要嵌入维数m大于某一临界值(即m＞2D+1，其中D表示该动力系统的自由度)，那么对任意的τ值，由公式(2)所描述的状态空间与该非线性动力系统真实的状态空间在拓扑上是等价的。这也就是说，可以通过研究由公式(2)重构出来的状态空间去了解系统的拓扑结构。但当观测变量含有噪声或者系统为随机动力系统时，时间延迟参数的选取对重构状态空间至关重要：在实际应用中，太小的τ会使相点之间不可分辨，从而产生大量冗余信息；太大的τ又会使相点之间丢失内蕴的因果关联。常用的选取时间延迟参数τ的方法包括自相关函数法和互信息法。但自相关函数方法并不适用于对非线性系统选取时间延迟参数，互信息法对数据样本要求较为严格，适用于大样本数据[9]。常用的选取嵌入维数的方法为伪临近点法，但该方法对伪临近点的判别存在较强的主观性。在本实施例中，为避免各种方法对选取重构参数的偏差，对τ和m选取合适的范围进行搜索，比如选取2≤τ≤26以及2≤m≤23，在此范围内对每一个τ和m计算目标统计量，从而得到关于目标统计量的一张曲面。显然，在2≤τ≤26和2≤m≤23的范围，仅有部分重构参数是合理的，因此在目标统计量的曲面上，也仅有部分值是合理的。对于不合理的重构参数，其对应的目标统计量应为统计不显著；对于合理的重构参数，其对应的目标统计量可能统计显著(说明系统之间存在因果关系)，也可能统计不显著(说明系统之间不存在因果关系)。为达到上述目标，本发明将引入自适应置信水平计算方法和有效面积指标对目标统计量曲面进行判别。

基于状态空间

和

在公式(3)中，对模态c_p(t)和模态d_q(t)选取的重构参数(τ和m)通常不同，因此其各自重构状态空间中的相点数(M和

)也不同，但此时问题的维度为四维，包含四个变量：c_p(t)的时间延迟和嵌入维数变量，以及d_q(t)的时间延迟和嵌入维数变量。为降低问题的复杂度，本实施例借鉴非线性动力学中庞加莱截面这一降维技术的思想，对c_p(t)和d_q(t)选取相同的重构参数(τ和m)，此时

问题简化为二维，包含c_p(t)和d_q(t)的时间延迟和嵌入维数两个变量。但即使不做简化，原四维问题(即对c_p(t)和d_q(t)选取不同的重构参数)也可在本发明提供的框架下求解。ε是预先设置的误差容限，通常取值满足重现率(recurrence rate，RR)不超过1％，重现率(RR)的定义见公式(4)：

取R_*(t_k，t_l)中右下角的*为模态c_p(t)为例：

描述了系统

中模态c_p(t)的递归属性，若状态

与状态

相似，则矩阵元素

取值为1，系统出现了递归，因此从递归矩阵可以看出系统

中的模态c_p(t)何时出现相似状态。同样可解释递归矩阵

接下来，再对模态c_p(t)和模态d_q(t)计算如公式(5)所示的联合递归矩阵

类似的，从联合递归矩阵能够看出：在何时，系统

的模态c_p(t)和系统

的模态d_q(t)同时出现了(各自的)相似状态：

根据

计算如公式(6)所示的递归条件概率(MCR)：

在公式(6)中，

表示条件概率，即在t＝t_i时刻，已知模态c_p(t)处于状态

时，模态d_q(t)处于状态

的概率。根据公式(6)可计算得到与模态c_p(t)与模态d_q(t)的因果关系为：

关系式(7)可以等价地表述如下：定义

ΔMCR＝MCR(d_q|c_p)-MCR(c_p|d_q)， (8)，

于是：

如步骤(三)所述，本实施例不再使用任何技巧去估计τ和m，而是让τ和m遍历某个范围，比如取τ＝2，…，26以及m＝2，…，23；再对该范围内的每一对(τ，m)取值，比如τ＝2，m＝3(记为(τ，m)＝(2，3))，依照公式(3)-(8)分别计算ΔMCR，于是ΔMCR在(τ，m)坐标系下张成一张曲面，是关于τ和m的函数，即ΔMCR＝ΔMCR(τ，m)。作为示例，在图2展示了图1中c₂(t)和d₅(t)之间的ΔMCR(τ，m)曲面。可以看到，该曲面上的值可正可负，这说明ΔMCR(τ，m)对不同的重构参数可以给出不同的结论。另一方面，并非所有的(τ，m)都是重构模态c_p(t)和d_q(t)状态空间的最优参数；在实际计算中，方法的系统误差、计算误差、机器误差以及噪声都会影响ΔMCR(τ，m)取值，因此需要对ΔMCR(τ，m)曲面做置信检验。在本实施例中，希望置信检验能够自动地排除不恰当重构参数计算的ΔMCR(τ，m)值，具体而言：对不恰当的重构参数，比如τ＝2，m＝3，置信检验能够判别出ΔMCR(2，3)的值在统计上不显著。但这一点在实际计算中很难达到，总是存在某些错误的重构参数，使得置信检验仍然会给出ΔMCR(τ，m)在统计上显著(也就是“两变量存在因果关系”)这一结论。在这里，本实施例创造性地提出了“有效面积”这一指标来解决上述困难，详细介绍参见步骤(五)。

(五)置信检验。

任何从数据计算得到的统计量，都应对其做置信检验以评估其可靠性。在这一步骤，主要描述本发明在置信检验中提出的判别指标：有效面积(记为S_ΔMCR(τ，m))。该指标是关于重构参数τ和m的函数，是判别两变量是否存在因果关系的最终依据。在步骤6)，我们将详细阐述置信检验中所使用的替代数据构造算法。这里，对ΔMCR(τ，m)取绝对值(即|ΔMCR(τ，m)|)，是因为根据公式(9)：ΔMCR(τ，m)可为大于零的正数，也可为小于零的负数，因此ΔMCR(τ，m)在统计上是否显著取决于其绝对值是否显著，而ΔMCR(τ，m)的正负符号仅用于判断因果方向，即是模态c_p(t)驱动d_q(t)，还是模态d_q(t)驱动模态c_p(t)。置信检验步骤可概括为如下三步：

(5.1)使用不同的算法对模态c_p(t)和模态d_q(t)分别生成K_s组替代数据(surrogate data)。在本发明中，我们同时使用了3种算法去生成替代数据，包括：改善的时间平移替代数据(improved-time-shifted surrogate，ITS)、保持幅度的块相位平移替代数据(preserve-amplitude-block-phase-shifted surrogate，PABPS)，以及调节幅度的傅里叶变换(Amplitude adjusted Fourier transform，AAFT)替代数据，其中ITS和PABPS替代数据的构造算法由本发明提出，AAFT替代数据的算法可参考文献(Surrogate data forhypothesis testing of physical systems.Physics Reports.)。记

和

(其中k＝1，2，...，K_s)表示原数据的ITS替代数据，

和

表示原数据的PABPS替代数据，

和

表示原数据的AAFT替代数据。

ITS替代数据的零假设为：两系统不存在非线性和线性关系。PABPS替代数据的零假设为：两系统的相位在动力学上相互独立。这两种替代数据的主要思想是按照一定的规则扰动原始数据的位置。AAFT替代数据的零假设是：系统的非线性是由观测手段导致，而并非系统内在属性。ITS替代数据生成算法是对时间平移替代数据(time-shiftedsurrogate)算法的改进，在ITS中，我们主要针对慢变模态的时间平移(time-shifted)替代数据做了改进。包含如下两个步骤：

第一步：对于快速变化的模态，采用时间平移替代数据生成ITS替代数据，具体算法可参考文献Surrogate data for hypothesis testing of physical systems.PhysicsReports.。

第二步：对于慢变模态，找出其极大和极小值点，将极大和极小值点随机打乱位置，再重新通过三次样条插值等插值方法插值得到ITS替代数据。这一算法是由本发明提出的。

对时间序列x(t)构造PABPS替代数据，包含如下三个步骤：

第一步：使用希尔伯特变换(Hilbert-transform，HT)提取x(t)的瞬时相位(用θ(t)表示)和瞬时振幅(用a(t)表示)；

第二步：保持a(t)不变，对θ(t)按照如下规则生成块相位平移的替代数据(block-phase-shifted surrogaet，用θ^BPS(t)表示)：计算θ(t)的极小值(或极大值)点个数，记为N_min个极小值点(或N_max个极大值点)；将θ^BPS(t)分为大约[N_min/2](或[N_max/2])个块(block)相位；将这[N_min/2](或[N_max/2])个块相位打乱生成θ^BPS(t)替代数据。这里[N_min/2]表示对N_min/2取整数。

第三步：将a(t)和θ^BTS(t)按照公式(10)合成为x(t)的PABPS数据(用x^PABPS(t))，

x^PABPS(t)＝a(t)cosθ^BPS(t)， (10)。

(5.2)给定τ和m，比如τ＝2，m＝3，对K_s组ITS替代数据

和

按照公式(3)-(8)分别计算平均递归条件概率之差(ΔMCR)，这里记为ΔMCR^ITS，k(2，3)，k＝1，2，...，K_s；通过ΔMCR^ITS，k(2，3)，k＝1，2，...，K_s可以计算得到它们的0.05-分位数，记为

类似的，对K_s组PABPS替代数据

和

计算ΔMCR^PABPS，k(2，3)，k＝1，2，...，K_s，从而得到

对K_s组AAFT替代数据

和

计算ΔMCR^AAFT，k(2，3)，k＝1，2，...，K_s，从而得到

用|θMCR_0.05(2，3)|表示|ΔMCR(2，3)|的95％-置信水平，于是|MCR_0.05(2，3)|定义为：

(5.2.2.1)若c_p(t)和d_q(t)存在因果关系，那么当嵌入维数大于某个临界值之后(我们记这个临界值为m₀)，对大于m₀的任意嵌入维数和任意的时间延迟参数，我们都能对c_p(t)和d_q(t)重构出质量较好的状态空间，从而得到它们的因果关系，也就是：|ΔMCR(τ，m)|＞|ΔMCR_0.05(τ，m)|，m＞m₀，τ＞0。但对于含噪声有限长的时间序列，其嵌入维数和时间延迟参数存在一个范围m₀≤m₁≤m≤m₂，τ₁≤τ≤τ₂(这里m₀，m₁，m₂，τ₀，τ₁表示某个正整数，比如m₀＝12等)，上述不等式仅在此范围内成立，也就是，

|ΔMCR(τ，m)|＞|ΔMCR_0.05(τ，m)|，m₁≤n≤m₂，τ₁≤τ≤τ₂，

这意味这|ΔMCR(τ，m)|和|ΔMCR_0.05(τ，m)|的相交曲面在(τ，m)-坐标系下是一块面积大于0的连续曲面(这里的连续是在自然数意义下的连续，而并非实数意义下的连续)，而并非面积测度为零的孤立点或者线。

(5.2.2.2)若c_p(t)和d_q(t)不存在因果关系，那么即使|ΔMCR(τ，m)|曲面和|ΔMCR_0.05(τ，m)|曲面相交，其曲面相交之处也表现为孤立的点或线，因而面积测度为零。

实施例2

本实施例是不依赖重构参数的自适应跨尺度因果分析方法的流程图。

S1步骤、收集数据，并将数据插值到均匀的时间格点上。

具体的，在本实施例中，我们以末次冰期的Dansgaard-Oeschger(DO)事件是否与黄赤交角(obliquity)存在因果关系来说明本发明技术。DO事件是指发生在末次冰期的一系列气候快速变化事件(约24次)，平均周期为约1470年(Grootes P M，MStuiver.1997.Oxygen 18/16 variability in Greenland snow and ice with 10-3-to-105 year time resolution.Journal of Geophysical Research，102(C12)：26455-26470；)。对于DO事件我们收集了一支来自格陵兰岛NGRIP冰芯的氧同位素(δ¹⁸O，单位为‰)数据(North Greenland Ice Core Project members.2004.High-resolution record ofNorthern Hemisphere climate extending into the last interglacialperiod.Nature，431：147-151.)。表征DO事件的δ¹⁸O时间序列和黄赤交角数据(Berger A，Loutre MF.Insolation values for the climate of the last 10 millionyears.Quaternary Sci.Rev.199110：297-317.)(单位为度，即°)可参见图5。在110-0千年至今这一时间段，黄赤交角的变化范围为约22°至24°之间，为方便对比，我们在图5中将黄赤交角减去某个常数值，平移地画到了时间(千年至今)-δ¹⁸O坐标系下。此外，由于δ¹⁸O数据分布在不均匀时间网格点上，因此我们使用分段埃米特插值将其插值到均匀网格时间节点上，时间步长为50年。

S2步骤、在图6中，我们使用集合经验模态分解(EEMD)将δ¹⁸O时间序列分解为如公式(1.1)所示的形式。因黄赤交角为规则的波形，因此我们仅需将其减去均值便可得到其本征模态d₁(t)。

S3步骤、我们关注尺度为约1470年的气候变化模态与黄赤交角的因果关系。计算图6中不同的本征模态平均频率，可知DO事件由c₃(t)(其平均周期为约1300年)表征。因此我们关注模态c₃(t)与d₁(t)的因果关系。

S4步骤、计算c₃(t)和d₁(t)的平均递归条件概率之差(ΔMCR(τ，m))曲面。同时，对c₃(t)和d₁(t)计算95％-置信曲面，即|ΔMCR_0.05(τ，m)|曲面。为得到|ΔMCR_0.05(τ，m)|曲面，首先计算

和

曲面，再对每一对重构参数(τ，m)取这三个置信水平中的最大值作为|ΔMCR_0.05(τ，m)|(参见公式(10))。已有模式方法证明在海洋同位素阶段3时期(Marine isotope stage 3，MIS3，约60-30千年至今)，在轨道变化作用下可形成DO事件。我们指出，古气候模式与现代气候模式的控制因子不完全相同，而这些控制因子在远久的地质时期通常是未知的，因此古气候模式存在较大的不确定性，并且能够模拟的地质时期也是有限的。这里，我们通过观测资料，重点分析DO事件与黄赤交角在67-0千年至今这一时期的因果关系。在图7中，我们对c₃(t)和d₁(t)仅展示了ITS，PABPS和AAFT替代数据中的其中一组。对这100组替代数据，可以计算得到100组|ΔMCR^ns，k(τ，m)|，|ΔMCR^FABPS，k(τ，m)|和|ΔMCR^AAFT，k(τ，m)|，其中k＝1，2，...，100，再分别计算其0.05-分位数从而得到

和

最后再通过公式(10)计算95％-置信曲面，即|ΔMCR_0.05(τ，m)|曲面。

S5步骤、为了方便判别因果方向，在图7中直接给出了模态c₃(t)和d₁(t)的ΔMCR(τ，m)曲面，以及+|ΔMCR_0.05(τ，m)|(正95％-置信曲面)曲面和-|ΔMCR_0.05(τ，m)|(负95％-置信曲面)。可以看到，ΔMCR(τ，m)曲面与+|ΔMCR_0.05(τ，m)|曲面出现了相交，并且有效面积为S_ΔMCR(τ，m)＝49，为大于零的正整数，因此c₃(t)和d₁(t)在67-0千年至今这一时期存在因果关系，因果方向为黄赤交角驱动DO事件(参见公式(9))。图7包含但不限于已发表的古气候模式结果(Zhang X.*，S.Barker，G.Knorr，G.Lohmann，R.Drysdale，Y.B.Sun，D.Hodelland F.H.Chen.2021：Direct astronomical influence on abrupt climatevariability.Nature Geosci 14：819-826.)：古气候模式仅模拟了DO事件在40-32千年至今这一时期与黄赤交角的因果关系，发现黄赤交角对DO事件存在驱动作用；但图8进一步指出：在67-40千年至今以及32-0千年至今，黄赤交角对DO事件也是存在驱动作用的。通过本发明的技术，这是第一次从观测数据中明确黄赤交角对DO事件有驱动作用，因此结论具有客观性。

本发明实施例提供的技术方案，是从观测数据中计算变量之间的因果关系，其结论不依赖对重构参数(τ和m)的选取，具有客观性和鲁棒性。该技术还可应用于大气和海洋观测数据，从观测数据中提取驱动因子和响应因子，从而改进气候系统模式和海洋环流数值模式中的物理过程；可应用于环境监测数据，从观测数据中判别环境事物的因果关系，为制定相关环保政策提供依据；可应用于医疗健康数据，从数据中判别出疾病的致病因素；可应用于经济或社会科学数据，从数据中揭示市场内部和社会的复杂性，为制定相关政策提供依据。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims

1.一种不依赖重构参数的自适应跨尺度因果分析方法，其特征在于所述方法具体步骤如下：

(一)数据搜集

对关注的两个系统

和系统

的观测变量为x(t)＝{x(t₁),x(t₂),…,x(t_N)}，系统

的观测变量为y(t)＝{y(t₁),y(t₂),…,y(t_N)}，其中t＝t_k表示观测时刻，且k＝1,2,...,N，若数据为非均匀时间格点上的数据，也就是Δt_k＝t_k-t_k-1为变量，则将数据插值到均匀时间格点上，使得Δt_k为常数；

(二)对搜集的数据进行尺度分解

对x(t)和y(t)分别做集合经验模态分解，得到两组如公式(1.1)-(1.2)所示的本征模态函数，简写IMF：

其中，c_j(t)＝{c_j(t₁),c_j(t₂),…,c_j(t_N)}和d_j(t)＝{d_j(t₁),d_j(t₂),…,d_j(t_N)}表示第j个IMF；r_x(t)＝{r_x(t₁),r_x(t₂),…,r_x(t_N)}和r_y(t)＝{r_y(t₁),r_y(t₂),…,r_y(t_N)}表示长期趋势项；

(三)对关注的模态重构状态空间

通过变量x(t)和变量y(t)的本征模态去重构系统

和系统

在对应尺度上的状态空间；对选定的模态c_p(t)和模态d_q(t)使用时间-延迟法重构状态空间,其中c_p(t)中1≤p≤J，d_q(t)中

给定时间延迟参数τ和嵌入维数m，其中τ和m为正整数，选取2≤τ≤40以及2≤m≤40，并在此范围内对每一个τ和m计算目标统计量，得到关于目标统计量的一张曲面；

(四)在重构的状态空间计算目标统计量：平均递归条件概率之差ΔMCR

基于状态空间

和

中的相点计算递归图RP，递归图的定义如公式(3)所示：

ε是预先设置的误差容限，取值满足重现率RR不超过1％；

根据

计算如公式(6)所示的递归条件概率(MCR)：

在公式(6)中，

表示条件概率，即在t＝t_i时刻，已知模态c_p(t)处于状态

时，模态d_q(t)处于状态

若MCR(d_q|c_p)＜MCR(c_p|d_q)，则

系统的模态c_p(t)驱动

系统的模态d_q(t)；

若MCR(d_q|c_p)＞MCR(c_p|d_q)，则

系统的模态d_q(t)驱动

系统的模态c_p(t)；(7)关系式(7)等价地表述如下：定义

ΔMCR＝MCR(d_q|c_p)-MCR(c_p|d_q), (8)，

于是：

若ΔMCR＜0，则

系统的模态c_p(t)驱动

系统的模态d_q(t)；

若ΔMCR＞0，则

系统的模态d_q(t)驱动

系统的模态c_p(t)；(9)

(五)置信检验

在置信检验中提出的判别指标：有效面积S_ΔMCR(τ,m)，所述指标是关于重构参数τ和m的函数，是判别两变量是否存在因果关系的最终依据，置信检验步骤如下：

(5.1)同时使用3种不同的算法对模态c_p(t)和模态d_q(t)分别生成K_s组替代数据，所述3种算法为改善的时间平移替代数据，简称ITS、保持幅度的块相位平移替代数据,简称PABPS，以及调节幅度的傅里叶变换替代数据，简称AAFT，用

和

表示原数据的ITS替代数据，

和

表示原数据的PABPS替代数据，

和

表示原数据的AAFT替代数据，其中k＝1,2,...,K_s；

(5.2)给定τ和m，τ＝a,m＝b，a和b为大于2的两个正整数，对K_s组ITS替代数据

和

计算平均递归条件概率之差，记为ΔMCR^ITS,k(a,b),k＝1,2,...,K_s；通过ΔMCR^ITS,k(a,b),k＝1,2,...,K_s计算得到它们的0.05-分位数，记为

对K_s组PABPS替代数据

和

计算ΔMCR^PABPS,k(a,b),k＝1,2,...,K_s，从而得到

对K_s组AAFT替代数据

和

计算ΔMCR^AAFT,k(a,b),k＝1,2,...,K_s，从而得到

(5.2.2)但若对某些τ和m，比如

存在|ΔMCR(a_i,b_j)|＞|ΔMCR_0.05(a_i,b_j)|，其中i＝1,2,...,N_τ,j＝1,2,...,N_m，存在以下几种情况：若

在(τ，m)坐标系下为离散的点或为线，此时有效面积等于零，即S_ΔMCR(τ,m)＝0，c_p(t)和d_q(t)的因果关系在统计上不显著；若

在(τ,m)坐标系下为连续的面，此时有效面积为大于零的某个正整数，即S_ΔMCR(τ,m)＞0，c_p(t)和d_q(t)存在因果关系，因果方向由ΔMCR(τ,m)在该连续曲面

内的正负符号决定，因果关系是指两过程是否存在某种非线性耦合作用，因果方向是哪个过程是对另一过程的响应，甄别驱动因子与响应因子。

2.根据权利要求1所述的方法，其特征在于所述步骤(三)中c_p(t)重构状态空间的具体步骤：给定时间延迟参数τ和嵌入维数m，其中，τ和m为正整数，按照公式(2)构造向量

由这组向量

构成的集合称之为模态c_p(t)的状态空间，记为

其中向量

称为状态空间

中的相点，表征模态c_p(t)的一个可能状态；M表示

中相点的总数；按照同样方法构造模态d_q(t)的状态空间

其相点记作

3.根据权利要求1所述的方法，其特征在于步骤(五)中的ITS生成数据算法包含如下两个步骤：

4.根据权利要求1所述的方法，其特征在于步骤(五)中对时间序列x(t)构造PABPS替代数据，包含如下三个步骤：

第一步：使用希尔伯特变换提取x(t)的瞬时相位θ(t)和瞬时振幅a(t)；

第二步：保持a(t)不变，对θ(t)按照如下规则生成块相位平移的替代数据θ^BPS(t)：计算θ(t)的极小值或极大值点个数，记为N_min或N_max个极小或极大值点；将θ^BPS(t)分为[N_min/2]或[N_max/2]个块；将这[N_min/2]或[N_max/2]打乱生成θ^BPS(t)替代数据；这里[N_min/2]表示对N_min/2取整数；

x^PABPS(t)＝a(t)cosθ^BPS(t)，(10)。