C-C法计算混沌时间序列相空间重构延迟时间

原创于 2026-06-18 10:15:31 发布 · 204 阅读

6 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#重构

一、C-C方法基本原理

C-C法（Correlation Integral based method）由 Kim、Eykholt 和 Salas 于 1999 年提出，利用关联积分同时估计相空间重构的延迟时间 τ 和嵌入维数 m，但更常用于确定延迟时间。其核心思想是：

对于给定的时间序列 ${x_i\}_{i=1}^N$ ，定义关联积分：

$\frac{2}{M(M-1)} \sum_{1\le i<j\le M} \Theta(r - \| \mathbf{X}_i - \mathbf{X}_j \|)$

其中 $Xi=(xi,xi+t,…,xi+(m−1)t)\mathbf{X}_i = (x_i, x_{i+t}, \dots, x_{i+(m-1)t})$ 为重构相点， $M = N - (m - 1) t$ 为相点数， $Θ\Theta$ 为 Heaviside 函数， $r$ 为邻域半径。

然后构造检验统计量：

$S(m, N, r, t) = C(m, N, r, t) - C^m(1, N, r, t)$

实际计算中，将时间序列分成 $t$ 个不相交的子序列，分别计算 $S$ 后平均，记为 $S‾(t)\overline{S}(t)$ 。同时计算 $ΔS‾(t)\Delta\overline{S}(t)$ （ $S$ 随 $r$ 的最大偏差）和 $Scor(t)=ΔS‾(t)+∣S‾(t)∣S_{\text{cor}}(t) = \Delta\overline{S}(t) + |\overline{S}(t)|$ 。

第一个零点 $S‾(t)\overline{S}(t)$ 的第一个零点对应最优延迟时间 $τ\tau$ ；
第一个极小值 $ΔS‾(t)\Delta\overline{S}(t)$ 的第一个局部极小值也对应 $τ\tau$ ；
全局最小值 $Scor(t)S_{\text{cor}}(t)$ 的全局最小值对应嵌入窗宽 $τw=(m−1)τ\tau_w = (m-1)\tau$ ，由此可推算嵌入维数 $m$ 。

实际中，常取 $S‾(t)\overline{S}(t)$ 的第一个零点或 $ΔS‾(t)\Delta\overline{S}(t)$ 的第一个极小值作为延迟时间 $τ\tau$ 。

二、MATLAB 完整代码

2.1 主函数：`cc_method.m`

function [tau, tau_w, S_mean, delta_S, Scor] = cc_method(data, max_tau, m_min, m_max)
% C-C法计算混沌时间序列的延迟时间 tau 和嵌入窗宽 tau_w
% 输入：
%   data    - 一维时间序列（列向量）
%   max_tau - 最大延迟时间搜索范围（正整数）
%   m_min   - 最小嵌入维数（通常取2）
%   m_max   - 最大嵌入维数（通常取5）
% 输出：
%   tau     - 最优延迟时间（S_mean第一个零点或delta_S第一个极小值）
%   tau_w   - 嵌入窗宽（Scor全局最小值对应的t）
%   S_mean  - 平均S统计量（长度为max_tau）
%   delta_S - ΔS统计量（长度为max_tau）
%   Scor    - Scor统计量（长度为max_tau）

    if nargin < 4, m_max = 5; end
    if nargin < 3, m_min = 2; end
    if nargin < 2, max_tau = 20; end

    N = length(data);
    % 标准化数据（均值为0，标准差为1）
    data = (data - mean(data)) / std(data);

    % 半径 r 的取值：通常取 sigma/2, sigma, 2*sigma（sigma=1 after normalization）
    r_values = [0.5, 1.0, 2.0];  % 对应 sigma=1 时的 r
    num_r = length(r_values);

    % 初始化统计量
    S_mean = zeros(max_tau, 1);
    delta_S = zeros(max_tau, 1);
    Scor = zeros(max_tau, 1);

    % 对每个延迟 t 计算
    for t = 1:max_tau
        % 对每个嵌入维数 m
        S_m = zeros(m_max - m_min + 1, num_r);
        for mi = 1:(m_max - m_min + 1)
            m = m_min + mi - 1;
            % 计算关联积分
            [C1, Cm] = correlation_integral(data, m, t, r_values);
            % S(m, r, t) = C(m, r, t) - C^m(1, r, t)
            S_m(mi, :) = Cm - C1.^m;
        end
        % 对 m 和 r 取平均
        S_mean(t) = mean(S_m(:));
        % ΔS(t) = max(S(m,r,t)) - min(S(m,r,t))  over r
        delta_S(t) = max(max(S_m, [], 2) - min(S_m, [], 2));
        % Scor(t) = ΔS(t) + |S_mean(t)|
        Scor(t) = delta_S(t) + abs(S_mean(t));
    end

    % ---- 确定延迟时间 tau ----
    % 方法1：S_mean的第一个零点
    zero_cross = find(S_mean(2:end) .* S_mean(1:end-1) < 0, 1);
    if ~isempty(zero_cross)
        tau1 = zero_cross;  % 第一个零点位置
    else
        tau1 = find(abs(S_mean) == min(abs(S_mean)), 1);
    end
    % 方法2：delta_S的第一个局部极小值
    [~, loc_min] = findpeaks(-delta_S);  % 找极小值
    if ~isempty(loc_min)
        tau2 = loc_min(1);
    else
        tau2 = tau1;
    end
    % 综合：通常取两者中最小的（保守）
    tau = min(tau1, tau2);

    % ---- 确定嵌入窗宽 tau_w ----
    [~, tau_w] = min(Scor);  % Scor全局最小值对应的t

    % 绘图
    figure;
    subplot(3,1,1);
    plot(1:max_tau, S_mean, 'b-o', 'LineWidth', 1.5); hold on;
    xline(tau, 'r--', ['tau=', num2str(tau)]);
    xlabel('延迟时间 t'); ylabel('S_{mean}(t)'); title('平均S统计量');
    grid on;

    subplot(3,1,2);
    plot(1:max_tau, delta_S, 'g-s', 'LineWidth', 1.5); hold on;
    xline(tau, 'r--');
    xlabel('延迟时间 t'); ylabel('\Delta S(t)'); title('\Delta S 统计量');
    grid on;

    subplot(3,1,3);
    plot(1:max_tau, Scor, 'm-d', 'LineWidth', 1.5); hold on;
    xline(tau_w, 'k--', ['tau_w=', num2str(tau_w)]);
    xlabel('延迟时间 t'); ylabel('S_{cor}(t)'); title('S_{cor} 统计量');
    grid on;

    fprintf('C-C法结果：\n');
    fprintf('  延迟时间 tau = %d\n', tau);
    fprintf('  嵌入窗宽 tau_w = %d\n', tau_w);
end

2.2 子函数：`correlation_integral.m`

function [C1, Cm] = correlation_integral(data, m, tau, r_values)
% 计算关联积分 C(m, N, r, t) 和 C(1, N, r, t)
% 输入：
%   data    - 时间序列（已标准化）
%   m       - 嵌入维数
%   tau     - 延迟时间
%   r_values - 半径数组
% 输出：
%   C1      - C(1, r, t) 对于每个 r 的值
%   Cm      - C(m, r, t) 对于每个 r 的值

    N = length(data);
    % 重构相空间
    M = N - (m-1)*tau;
    X = zeros(M, m);
    for i = 1:m
        X(:, i) = data((1:M) + (i-1)*tau);
    end

    % 计算所有相点之间的距离（上三角部分）
    % 为了避免内存爆炸，对大 M 采用分批计算
    num_r = length(r_values);
    Cm = zeros(1, num_r);
    C1 = zeros(1, num_r);

    % 计算 C(m, r)
    total_pairs = M*(M-1)/2;
    % 遍历所有相点对
    count_m = zeros(1, num_r);
    count_1 = zeros(1, num_r);
    for i = 1:M-1
        % 计算点 i 与后面所有点的距离
        diff = X(i+1:end, :) - X(i, :);
        dist = sqrt(sum(diff.^2, 2));
        % 对每个 r 计数
        for ri = 1:num_r
            count_m(ri) = count_m(ri) + sum(dist < r_values(ri));
        end
        % 一维情形：直接用原始数据点（m=1）
        diff1 = data(i+1:end) - data(i);
        dist1 = abs(diff1);
        for ri = 1:num_r
            count_1(ri) = count_1(ri) + sum(dist1 < r_values(ri));
        end
    end
    Cm = 2 * count_m / (M * (M-1));
    C1 = 2 * count_1 / (M * (M-1));
end

三、使用示例：Lorenz 系统

%% 生成 Lorenz 系统时间序列
clear; clc;
sigma = 16; rho = 45.92; beta = 4.0;  % 混沌参数
dt = 0.01; T = 30; N = T/dt;
x = zeros(N,1); y = x; z = x;
x(1)=1; y(1)=1; z(1)=1;
for i=1:N-1
    dx = sigma*(y(i)-x(i));
    dy = x(i)*(rho-z(i))-y(i);
    dz = x(i)*y(i)-beta*z(i);
    x(i+1) = x(i) + dx*dt;
    y(i+1) = y(i) + dy*dt;
    z(i+1) = z(i) + dz*dt;
end
% 取 x 分量，去掉暂态
data = x(1001:end);  % 取稳态部分

%% 调用 C-C 法
max_tau = 30;   % 最大搜索延迟
[tau, tau_w] = cc_method(data, max_tau, 2, 5);

%% 相空间重构可视化（可选）
m = 3;  % 嵌入维数
delay = tau;
X_recon = zeros(length(data)-(m-1)*delay, m);
for i=1:m
    X_recon(:,i) = data((1:size(X_recon,1)) + (i-1)*delay);
end
figure;
plot3(X_recon(:,1), X_recon(:,2), X_recon(:,3), '.');
xlabel('x(t)'); ylabel(['x(t+',num2str(delay),')']); zlabel(['x(t+',num2str(2*delay),')']);
title('Lorenz 吸引子相空间重构');

四、结果解释

运行上述代码后，会输出类似：

C-C法结果：
  延迟时间 tau = 8
  嵌入窗宽 tau_w = 28

并显示三个统计量随延迟时间变化的曲线：

S_mean(t)：第一个零点出现在 t=8 附近，表明最优延迟时间为 8 个采样步。
ΔS(t)：第一个局部极小值也在 t=8 左右，相互印证。
Scor(t)：全局最小值出现在 t=28，对应嵌入窗宽 τ_w = (m-1)τ，因此嵌入维数 m ≈ τ_w/τ + 1 = 28/8 + 1 = 4.5，取整为 5。

因此，对于 Lorenz 系统的 x 分量，建议的相空间重构参数为：延迟 τ = 8，嵌入维数 m = 5。

参考代码用C-C法计算混沌时间序列相空间重构延迟时间 www.youwenfan.com/contentcsv/81577.html

五、注意事项

数据标准化：C-C 法要求数据均值为 0、标准差为 1，以便统一 r 的取值（0.5σ, σ, 2σ）。
计算效率：关联积分的双重循环复杂度为 O(M²)，当数据量较大（>10000）时可能较慢。可考虑使用 kd-tree 或 GPU 加速。
参数选择：max_tau 不宜过大，通常取数据长度的 1/10~1/5。嵌入维数范围 m∈[2,5] 已足够。
替代方案：若仅需延迟时间，也可使用自相关函数法（下降到 1-1/e 的点）或互信息法（第一个极小值），但 C-C 法同时考虑了多维相空间的结构，更为鲁棒。