DCCA去趋势互相关分析算法详解与Python实现

简介：面向时间序列分析研究者和学生，提供去趋势互相关分析（DCCA）算法的MATLAB实现，用于量化两组非平稳信号间的长期幂律相关性，可广泛应用于气象、金融、生理信号等领域。压缩包共4个文件，包含3个.m源码文件和一个readme.txt说明文档，核心脚本分别覆盖DCCA主函数、综合示例与独立测试用例，txt文件辅助讲解调用方式与参数设置，整体压缩后仅2KB，轻量易用。目前已有995人学习下载，适合需要快速上手互相关分析、开展实验对比或复现论文结果的算法学习者与科研人员。通过阅读源码可清晰理解去趋势协方差计算、波动函数拟合等关键步骤，代码结构简洁，便于在此基础上扩展为多尺度或滑动窗口版本，能有效支持后续研究与实践应用。

1. 去趋势互相关分析的 DCCA 算法到底解决什么问题

两个时间序列在短窗口里看不出同步，在长窗口里却稳定同涨同跌，这种情况在金融、气象和生理信号里都很多见。常见做法是先算 Pearson 相关系数，但 Pearson 对非平稳序列很敏感：只要两个序列各自带趋势，哪怕噪声互不相关，相关系数也会被趋势带偏。DCCA 算法的处理方式是先对每条序列做去均值累加，再按不同尺度切段，分段去趋势后计算残差乘积的平均大小。它把“相关”和“趋势”分开，适合回答两个非平稳序列在特定时间尺度上是否存在交叉相关性。如果你需要可复现的 Python 示例代码，下面这套从核心实现到显著性检验的路径可以直接拿去改。

2. DCCA 算法的去趋势步骤与尺度参数选择

2.1 去趋势互相关分析为什么先要累加

DCCA 的起点是 DFA（去趋势波动分析）。DFA 面对单个序列时，先把原始序列减去均值后做累积和，得到 profile。这个累加操作会把局部波动累积成一条随机游走般的曲线，局部趋势在曲线里变成明显的斜率或弯曲，之后就能用分段多项式去掉。DCCA 把同一套操作同时施加到两个序列上：两个序列各自得到 profile，再按同一套切段方式分段，分别拟合趋势，最后计算残差乘积。

直接对原始序列做 Pearson 相关的问题在于，趋势和噪声混在一起，相关系数往往只反映趋势方向是否一致。DCCA 去趋势后，残差表示的是“去掉局部趋势后还剩下的波动”，残差乘积能反映两个序列在波动层面是否协同。常见做法是把 DCCA 和 DFA 对比着看：如果 DCCA 的波动函数显著大于两个序列各自 DFA 波动函数的乘积期望，就说明除了自身长记忆之外，序列之间存在额外的交叉相关。

2.2 DCCA 的完整算法流程与窗口内计算

给定长度 N 的序列 x 和 y，先算累计和序列 X_k = sum(x_i - mean(x))，Y_k = sum(y_i - mean(y))。对每个尺度 s，把 X 和 Y 切成不重叠的窗口，每段长度 s。每个窗口内用最小二乘多项式分别拟合 X 和 Y 的局部趋势，阶数为 order，然后计算窗口内每个点的残差乘积均值，得到该窗口的协方差贡献 f2(s, v)。把所有窗口的 f2 取平均再开根号，就得到该尺度下的 DCCA 波动函数 F_xy(s)。

具体算法流程按下面几步走：

1. 输入 x、y、scales 列表、order 阶数；
2. 分别对 x、y 去均值并累加，得到 profile X、Y；
3. 对每个尺度 s，从序列开头切出 floor(N / s) 段；
4. 如果 N 除以 s 有余数，再从序列尾部反向切出同样数量的段，把尾部样本用上；
5. 每段内用 np.polyfit 拟合趋势，顺序为 order，残差乘积均值存入列表；
6. 所有段的 f2 取平均后开根号，得到 F_xy(s)；
7. 把 F_xy(s) 对 s 做双对数回归，斜率就是交叉标度指数。

第 4 步是 DCCA 和普通滑动窗口相关最大的区别。反向分段不是为了增加精度，而是为了不浪费尾部余数。如果 N 能被 s 整除，反向段会和正向段重复，可以只算正向段，不影响均值。

2.3 尺度 s、去趋势阶数 order 和分段次数的约定

DCCA 对参数的要求比 Pearson 苛刻，常见问题都出在参数选择上。scales 通常用对数等距生成，因为最终分析依赖 log(s) 对 log(F_xy(s)) 的线性关系。s 太小时，每段内只有几个点，多项式拟合不稳定；s 太大时，有效段数太少，F_xy(s) 的方差变大，回归斜率也失去意义。

下面是项目里常用的默认参数表，可以直接当作起点：

去趋势阶数 order 是容易理解错的点。order=1 表示用直线拟合窗口内的 profile，order=2 表示用二次曲线。金融价格序列的累积和往往接近局部线性，order=1 常用；心电、气候数据里常有缓慢的二次漂移，order=2 会更稳。如果 order 低于真实趋势的阶数，残差里会残留趋势成分，DCCA 会出现虚假的高相关；如果 order 太高，会把有用的低频相关也当成趋势去掉，相关系数偏保守。一般做法是把 order=1 和 order=2 各跑一遍，结果差异很大时说明趋势阶数没选对。

2.4 参数之间的联动关系

scales 和段数不是独立的。尺度 s 越大，段数越少。即使反向分段让段数翻倍，当 s 超过 N/4 时，有效段数通常只有个位数，log-log 回归不可靠。建议先按段数过滤尺度：只保留满足 2 * floor(N / s) >= 10 的 s，再参与回归。另一个联动点是 s 要和 order 匹配。s=8 时用 order=2 拟合一条二次曲线，自由度过低，拟合会完全穿过残差点；所以小尺度下应把 order 控制在 1。算法实现里可以加一个过滤条件，比如 s <= order + 1 时直接跳过，避免 polyfit 报错或返回无意义结果。

3. 用 Python 实现 DCCA：示例代码与参数说明

3.1 最小可运行的 numpy 实现

下面这份代码只依赖 numpy，核心函数很短，适合作为 DCCA 算法的基础版本。代码里把 DFA 和 DCCA 合在同一个函数里，因为当 x 和 y 是同一个序列时，DCCA 就退化为 DFA。

import numpy as np def _profile(x): # 去均值后累加，得到 DFA/DCCA 需要的 profile return np.cumsum(x - np.mean(x)) def _fit(seg, order): # 对窗口内的 profile 做多项式拟合，返回拟合后的趋势值 t = np.arange(len(seg)) coeffs = np.polyfit(t, seg, order) return np.polyval(coeffs, t) def dcca(x, y, scales, order=1): """每个尺度 s 对应的 DCCA 波动函数 F_xy(s)。""" n = len(x) prof_x = _profile(x) prof_y = _profile(y) result = [] for raw_s in scales: s = int(raw_s) if s <= order + 1: # 点太少，拟合没有自由度 continue n_seg = n // s if n_seg < 1: continue f2_list = [] # 正向分段 for v in range(n_seg): sl = slice(v * s, (v + 1) * s) trend_x = _fit(prof_x[sl], order) trend_y = _fit(prof_y[sl], order) f2_list.append(np.mean((prof_x[sl] - trend_x) * (prof_y[sl] - trend_y))) # 反向分段，利用末尾余数样本 for v in range(n_seg): sl = slice(n - (v + 1) * s, n - v * s) trend_x = _fit(prof_x[sl], order) trend_y = _fit(prof_y[sl], order) f2_list.append(np.mean((prof_x[sl] - trend_x) * (prof_y[sl] - trend_y))) result.append(np.sqrt(np.mean(f2_list))) return np.array(result) def dcca_coef(x, y, s, order=1): """同一尺度下去趋势互相关标准化系数，取值范围接近 [-1, 1]。""" s = int(s) fxy = dcca(x, y, [s], order)[0] fxx = dcca(x, x, [s], order)[0] fyy = dcca(y, y, [s], order)[0] return fxy / np.sqrt(fxx * fyy)

3.2 代码分块说明与参数含义

_profile里的去均值非常关键。不去均值时，累积和会带一个随序列长度增加的常数偏移，虽然窗口内拟合的截距能吸收这个偏移，但残差乘积的数值稳定性会变差。去均值之后，profile 起始点归零，不同序列之间的比较更公平。

dcca内部对每个尺度单独计算。f2_list 里每一项都是一个窗口的残差乘积均值，最后对窗口列表取平均再开根号。这里不是先算所有窗口的总和再除以 s，而是每个窗口内部先除以 s，避免不同窗口因为段首段尾位置不同产生权重偏差。

dcca_coef使用同一尺度下的 F_xy(s)、F_xx(s)、F_yy(s) 做标准化。F_xx(s) 就是 DFA 的波动函数，所以代码里直接调用dcca(x, x, [s], order)。这个标准化系数相当于“去掉各自趋势后的相关系数”，数值范围一般落在 -1 到 1 之间，但小样本下偶尔会略超出，不要惊慌。

容易犯的错是把 scales 设成整数列表后忘记转 int。numpy 生成的 logspace 是 float，如果直接传给range会报错。代码里用int(raw_s)做了转换，但如果你自己写循环，也要注意这一点。

3.3 用合成数据验证 DCCA 实现是否正确

验证代码最直接的方法是跑两个 sanity check：第一个是相同序列对相同序列，DCCA 系数必须等于 1；第二个是独立白噪声对独立白噪声，DCCA 系数应该在 0 附近。另外，白噪声对自身做 DCCA 等价于 DFA，log-log 斜率应该接近 0.5。

rng = np.random.default_rng(42) n = 800 a = rng.standard_normal(n) b = rng.standard_normal(n) # 同一序列的 DCCA 系数必然为 1 print("same series rho:", dcca_coef(a, a, 64)) # 独立白噪声的 DCCA 系数应该在 0 附近 print("independent rho:", dcca_coef(a, b, 64)) # 白噪声的 DFA 斜率约 0.5，说明尺度计算没写错 scales = np.unique( np.logspace(np.log10(10), np.log10(n // 4), 15).astype(int) ) f_aa = dcca(a, a, scales) slope = np.polyfit(np.log(scales), np.log(f_aa), 1)[0] print("white noise DFA slope:", round(slope, 3))

独立白噪声的 DCCA 系数不会严格等于 0，因为样本有限，残差乘积均值会有随机波动，落在 -0.2 到 0.2 之间都正常。如果独立序列跑出来接近 0.9，优先检查是不是切段索引写错，或者把同一个序列传进了两个参数。DFA 斜率接近 0.5 可以确认累加、分段、去趋势和平均这四步都没有结构性问题。

4. 批量计算、显著性判断与 DCCA 的常见坑

4.1 多序列批量计算 DCCA 系数矩阵

实际项目里通常要算一个资产池或者一组脑电通道两两之间的 DCCA 相关系数。写一个矩阵函数会很方便，输入形状是 M 行 N 列的数据，输出对称矩阵。

def rho_matrix(data, s, order=1): """data: (M, N)，M 个序列，每个长度 N。""" m = len(data) rho = np.ones((m, m)) for i in range(m): for j in range(i + 1, m): r = dcca_coef(data[i], data[j], s, order) rho[i, j] = rho[j, i] = r return rho

矩阵对角线置为 1，因为序列与自身的 DCCA 系数恒为 1。嵌套循环的时间复杂度是 O(M^2 * N * number_of_segments)，序列数量超过 50 时会明显变慢。如果只需要所有序列在少数几个尺度上的矩阵，可以预先算好每个序列的 DFA 波动函数，再在组合时复用 F_xx，避免 dcca_coef 内部重复计算。

4.2 显著性判断：用相位随机化替代简单置换

DCCA 系数没有现成的置信区间，通常用替代数据法做显著性检验。简单置换测试会把 y 随机打乱，但打乱会破坏 y 自身的自相关结构，导致检验过于乐观。更常用的是相位随机化：对 y 做 FFT，保持幅值谱不变，把相位换成随机值，再逆变换回去。这样生成的替代序列与 y 有相同的功率谱，但交叉相关结构被破坏。

def phase_surrogate(y, rng): """保留 y 的功率谱，打乱相位。""" f = np.fft.rfft(y) phase = rng.uniform(-np.pi, np.pi, size=f.size) return np.fft.irfft(f * np.exp(1j * phase), n=len(y)).real def dcca_pvalue(x, y, s, order=1, n_surr=199, seed=0): """返回 DCCA 系数对应的置换 p 值。""" rng = np.random.default_rng(seed) obs = dcca_coef(x, y, s, order) count = 0 for _ in range(n_surr): y_surr = phase_surrogate(y, rng) r = dcca_coef(x, y_surr, s, order) if abs(r) >= abs(obs): count += 1 return (count + 1) / (n_surr + 1)

p 值用 count + 1 除以 n_surr + 1，是常见的保守估计，因为观测样本本身也算一次。相位随机化不适合极短序列，因为 FFT 在长度很短时无法保留稳定的功率谱。n_surr 至少设 199，如果要发布论文，建议 999 或更高。

4.3 DCCA 实现的常见坑和排查表

下面的排查表来自实际跑数据时最容易踩到的地方，可以直接对照症状找原因。

除了表格里的问题，还有一个隐蔽点：两个序列的长度必须一致，且采样点必须对齐。DCCA 的切段是同时作用于 X 和 Y 的，如果 y 比 x 短，正向分段会错误地把不同时间的点当成同一段。遇到缺失值时，不要直接删掉一侧的数据再拼接，因为拼接点会产生人造趋势。常见做法是先用插值补齐，再做 DCCA，并把有长段缺失的尺度剔除。

5. 用 DCCA 系数矩阵收尾：可视化与参数固化

5.1 把系数矩阵画成热力图

DCCA 系数矩阵结果适合用热力图展示。颜色用红蓝发散色，因为系数可正可负，零值用白色显示最直观。

import matplotlib.pyplot as plt names = ["A", "B", "C", "D"] rho = rho_matrix(data, s=64, order=1) fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(rho, cmap="RdBu_r", vmin=-1, vmax=1) ax.set_xticks(range(len(names))) ax.set_yticks(range(len(names))) ax.set_xticklabels(names) ax.set_yticklabels(names) for i in range(len(names)): for j in range(len(names)): ax.text(j, i, f"{rho[i, j]:.2f}", ha="center", va="center", fontsize=8) plt.colorbar(im, ax=ax) plt.savefig("dcca_rho.png", dpi=200, bbox_inches="tight")

热力图上的每个格子代表一个固定尺度 s 下两个序列的去趋势互相关强度。多尺度对比时，建议对每个尺度单独出一张图，而不是把不同尺度混在一张图里。

5.2 用配置字典固化算法参数

DCCA 参数关联性强，复现时必须记录 scales、order、最少段数、随机种子这几个值。我一般会把参数集中放在一个字典里，每次跑完原样存进文件。

cfg = { "sequence_length": data.shape[1], "scale": 64, "order": 1, "min_segments": 10, "n_surr": 199, "seed": 0, "scales": scales.tolist(), }

这样下次拿到数据时，先检查cfg["scales"]是否仍然满足2 * (N // s) >= min_segments，不满足就先过滤再算。把参数和结果分开保存，后面改 order 做稳健性检验时，能直接知道上一版结果是在什么条件下得到的。

5.3 输出 CSV 时保留尺度字段

批量分析经常需要把 DCCA 结果与其它特征写到一张表里。保存时不要只写相关系数，把尺度列也写进去，避免之后把不同尺度的数字混在一起比较。一个简单的输出方式是：

out = np.column_stack(( scales, dcca(a, a, scales), dcca(a, b, scales), )) np.savetxt( "dcca_output.csv", out, delimiter=",", header="scale,F_aa,F_ab", comments="", )

CSV 第一列是尺度，第二列是序列 a 自身的 DFA 波动函数，第三列是 a 和 b 的 DCCA 波动函数。读取时用np.loadtxt的unpack=True拆分三列即可。header 里一定要写明scale而不是s，避免和序列长度缩写混淆。可视化、参数、CSV 三样东西齐了，DCCA 结果就可以交给下游直接使用。

本文还有配套的精品资源，点击获取