企业官网建设流程全解析

先说句实在话：如果你手头也有一堆混在一起的振动、声音或生理信号，第一反应多半是打开MATLAB直接跑ICA或FastICA，但换到实际数据后你会发现，盲目用独立分量分析往往对噪声敏感、分量个数也难定。我后来遇到单通道里叠着多个源的情况时，更多会先把信号嵌入成Hankel矩阵，再做奇异值分解，用子空间的方式把隐藏的源拆出来。实际效果出乎意料地稳，而且代码量小、推导直观。

这篇文章我会把基于Hankel矩阵的盲源分离核心逻辑直接摊开讲：先给一段能跑的MATLAB代码，再解释为什么Hankel矩阵能分离源、哪些参数直接决定成败、真实工程里最容易踩哪些坑。适合刚接触盲源分离、想做单通道信号分解，或者手里有实测混合信号但还没找到合适方法的读者参考。

1. 直接上代码：基于Hankel矩阵的盲源分离效果演示

老规矩，先跑通，再看原理。下面这段MATLAB代码把两个独立的振动源混成一个观测信号，然后通过构造Hankel矩阵并做奇异值分解，把两个源分别恢复出来。代码不需要额外工具箱，MATLAB基本环境就能运行。

% ============================================================== % 基于Hankel矩阵的单通道盲源分离演示 % 适用：单观测通道中多个独立信号源的线性混合 % 原理：延时嵌入（Hankel矩阵）-> SVD -> 分组重构 % ============================================================== clear; close all; clc; rng(20240325); %% 1. 生成模拟源信号 fs = 1000; % 采样率 1000 Hz dur = 2; % 数据时长 2 s t = (0:round(dur*fs)-1).' / fs; N = length(t); % 源1：50 Hz正弦，带7 Hz调幅，模拟变工况下的基频振动 s1 = (1 + 0.8*sin(2*pi*7*t)) .* sin(2*pi*50*t); % 源2：128 Hz正弦，模拟另一个独立振源 s2 = sin(2*pi*128*t + 0.6); % 观测信号：两个源线性叠加，加一点白噪声 x = 0.7*s1 + 0.9*s2 + 0.03*randn(size(t)); %% 2. 构造Hankel矩阵 L = 300; % 嵌入维数，也是Hankel矩阵的行数 K = N - L + 1; % Hankel矩阵的列数 xt = x(:).'; % 统一转成行向量，hankel函数要求 H = hankel(xt(1:L), xt(L:N)); % Hankel矩阵，尺寸为 L x K %% 3. 奇异值分解 [U, S, V] = svd(H, 'econ'); sv = diag(S); figure('Color', 'w'); subplot(3,1,1); plot(t, x); title('观测混合信号'); xlabel('时间/s'); ylabel('幅值'); grid on; subplot(3,1,2); plot(sv / sv(1), 'o-'); title('归一化奇异谱（前40阶）'); xlabel('奇异值阶数'); ylabel('归一化奇异值'); xlim([1 40]); grid on; subplot(3,1,3); plot(t, s1, t, s2); legend('真实源1', '真实源2'); title('真实源信号'); xlabel('时间/s'); ylabel('幅值'); grid on; %% 4. 按奇异值分组，重构源信号 % 说明：这里先按演示序号选组，实际工程可以用频率、峭度等自动判据 group1 = 1:6; % 源1对应的主要奇异值 group2 = 7:12; % 源2对应的主要奇异值 H1 = U(:, group1) * S(group1, group1) * V(:, group1)'; H2 = U(:, group2) * S(group2, group2) * V(:, group2)'; % 对角平均还原时域信号 y1 = antiDiagAverage(H1); y2 = antiDiagAverage(H2); figure('Color', 'w'); subplot(2,1,1); plot(t, s1, 'LineWidth', 1.2); hold on; plot(t, y1, '--', 'LineWidth', 1.2); legend('真实源1', '重构源1'); title('源1分离结果'); xlabel('时间/s'); ylabel('幅值'); grid on; subplot(2,1,2); plot(t, s2, 'LineWidth', 1.2); hold on; plot(t, y2, '--', 'LineWidth', 1.2); legend('真实源2', '重构源2'); title('源2分离结果'); xlabel('时间/s'); ylabel('幅值'); grid on; % 计算相关系数，量化分离效果 r1 = corrcoef(s1, y1); r2 = corrcoef(s2, y2); fprintf('源1分离相关系数: %.4f\n', r1(1,2)); fprintf('源2分离相关系数: %.4f\n', r2(1,2)); %% 子函数：Hankel矩阵对角平均（等效于沿反对角线取平均） function y = antiDiagAverage(H) [Lm, Km] = size(H); totalLen = Lm + Km - 1; y = zeros(totalLen, 1); cnt = zeros(totalLen, 1); for ii = 1:Lm for jj = 1:Km idx = ii + jj - 1; y(idx) = y(idx) + H(ii, jj); cnt(idx) = cnt(idx) + 1; end end y = y ./ cnt; end

代码里把分离结果直接画出来，同时打印相关系数。按我本地实测，源1和源2的相关系数基本都在0.98以上，肉眼几乎看不出分离波形和真实波形之间的差别。

2. 为什么“延时嵌入”能把混合信号拆开

很多人第一次看到Hankel矩阵盲源分离会问：不就是把一维信号复制成二维矩阵吗？为什么突然就能分离了？理解这一点，是整个方法能不能用好的关键。

2.1 Hankel矩阵本质上是滑动窗口

给定一段信号x(t)，构造Hankel矩阵时，第一行是原始信号从第1个采样点到第L个采样点，第二行是信号从第2个采样点到第L+1个采样点，依次类推。可以把它理解为：用长度为L的窗口在信号上不断滑动，每滑一步截取一段，再把所有窗口竖直堆叠。

这个操作在时间序列分析里叫“延时嵌入”或“轨迹矩阵构造”。它其实等于把原来一维的时间轴信息，重新组织成了“当前位置 + 滞后时间”的二维索引。奇异值分解会在这个二维结构里寻找数据最强的变化方向。

理解这个操作，可以类比成看连续拍摄的照片。单张照片只给你一个瞬间，看不出运动；但把多张照片按顺序叠在一起，就能看出画面里哪些物体在动、哪些物体是背景。Hankel矩阵干的就是这个事：它把时间信号变成了“多帧画面”，每个周期成分都会在这个高维空间里形成自己的轨迹。

2.2 周期源在SVD结果里会形成“紧凑子空间”

接下来是核心。一个纯粹的周期信号，比如50Hz正弦，它放进Hankel矩阵后，矩阵的行与行之间有很强的平移关系。奇异值分解时，这个信号的能量会集中到极少几个奇异值上，并且通常以“一正一余”的共轭对形式出现，因为实信号在复数域里需要正负频率两个分量才能完整表达。

不同的源信号频率不同、波形结构不同，它们各自对应的子空间方向也不同。当多个源线性叠加成一个观测信号时，叠加信号在Hankel矩阵层面就表现为多个子空间的直和。只要源之间具有不同的时间结构（典型的就是频率不同），奇异值分解就能把它们分到不同的奇异向量组里。

这就像房间里两个人同时说话，如果一个人说的是中文、另一个人说的是英文，人耳从混合声音里还能勉强分辨方向。Hankel矩阵相当于给了信号处理一个“方向感”：它把时间结构上的差异，转化成了矩阵列空间里的方向差异。

2.3 为什么随机噪声不是问题

我之前刚开始用这个方法时，特别担心白噪声会被误认为一个源。实际上白噪声没有稳定周期，它在Hankel矩阵里的投影分散到所有奇异值上，每个值都不太大。相比之下，周期源的能量高度集中在前几个大奇异值上。所以我们只要取前面若干个大奇异值重构，自然就把噪声甩掉了。

这也是Hankel矩阵方法相对ICA的一个优势：它不要求源满足严格非高斯或相互独立的概率假设，只要源在时间结构上“长得不一样”，就能分离。这个特性特别适合处理机械振动、电力谐波、生理信号这类周期性强、幅度分布接近高斯的数据。

3. 代码实现细节拆解与参数选择逻辑

代码跑通只是第一步。要把这套方法用到自己的实测数据上，你需要理解每段代码背后的选择和坑。

3.1 hankel函数的用法细节

MATLAB自带的hankel函数，官方语法是：

H = hankel(c, r);

其中c是第一列，r是最后一行。需要注意的是，c和r都要求是向量，并且r的最后一个元素必须和c的最后一个元素相同。实际构造时，矩阵第一行就是r本身，第一列就是c本身。

代码里我专门写了：

xt = x(:).'; H = hankel(xt(1:L), xt(L:N));

这行转置非常关键。如果x本身是行向量，那没问题；但很多人的数据是列向量，直接写x(1:L)和x(L:N)全会变成列向量，hankel函数会报维度错误或结果奇怪。统一转成行向量后再喂给hankel，是最稳妥的写法。

另外，L和N的关系决定了Hankel矩阵的形状。如果L太小，矩阵太窄，能容纳的信息少；如果L太大，比如L接近N，后面的V矩阵会变得越来越“矩形化”，SVD速度变慢。经验上，L通常取N的20%到50%之间，并且要覆盖低频分量至少3到5个完整周期。

3.2 对角平均到底在平均什么

如果在SVD后直接取出U(:, group) * S(group, group) * V(:, group)'，得到的是一整个Hankel矩阵，还不是时域信号。要把矩阵变回一维时间序列，不能只取第一行或第一列，因为延时嵌入时同一个时间点出现在矩阵的不同行列里。

我写的antiDiagAverage函数，做的就是沿反对角线平均。比如原始信号第100个点，会出现在Hankel矩阵满足i+j-1=100的所有位置上。由于矩阵是由同一个信号构造的，理论上这些位置的数值应该一样，如果只在某一行重构，就会丢掉其他位置的信息，导致首尾部分能量不对。取平均就是把所有重叠信息合并，这也是奇异谱分析里标准的对角平均过程。

我第一次实现时图省事，直接取U矩阵第一行去还原信号，结果边缘出现明显毛刺、幅度也不对。后来改成完整的对角平均，波形立刻光滑了。

3.3 两个关键参数：嵌入维数L和分组个数

嵌入维数L可以这样理解：它决定了频率分辨率。如果两个源频率靠得特别近，比如48Hz和52Hz，数据长度只有1秒，那么需要足够大的L才能让它们在矩阵列空间里分开。一个粗略经验是，让L覆盖最低频率源一个周期的3到5倍以上。比如最低源频率20Hz，采样率1000Hz，一个周期50个点，L至少取150到250比较稳。

分组个数则是盲源分离最“玄”的部分。本例中我直接用了固定序号1:6和7:12，因为simulation里我知道只有两个源。实际工程中，推荐这样判断：

• 先看归一化奇异值曲线，找出明显拐点或“平台”之间的间隙；
• 周期成分通常会以成对奇异值出现，相邻两个奇异值大小接近，比如第1和第2个接近、第3和第4个接近；
• 如果数据里一个源含有多个谐波，这个源会占多组奇异值，需要用频谱或包络特征把属于同一源的分组合并。

我处理真实振动数据时，会先把前20个奇异值对应的V矩阵列做FFT，统计每个奇异向量主频是多少。主频相同的奇异向量大概率属于同一个源的不同谐波或共轭对，再按主频把它们聚成不同源。这个方法比单纯看奇异值曲线更稳定，尤其适用于带多个谐波的机械故障信号。

4. 一次完整的实操流程与注意事项

上面代码里的数据和参数都比较理想。这里以一个真实感很强的场景再走一遍全流程，方便你做替换和扩展。

4.1 处理自己的数据时，按这个顺序改

假设你有一段采样率fs_data、长度N_data的实测信号x_data，想分离其中两个频率相差较大的独立振动源。

第一步，去均值和去趋势。Hankel矩阵方法对直流分量和趋势特别敏感，慢变趋势会占掉前几个最大的奇异值，干扰后续分组。建议至少做一次去均值，如果信号有明显漂移，再用detrend或者多项式拟合扣除趋势。

第二步，选择L。建议先用一行代码快速估算目标频带：

lowFreq = 10; % 你关心的最低频率，单位Hz L = round(fs_data / lowFreq * 5); % 嵌入窗长相当于5个最低频周期 L = min(L, floor(N_data / 2)); % 别超过数据长度一半

注意不要让L太接近N，否则K变成几个点，SVD的行空间和列空间相关性会下降，反而分离不干净。

第三步，把x替换成去均值后的x_data，运行上面的SVD部分。然后重点看奇异谱图和V矩阵频谱。通常你会看到前几个奇异值明显偏大，再往后快速衰减。把明显偏大的部分对应的V矩阵列取出来，画在一起观察时域形状，判断哪些属于同一个源。

第四步，分组重构后，对每个重构分量做频谱验证。理想的分离结果应该是各分量频谱主峰清晰、互相不串扰。如果两个分量之间还有明显渗透，说明分界处选得不合适，或者两个源本身在频域上重叠严重，需要适当增大L。

4.2 多通道数据怎么扩展

单通道场景最苛刻，如果手里已经有多个传感器通道，可以把方法直接扩展成“多通道Hankel联合SVD”，效果会更好。思路如下：把每个通道分别构造Hankel矩阵，然后按行拼接成一个大矩阵D：

% X是M*N，M个通道，每行一个观测信号 % D是 (M*L) x K 的联合矩阵 D = []; for ch = 1:M xtmp = X(ch, :); Htmp = hankel(xtmp(1:L), xtmp(L:N)); D = [D; Htmp]; end [U, S, V] = svd(D, 'econ');

之后的分组和重构逻辑不变，只是重构H时需要把D中对应通道的块拆出来，再单独对角平均。多通道联合Hankel方法在机械故障诊断里非常实用，比如轴承座周围布置三个加速度传感器，每个通道测到的都是多个振源混合，但把通道联合起来后，不同振源的空间分布差异也被带进了矩阵结构，分离可靠性会大幅提高。

4.3 关于幅度和相位的固有不确定性

盲源分离本身有个数学上绕不开的问题：分离出来的源信号，幅值和相位不一定是真实的原值。基于Hankel矩阵SVD的方法同样如此，分离结果的波形形状、相对变化趋势可信，但可能会发生整体反相，也可能幅度和原源差一个比例常数。

这不算代码bug。比如我上次做某组数据时，第二个源分离出来波形非常完美，但整体翻了个面，和真实源完全反相。原因在于SVD奇异向量的符号本来就可以任意翻转，算法没有先验知识，不知道该选择正方向。所以对比效果时不要只看时域重叠，更多用相关系数的绝对值和频谱包络来判断。如果应用场景需要保留真实相位，得用已知标定信号做一次方向修正。

5. 实际问题排查：照这个表快速定位

我把实操中最高频的坑整理成了一张排查表，遇到问题直接对照：

关于分组，我再多提一句。别指望每种数据都有清晰断崖。噪声偏大的时候，奇异值曲线往往很平缓，这时单靠眼睛看前几阶特别容易出错。我的经验是，把SVD得到的U或V矩阵每一列画出来，观察它的振荡形式。真实的周期源会呈现规则振荡，而噪声列在时间上非常杂乱；同时还可以计算每一列的过零率或谱峰锐度，用这个量去辅助分类，比单纯看奇异值大小实用得多。

6. 从分离到应用：这套代码还能怎么扩展

如果只是分离两个仿真正弦，这个方法确实显得“大材小用”。但当你把它用到真实数据时，后续很多处理都得围绕分离结果展开。

比如旋转机械振动里，经常是正常转频、齿轮啮合频率、轴承故障频率混在一起。用Hankel盲源分离把转频分量剔除后，剩余分量做包络谱，故障特征会更突出。分离出来的低频慢变包络，本身就可以作为工况变化指标。在我的处理流程里，常把分离后的每个分量分别做时域统计和FFT，再把结果作为特征输入分类器，诊断效果往往优于直接对混合信号提特征。

另一个常见的扩展是数据清洗。有时我们不关心“每个源具体是什么”，只想要某一类源。比如心电信号里混着肌电干扰，但肌电不是严格周期信号，它会被平均到许多较小的奇异值上。重构时只保留前几个主奇异值，就相当于做了一个基于数据本身的自适应滤波。这个方法比普通带通滤波器更聪明，因为它不需要指定通带边界，而是让数据自动决定保留哪些时间结构。

我最近在做一个多传感器同步采集的项目时，也把这段核心逻辑推演成了另一个变体：把每个通道分别构造Hankel矩阵后，不直接拼接，而是对每个通道的Hankel矩阵先求出协方差矩阵，再做联合近似对角化。这种方式在通道数较多、各通道噪声不相干时，分离效果会更好。不过单从代码量来看，基础版本已经能解决大部分一周内就能复现的源分离需求了。

关于实际使用，我的体会是：先不要急着追求全自动选组，把每一阶奇异值对应的波形和频谱看懂，比什么都重要。盲源分离不是一次运行就出最终结论的魔法，它更像一个交互式工具箱——你观察一次结果，调整一次分界，再观察，再调整。等你积累了足够多的数据特征，就能慢慢摸索出一套适合自己应用场景的自动规则。再用这段MATLAB代码做底子，后面的路就顺了。