混沌时间序列参数计算工具箱
混沌时间序列参数计算工具箱，包含时间延迟、嵌入维数、关联维等关键参数的计算程序。
混沌时间序列参数计算工具箱
1. 主函数 - 混沌参数综合分析
function chaos_analysis_results = chaotic_time_series_analysis(time_series, varargin) % 混沌时间序列参数综合分析 % 输入: % time_series - 时间序列数据 % varargin - 可选参数 % 输出: % chaos_analysis_results - 包含所有分析结果的结构体p = inputParser; addRequired(p, 'time_series', @isnumeric); addParameter(p, 'max_delay', 50, @isnumeric); addParameter(p, 'max_dim', 10, @isnumeric); addParameter(p, 'radius_range', logspace(-3, 0, 50), @isnumeric); addParameter(p, 'show_plots', true, @islogical); parse(p, time_series, varargin{:}); fprintf('开始混沌时间序列分析...\n'); fprintf('数据长度: %d\n', length(time_series)); % 初始化结果结构 results = struct(); %% 1. 时间延迟计算 fprintf('计算时间延迟...\n'); [results.tau_autocorr, results.autocorr_func] = calculate_time_delay_autocorr(time_series, p.Results.max_delay); [results.tau_mutual, results.mutual_info] = calculate_time_delay_mutual(time_series, p.Results.max_delay); %% 2. 嵌入维数计算 fprintf('计算嵌入维数...\n'); [results.embedding_dim, results.fnn_ratio] = calculate_embedding_dimension(time_series, results.tau_mutual, p.Results.max_dim); %% 3. 关联维计算 fprintf('计算关联维...\n'); [results.correlation_dim, results.correlation_integral] = calculate_correlation_dimension(time_series, results.tau_mutual, results.embedding_dim, p.Results.radius_range); %% 4. Lyapunov指数计算 fprintf('计算最大Lyapunov指数...\n'); [results.lyapunov_exponent, results.lyapunov_curve] = calculate_lyapunov_exponent(time_series, results.tau_mutual, results.embedding_dim); %% 5. 其他混沌指标 fprintf('计算其他混沌指标...\n'); results.hurst_exponent = calculate_hurst_exponent(time_series); [results.approximate_entropy, results.sample_entropy] = calculate_entropy_measures(time_series); %% 6. 可视化结果 if p.Results.show_plots plot_chaos_analysis_results(results, time_series); end %% 7. 生成报告 generate_analysis_report(results); chaos_analysis_results = results; fprintf('混沌分析完成!\n');
end
2. 时间延迟计算
2.1 自相关函数法
function [optimal_tau, autocorr_values] = calculate_time_delay_autocorr(time_series, max_delay) % 使用自相关函数法计算时间延迟 % 选择自相关函数第一次过零或下降到1-1/e时的延迟n = length(time_series); autocorr_values = zeros(1, max_delay); % 计算自相关函数 for tau = 1:max_delay autocorr_values(tau) = abs(corr(time_series(1:end-tau)', time_series(1+tau:end)')); end % 寻找最优延迟 threshold = 1 - 1/exp(1); % 约0.6321 % 方法1: 第一次低于阈值 below_threshold = find(autocorr_values < threshold, 1); if ~isempty(below_threshold) optimal_tau = below_threshold; else % 方法2: 选择最小值 [~, optimal_tau] = min(autocorr_values); end fprintf('自相关函数法 - 最优时间延迟: %d\n', optimal_tau);
end
2.2 互信息法
function [optimal_tau, mutual_info_values] = calculate_time_delay_mutual(time_series, max_delay) % 使用互信息法计算时间延迟 % 选择互信息第一次达到局部最小值时的延迟n = length(time_series); mutual_info_values = zeros(1, max_delay); % 将时间序列离散化 num_bins = 10; % 直方图箱数 [~, edges] = histcounts(time_series, num_bins); % 计算原始序列的概率分布 px = histcounts(time_series, edges, 'Normalization', 'probability'); for tau = 1:max_delay % 延迟序列 x_t = time_series(1:end-tau); x_t_tau = time_series(1+tau:end); % 计算联合概率分布 joint_hist = histcounts2(x_t, x_t_tau, edges, edges, 'Normalization', 'probability'); py = histcounts(x_t_tau, edges, 'Normalization', 'probability'); % 计算互信息 mi = 0; for i = 1:num_bins for j = 1:num_bins if joint_hist(i,j) > 0 && px(i) > 0 && py(j) > 0 mi = mi + joint_hist(i,j) * log(joint_hist(i,j) / (px(i) * py(j))); end end end mutual_info_values(tau) = mi; end % 寻找第一个局部最小值 optimal_tau = find_first_local_minimum(mutual_info_values); if isempty(optimal_tau) optimal_tau = 1; % 默认值 end fprintf('互信息法 - 最优时间延迟: %d\n', optimal_tau);end
function min_idx = find_first_local_minimum(data)
% 寻找第一个局部最小值
for i = 2:length(data)-1
if data(i) < data(i-1) && data(i) < data(i+1)
min_idx = i;
return;
end
end
min_idx = [];
end
3. 嵌入维数计算
function [optimal_dim, fnn_ratio] = calculate_embedding_dimension(time_series, tau, max_dim) % 使用伪最近邻法计算嵌入维数 % 当伪最近邻比例低于阈值时的最小维数即为最优嵌入维数n = length(time_series); fnn_ratio = zeros(1, max_dim); % 计算每个嵌入维数的伪最近邻比例 for m = 1:max_dim % 重构相空间 embedded_data = phase_space_reconstruction(time_series, m, tau); % 寻找最近邻 [fnn_ratio(m), ~] = calculate_fnn_ratio(embedded_data, m); fprintf('嵌入维数 %d: 伪最近邻比例 = %.4f\n', m, fnn_ratio(m)); end % 选择伪最近邻比例低于阈值的最小维数 threshold = 0.1; % 常用阈值 below_threshold = find(fnn_ratio < threshold, 1); if ~isempty(below_threshold) optimal_dim = below_threshold; else % 选择变化平缓的点 diff_ratio = diff(fnn_ratio); if any(diff_ratio < 0.01) optimal_dim = find(diff_ratio < 0.01, 1) + 1; else optimal_dim = max_dim; end end fprintf('最优嵌入维数: %d\n', optimal_dim);end
function embedded_data = phase_space_reconstruction(time_series, m, tau)
% 相空间重构
n = length(time_series);
embedded_length = n - (m-1)*tau;embedded_data = zeros(embedded_length, m); for i = 1:embedded_length for j = 1:m embedded_data(i, j) = time_series(i + (j-1)*tau); end endend
function [fnn_ratio, neighbors] = calculate_fnn_ratio(embedded_data, m)
% 计算伪最近邻比例[n_points, ~] = size(embedded_data); % 寻找每个点的最近邻 distances = zeros(n_points, 1); neighbor_indices = zeros(n_points, 1); for i = 1:n_points % 计算到其他点的距离（排除自身） dist = zeros(n_points, 1); for j = 1:n_points if i ~= j dist(j) = norm(embedded_data(i, :) - embedded_data(j, :)); else dist(j) = inf; end end [min_dist, min_idx] = min(dist); distances(i) = min_dist; neighbor_indices(i) = min_idx; end % 计算伪最近邻比例 fnn_count = 0; threshold = 10; % 常用阈值 for i = 1:n_points current_point = embedded_data(i, :); neighbor_point = embedded_data(neighbor_indices(i), :); % 计算距离比 if m > 1 % 在m-1维空间中的距离 dist_m_minus_1 = norm(current_point(1:end-1) - neighbor_point(1:end-1)); dist_m = distances(i); if dist_m_minus_1 > 0 ratio = abs(dist_m - dist_m_minus_1) / dist_m_minus_1; if ratio > threshold fnn_count = fnn_count + 1; end end end end fnn_ratio = fnn_count / n_points; neighbors = neighbor_indices;
end
4. 关联维计算
function [correlation_dim, corr_integral] = calculate_correlation_dimension(time_series, tau, m, radius_range) % 使用G-P算法计算关联维% 重构相空间 embedded_data = phase_space_reconstruction(time_series, m, tau); [n_points, ~] = size(embedded_data); corr_integral = zeros(length(radius_range), 1); % 计算关联积分 for r_idx = 1:length(radius_range) r = radius_range(r_idx); count = 0; % 使用随机采样来加速计算 sample_size = min(1000, n_points); sample_indices = randperm(n_points, sample_size); for i = 1:sample_size idx_i = sample_indices(i); for j = i+1:sample_size idx_j = sample_indices(j); distance_ij = norm(embedded_data(idx_i, :) - embedded_data(idx_j, :)); if distance_ij < r count = count + 1; end end end % 计算关联积分 corr_integral(r_idx) = 2 * count / (sample_size * (sample_size - 1)); end % 寻找线性区域并计算关联维 [correlation_dim, scaling_region] = estimate_correlation_dimension(radius_range, corr_integral); fprintf('关联维数: %.4f\n', correlation_dim);end
function [correlation_dim, scaling_region] = estimate_correlation_dimension(radius, corr_integral)
% 估计关联维数 - 通过线性区域的斜率% 对数变换 log_r = log10(radius); log_c = log10(corr_integral); % 移除无效点（零或负值） valid_idx = isfinite(log_r) & isfinite(log_c) & log_c > -inf; log_r_valid = log_r(valid_idx); log_c_valid = log_c(valid_idx); if length(log_r_valid) < 10 correlation_dim = 0; scaling_region = []; return; end % 寻找线性区域（通过滑动窗口） window_size = 10; max_slope = -inf; best_start = 1; for i = 1:length(log_r_valid)-window_size+1 window_r = log_r_valid(i:i+window_size-1); window_c = log_c_valid(i:i+window_size-1); % 线性拟合 p = polyfit(window_r, window_c, 1); slope = p(1); % 检查拟合质量 residuals = window_c - polyval(p, window_r); r_squared = 1 - sum(residuals.^2) / sum((window_c - mean(window_c)).^2); if slope > max_slope && r_squared > 0.9 max_slope = slope; best_start = i; end end % 使用最佳线性区域计算关联维 scaling_region = best_start:best_start+window_size-1; p_final = polyfit(log_r_valid(scaling_region), log_c_valid(scaling_region), 1); correlation_dim = p_final(1);
end
5. Lyapunov指数计算
function [lyapunov_exp, divergence_curve] = calculate_lyapunov_exponent(time_series, tau, m) % 计算最大Lyapunov指数 - 使用小数据量法% 重构相空间 embedded_data = phase_space_reconstruction(time_series, m, tau); [n_points, ~] = size(embedded_data); % 参数设置 min_neighbors = 5; evolution_steps = min(50, floor(n_points/10)); % 为每个点寻找最近邻 divergence = zeros(evolution_steps, 1); count = zeros(evolution_steps, 1); for i = 1:n_points-evolution_steps % 当前点 current_point = embedded_data(i, :); % 寻找最近邻（排除太近的点） distances = zeros(n_points, 1); for j = 1:n_points if abs(i - j) > evolution_steps % 避免时间相关 distances(j) = norm(current_point - embedded_data(j, :)); else distances(j) = inf; end end [~, neighbor_indices] = mink(distances, min_neighbors); % 追踪每个最近邻的演化 for k = 1:min_neighbors neighbor_idx = neighbor_indices(k); initial_distance = norm(embedded_data(i, :) - embedded_data(neighbor_idx, :)); if initial_distance > 0 for step = 1:min(evolution_steps, n_points-max(i, neighbor_idx)) d_i = i + step; d_j = neighbor_idx + step; if d_j <= n_points current_distance = norm(embedded_data(d_i, :) - embedded_data(d_j, :)); if current_distance > 0 divergence(step) = divergence(step) + log(current_distance / initial_distance); count(step) = count(step) + 1; end end end end end end % 计算平均发散率 valid_steps = count > 0; divergence_curve = zeros(size(divergence)); divergence_curve(valid_steps) = divergence(valid_steps) ./ count(valid_steps); % 线性拟合得到Lyapunov指数 time_steps = (1:evolution_steps)'; valid_range = find(valid_steps & time_steps > 0 & time_steps < evolution_steps/2); if length(valid_range) > 5 p = polyfit(time_steps(valid_range), divergence_curve(valid_range), 1); lyapunov_exp = p(1); else lyapunov_exp = 0; end fprintf('最大Lyapunov指数: %.6f\n', lyapunov_exp);
end
6. 其他混沌指标
function hurst_exp = calculate_hurst_exponent(time_series) % 计算Hurst指数 - 重标极差法(R/S方法)n = length(time_series); max_window = floor(n/10); window_sizes = unique(round(logspace(log10(10), log10(max_window), 20))); rs_values = zeros(length(window_sizes), 1); for i = 1:length(window_sizes) window_size = window_sizes(i); num_windows = floor(n / window_size); window_rs = zeros(num_windows, 1); for j = 1:num_windows start_idx = (j-1)*window_size + 1; end_idx = j*window_size; window_data = time_series(start_idx:end_idx); % 计算累积偏差 mean_val = mean(window_data); deviations = window_data - mean_val; cumulative_deviations = cumsum(deviations); % 范围 R = max(cumulative_deviations) - min(cumulative_deviations); S = std(window_data); if S > 0 window_rs(j) = R / S; else window_rs(j) = 0; end end rs_values(i) = mean(window_rs(window_rs > 0)); end % 线性拟合得到Hurst指数 valid_idx = rs_values > 0 & window_sizes' > 0; if sum(valid_idx) > 3 p = polyfit(log(window_sizes(valid_idx)), log(rs_values(valid_idx)), 1); hurst_exp = p(1); else hurst_exp = 0.5; end fprintf('Hurst指数: %.4f\n', hurst_exp);end
function [approx_entropy, sample_entropy] = calculate_entropy_measures(time_series)
% 计算近似熵和样本熵m = 2; % 嵌入维数 r = 0.2 * std(time_series); % 容忍度 % 近似熵 approx_entropy = approximate_entropy(time_series, m, r); % 样本熵 sample_entropy = sample_entropy_calc(time_series, m, r); fprintf('近似熵: %.4f, 样本熵: %.4f\n', approx_entropy, sample_entropy);end
function apen = approximate_entropy(data, m, r)
% 计算近似熵
n = length(data);% 重构相空间 embedded_m = zeros(n-m+1, m); for i = 1:n-m+1 embedded_m(i, :) = data(i:i+m-1); end % 计算距离矩阵 phi_m = calculate_phi(embedded_m, n-m+1, r); % 对于m+1 embedded_m1 = zeros(n-m, m+1); for i = 1:n-m embedded_m1(i, :) = data(i:i+m); end phi_m1 = calculate_phi(embedded_m1, n-m, r); % 计算近似熵 apen = log(phi_m / phi_m1);end
function phi = calculate_phi(embedded_data, n, r)
% 计算Phi函数
count = zeros(n, 1);for i = 1:n for j = 1:n if i ~= j if max(abs(embedded_data(i, :) - embedded_data(j, :))) <= r count(i) = count(i) + 1; end end end end C = count / (n-1); phi = mean(C);end
function sampen = sample_entropy_calc(data, m, r)
% 计算样本熵
n = length(data);% 重构相空间 - m维 embedded_m = zeros(n-m, m); for i = 1:n-m embedded_m(i, :) = data(i:i+m-1); end % 计算匹配数 - m维 A = 0; for i = 1:n-m for j = i+1:n-m if max(abs(embedded_m(i, :) - embedded_m(j, :))) <= r A = A + 1; end end end % 重构相空间 - m+1维 embedded_m1 = zeros(n-m-1, m+1); for i = 1:n-m-1 embedded_m1(i, :) = data(i:i+m); end % 计算匹配数 - m+1维 B = 0; for i = 1:n-m-1 for j = i+1:n-m-1 if max(abs(embedded_m1(i, :) - embedded_m1(j, :))) <= r B = B + 1; end end end % 计算样本熵 if B > 0 && A > 0 sampen = -log(B / A); else sampen = 0; end
end
7. 可视化函数
function plot_chaos_analysis_results(results, time_series) % 绘制混沌分析结果figure('Position', [100, 100, 1400, 1000]); % 原始时间序列 subplot(3, 4, 1); plot(time_series, 'b-', 'LineWidth', 1); title('原始时间序列'); xlabel('时间'); ylabel('幅值'); grid on; % 自相关函数 subplot(3, 4, 2); plot(1:length(results.autocorr_func), results.autocorr_func, 'ro-', 'LineWidth', 1.5); hold on; plot([results.tau_autocorr, results.tau_autocorr], [0, 1], 'r--', 'LineWidth', 2); title('自相关函数法'); xlabel('延迟'); ylabel('自相关'); legend('自相关', '最优延迟', 'Location', 'best'); grid on; % 互信息 subplot(3, 4, 3); plot(1:length(results.mutual_info), results.mutual_info, 'go-', 'LineWidth', 1.5); hold on; plot([results.tau_mutual, results.tau_mutual], [0, max(results.mutual_info)], 'g--', 'LineWidth', 2); title('互信息法'); xlabel('延迟'); ylabel('互信息'); legend('互信息', '最优延迟', 'Location', 'best'); grid on; % 伪最近邻 subplot(3, 4, 4); plot(1:length(results.fnn_ratio), results.fnn_ratio, 'bo-', 'LineWidth', 1.5); hold on; plot([results.embedding_dim, results.embedding_dim], [0, max(results.fnn_ratio)], 'b--', 'LineWidth', 2); title('伪最近邻法'); xlabel('嵌入维数'); ylabel('伪最近邻比例'); legend('FNN比例', '最优维数', 'Location', 'best'); grid on; % 关联维 subplot(3, 4, 5); loglog(results.correlation_integral.radius, results.correlation_integral.C, 'mo-', 'LineWidth', 1.5); title('关联积分'); xlabel('距离 r'); ylabel('关联积分 C(r)'); grid on; % 双对数图 subplot(3, 4, 6); plot(log10(results.correlation_integral.radius), log10(results.correlation_integral.C), 'mo-', 'LineWidth', 1.5); title('关联维计算'); xlabel('log(r)'); ylabel('log(C(r))'); grid on; % Lyapunov指数 subplot(3, 4, 7); plot(1:length(results.lyapunov_curve), results.lyapunov_curve, 'co-', 'LineWidth', 1.5); title('Lyapunov指数计算'); xlabel('演化步数'); ylabel('平均发散率'); grid on; % 相空间重构 (3D) if results.embedding_dim >= 3 subplot(3, 4, 8); embedded_3d = phase_space_reconstruction(time_series, 3, results.tau_mutual); plot3(embedded_3d(:,1), embedded_3d(:,2), embedded_3d(:,3), 'b-', 'LineWidth', 0.5); title('相空间重构 (3D)'); grid on; end % 结果汇总 subplot(3, 4, [9, 10, 11, 12]); text(0.1, 0.9, sprintf('混沌参数分析结果汇总:'), 'FontSize', 12, 'FontWeight', 'bold'); text(0.1, 0.7, sprintf('时间延迟 (自相关法): %d', results.tau_autocorr), 'FontSize', 10); text(0.1, 0.6, sprintf('时间延迟 (互信息法): %d', results.tau_mutual), 'FontSize', 10); text(0.1, 0.5, sprintf('嵌入维数: %d', results.embedding_dim), 'FontSize', 10); text(0.1, 0.4, sprintf('关联维数: %.4f', results.correlation_dim), 'FontSize', 10); text(0.1, 0.3, sprintf('最大Lyapunov指数: %.6f', results.lyapunov_exponent), 'FontSize', 10); text(0.1, 0.2, sprintf('Hurst指数: %.4f', results.hurst_exponent), 'FontSize', 10); text(0.1, 0.1, sprintf('近似熵: %.4f, 样本熵: %.4f', results.approximate_entropy, results.sample_entropy), 'FontSize', 10); % 混沌特性判断 if results.lyapunov_exponent > 0 chaos_status = '具有混沌特性'; else chaos_status = '无显著混沌特性'; end text(0.5, 0.9, sprintf('系统状态: %s', chaos_status), 'FontSize', 12, 'Color', 'red', 'FontWeight', 'bold'); axis off;end
function generate_analysis_report(results)
% 生成分析报告fprintf('\n=== 混沌时间序列分析报告 ===\n'); fprintf('1. 时间延迟参数:\n'); fprintf(' - 自相关函数法: τ = %d\n', results.tau_autocorr); fprintf(' - 互信息法: τ = %d\n', results.tau_mutual); fprintf('2. 嵌入参数:\n'); fprintf(' - 最优嵌入维数: m = %d\n', results.embedding_dim); fprintf('3. 分形特性:\n'); fprintf(' - 关联维数: D₂ = %.4f\n', results.correlation_dim); fprintf(' - Hurst指数: H = %.4f\n', results.hurst_exponent); fprintf('4. 动力学特性:\n'); fprintf(' - 最大Lyapunov指数: λ₁ = %.6f\n', results.lyapunov_exponent); fprintf(' - 近似熵: %.4f\n', results.approximate_entropy); fprintf(' - 样本熵: %.4f\n', results.sample_entropy); % 混沌判断 fprintf('5. 混沌特性判断:\n'); if results.lyapunov_exponent > 0 fprintf(' - 系统具有混沌特性 (λ₁ > 0)\n'); if results.correlation_dim > 0 && ~isinf(results.correlation_dim) fprintf(' - 系统具有分形吸引子\n'); end else fprintf(' - 系统无显著混沌特性\n'); end fprintf('================================\n\n');
end
8. 使用示例和测试
function chaos_analysis_demo() % 混沌时间序列分析演示fprintf('混沌时间序列参数计算工具箱演示\n'); % 生成测试数据 fprintf('生成测试时间序列...\n'); % 1. Lorenz系统 lorenz_data = generate_lorenz_system(1000); % 2. 混沌逻辑映射 logistic_data = generate_logistic_map(0.4, 3.9, 1000); % 3. 随机序列 random_data = randn(1000, 1); % 分析Lorenz系统 fprintf('\n分析Lorenz系统:\n'); results_lorenz = chaotic_time_series_analysis(lorenz_data(:,1), 'max_delay', 30, 'max_dim', 8); % 分析逻辑映射 fprintf('\n分析混沌逻辑映射:\n'); results_logistic = chaotic_time_series_analysis(logistic_data, 'max_delay', 20, 'max_dim', 6); % 分析随机序列 fprintf('\n分析随机序列:\n'); results_random = chaotic_time_series_analysis(random_data, 'max_delay', 20, 'max_dim', 6); % 比较结果 compare_results(results_lorenz, results_logistic, results_random);end
function lorenz_data = generate_lorenz_system(n_points)
% 生成Lorenz系统数据
sigma = 10;
rho = 28;
beta = 8/3;dt = 0.01; x = zeros(n_points, 1); y = zeros(n_points, 1); z = zeros(n_points, 1); % 初始条件 x(1) = 1; y(1) = 1; z(1) = 1; % 数值积分 for i = 1:n_points-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 lorenz_data = [x, y, z];end
function logistic_data = generate_logistic_map(x0, r, n_points)
% 生成逻辑映射数据
logistic_data = zeros(n_points, 1);
logistic_data(1) = x0;for i = 2:n_points logistic_data(i) = r * logistic_data(i-1) * (1 - logistic_data(i-1)); endend
function compare_results(results1, results2, results3)
% 比较不同系统的分析结果fprintf('\n=== 不同系统混沌特性比较 ===\n'); fprintf('系统\t\tLyapunov指数\t关联维数\tHurst指数\t混沌特性\n'); fprintf('----------------------------------------------------------------\n'); systems = {'Lorenz系统', '逻辑映射', '随机序列'}; results = {results1, results2, results3}; for i = 1:3 lyap = results{i}.lyapunov_exponent; corr_dim = results{i}.correlation_dim; hurst = results{i}.hurst_exponent; if lyap > 0 chaos_char = '是'; else chaos_char = '否'; end fprintf('%-12s\t%.6f\t\t%.4f\t\t%.4f\t\t%s\n', ... systems{i}, lyap, corr_dim, hurst, chaos_char); end
end
参考工具箱 含有各种混沌时间序列参数（时间延迟，嵌入维数，关联维等）计算程序 www.3dddown.com/cna/82714.html
工具箱特性总结
| 功能模块 | 计算方法 | 输出参数 |
|---|---|---|
| 时间延迟 | 自相关函数法、互信息法 | 最优延迟τ |
| 嵌入维数 | 伪最近邻法 | 最优嵌入维数m |
| 关联维数 | G-P算法 | 关联维D₂ |
| Lyapunov指数 | 小数据量法 | 最大Lyapunov指数λ₁ |
| 其他指标 | R/S法、熵计算 | Hurst指数、近似熵、样本熵 |
浙公网安备 33010602011771号