1. 从《三体》到现实为什么我们需要“熵”来监听宇宙想象一下你坐在一个巨大的射电望远镜控制室里面前是几十块屏幕上面滚动着永不停歇、看似毫无规律的波形曲线。你的工作就是盯着这些来自宇宙深处的“噪声”日复一日年复一年试图从中发现一丝不同寻常的“规律”——那可能是一个遥远文明发出的问候也可能是某种未知天体物理现象的信号。这听起来是不是很像《三体》里那位可怜的1379号监听员他凭借上千年的经验最终用肉眼察觉到了“波形似乎多了某种说不出来的东西”。但在现实中我们等不了上千年也几乎不可能训练出拥有如此直觉的科学家。我们需要一种更可靠、更自动化的方法来替我们“盯”着这些海量数据。这就是“信息熵”大显身手的地方。简单来说信息熵是衡量一段信号“混乱程度”或“不确定性”的数学工具。一个完全随机、没有任何规律的白噪声它的信息熵是最高的因为它的“不确定性”最大你永远猜不到下一个数据点会是什么值。相反一段被“智能调剂”过的信号比如一段有规律的无线电广播、一个重复的脉冲或者一个叠加了特定频率的正弦波它的内部是有结构和秩序的。这种秩序会降低信号的“不确定性”从而使其信息熵的值降低。所以监听宇宙信号的逻辑就变得异常清晰了我们不再需要肉眼去分辨波形那细微的“灵魂”而是编写一段程序让它持续计算接收到的信号的信息熵。当熵值长期稳定在一个较高的水平时我们知道那是宇宙背景噪声可以安心喝咖啡。一旦熵值突然、显著地下降程序就会立刻报警“注意这里可能有‘东西’” 这种方法不仅解放了人力而且更加客观、可量化、可重复是迈向自动化信号检测的关键一步。但直接对原始的时域信号计算信息熵效果往往不尽如人意。因为很多有用的信息隐藏在信号的频率结构、内在模式或能量分布中。这就引出了我们今天要深入探讨的三种强大的“熵武器”功率谱熵、奇异谱熵和能量熵。它们就像三副不同的“眼镜”能让我们从频域、矩阵分解和信号分解的角度更清晰地看到信号的本质从而更精准地从噪声的海洋中捞出我们想要的“珍珠”。接下来我就结合自己处理天文和通信信号的经验带你一步步理解这三种方法并用MATLAB手把手实现它们。2. 理解信息熵不只是数学公式在深入那三种具体的熵之前我们得先把信息熵这个基础概念吃透。很多教程一上来就扔公式看得人头大。咱们换个方式用个更生活的例子来理解。假设你是一家公司的老板每天下午都要猜你的员工小明为什么没来上班。日复一日你积累了经验80%的情况他是睡过头了15%是生病了还有5%是……被外星人抓走了纯粹举例。现在有人跑来告诉你“老板小明今天又没来” 这句话带给你的“信息量”有多大呢其实不大因为大概率就是睡过了你几乎能猜到。但如果有人告诉你“老板小明今天没来原因是被外星人抓走了” 这句话的信息量就爆炸了因为这是一个极小概率事件完全出乎你的意料。香农的信息熵本质上就是对这种“平均惊喜度”的量化。它计算的是当你得知小明没来上班的“原因”时平均能获得多少信息量。如果原因总是确定的比如100%睡过头那你就毫无惊喜信息熵为0。如果各种原因概率均等睡过头、生病、被外星人抓走各占三分之一那你每次猜测的悬念最大信息熵也达到最高值。把这个逻辑平移到信号上一段完全随机的白噪声每个样本点都像是一个独立且概率均等的事件你永远无法预测下一个点所以它的“平均惊喜度”很高信息熵很大。而一段包含固定频率正弦波的信号它的值在-1和1之间有规律地摆动可预测性很强你的“惊喜度”很低所以信息熵就小。因此信息熵越低信号越有序越可能包含我们感兴趣的“信息”。这就是我们用它来检测非随机信号的理论基石。2.1 从理论到感知一个简单的MATLAB实验光说不练假把式我们立刻用MATLAB来感受一下。我们来生成两段信号一段纯白噪声一段是白噪声叠加了一个明显的正弦波。% 生成示例信号 fs 1000; % 采样率 1000 Hz t 0:1/fs:1-1/fs; % 1秒时间向量 N length(t); % 信号长度 1000 % 生成高斯白噪声 rng(default); % 保证可重复性 white_noise 0.5 * randn(1, N); % 生成一个10Hz的正弦波模拟“智能信号” sine_wave 0.8 * sin(2*pi*10*t); % 合成信号纯噪声 和 噪声信号 signal_pure_noise white_noise; signal_with_info white_noise sine_wave; % 绘制时域图 figure; subplot(2,1,1); plot(t, signal_pure_noise); title(纯白噪声信号); xlabel(时间 (s)); ylabel(幅值); grid on; subplot(2,1,2); plot(t, signal_with_info); title(包含10Hz正弦波的信号); xlabel(时间 (s)); ylabel(幅值); grid on;运行这段代码你会看到两张图。坦白说仅凭肉眼你能非常肯定地区分哪一段包含了规律性的正弦波吗可能第二段信号的波形看起来“胖”了一点规律性隐约可见但并不绝对明显尤其是当噪声强度变化时。我们的眼睛和大脑对时域波形中的周期性并不总是那么敏感。接下来我们就请出第一位“熵侦探”——功率谱熵让它来告诉我们答案。3. 功率谱熵在频率的王国里寻找秩序功率谱熵是我在分析周期性或谐振类信号时最常用的工具。它的核心思想是不去直接看信号随时间怎么变而是看它的能量在不同频率上是如何分布的。一个纯噪声它的能量会平均地分布在所有频率上像一片平坦的沙滩。而一个有规律的信号比如我们的10Hz正弦波它的能量会高度集中在10Hz这个频率点附近就像沙滩上隆起的一座沙堡。功率谱熵的计算就是先画出这片“能量沙滩”的等高线图即功率谱然后计算这些能量分布的“不均匀程度”。能量分布越集中沙堡越高越尖不确定性越低熵值就越小能量分布越分散沙滩越平坦不确定性越高熵值就越大。3.1 手把手计算功率谱熵我们来一步步拆解计算过程并附上详细的MATLAB代码和注释。第一步计算信号的功率谱密度PSD。这是将信号从时域变换到频域的关键步骤。MATLAB里可以用pwelch函数它使用Welch方法能有效减少频谱估计的方差结果更平滑稳定。function [pxx, f] compute_psd(signal, fs) % 计算信号的功率谱密度PSD % 输入 % signal: 输入时域信号 % fs: 采样频率 % 输出 % pxx: 功率谱密度估计值 % f: 对应的频率向量 window min(256, floor(length(signal)/4)); % 窗长取信号长度的1/4或256点中的较小值 noverlap floor(window/2); % 重叠率通常为窗长的一半 nfft max(256, 2^nextpow2(window)); % FFT点数保证足够分辨率 [pxx, f] pwelch(signal, window, noverlap, nfft, fs); end第二步将功率谱归一化得到概率分布。我们把每个频率点上的功率值看作是能量分布在该频率上的“概率”。需要将所有功率值加起来然后计算每个值占总和的比例。function prob_dist normalize_psd(pxx) % 将功率谱密度归一化为概率分布 % 输入pxx - 功率谱密度向量 % 输出prob_dist - 概率分布向量所有元素和为1 total_power sum(pxx); if total_power 0 error(信号总功率为零无法归一化。); end prob_dist pxx / total_power; end第三步将概率分布代入信息熵公式计算。这就是标准的香农熵计算了。注意我们通常约定0 * log2(0) 0所以代码中需要处理概率为零的情况。function entropy_val shannon_entropy(prob_dist) % 计算香农熵 % 输入prob_dist - 概率分布向量所有元素和为1 % 输出entropy_val - 计算得到的信息熵值单位比特 % 移除概率为零的元素避免log2(0)产生-Inf prob_dist prob_dist(prob_dist 0); entropy_val -sum(prob_dist .* log2(prob_dist)); end第四步整合成功率谱熵计算函数。现在我们把前三步组合起来形成一个完整的功率谱熵计算函数。function psd_entropy compute_psd_entropy(signal, fs) % 计算信号的功率谱熵 % 输入 % signal: 输入时域信号行向量或列向量 % fs: 采样频率 % 输出 % psd_entropy: 功率谱熵值 % 1. 计算功率谱密度 [pxx, ~] compute_psd(signal, fs); % 2. 归一化为概率分布 prob normalize_psd(pxx); % 3. 计算香农熵 psd_entropy shannon_entropy(prob); end3.2 实战检验找出隐藏的正弦波让我们用这个函数来检验之前生成的两段信号。% 计算两段信号的功率谱熵 psdE_noise compute_psd_entropy(signal_pure_noise, fs); psdE_with_sine compute_psd_entropy(signal_with_info, fs); fprintf(纯噪声信号的功率谱熵: %.4f\n, psdE_noise); fprintf(含正弦波信号的功率谱熵: %.4f\n, psdE_with_sine); % 可视化功率谱 [pxx1, f1] compute_psd(signal_pure_noise, fs); [pxx2, f2] compute_psd(signal_with_info, fs); figure; subplot(2,1,1); plot(f1, 10*log10(pxx1)); % 用分贝(dB)显示更直观 title([纯噪声功率谱 (PSD Entropy , num2str(psdE_noise, %.2f), )]); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on; xlim([0, fs/2]); subplot(2,1,2); plot(f2, 10*log10(pxx2)); title([含正弦波信号功率谱 (PSD Entropy , num2str(psdE_with_sine, %.2f), )]); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on; xlim([0, fs/2]); hold on; % 标记出10Hz的位置 [~, idx] min(abs(f2-10)); plot(f2(idx), 10*log10(pxx2(idx)), r*, MarkerSize, 10); text(f2(idx)5, 10*log10(pxx2(idx)), 10Hz Peak, Color, red);运行这段代码你会看到两个鲜明的对比。纯噪声的功率谱像一条起伏的“毛毯”在各个频率上都有能量因此计算出的功率谱熵值较高例如在4.6左右。而包含了正弦波的信号其功率谱在10Hz处出现了一个尖锐的峰值大部分能量都集中于此其他频率的能量相对很少因此其功率谱熵值会显著降低例如在3.2左右。这个差异是如此明显以至于你可以轻松地设定一个阈值比如4.0一旦熵值低于它就触发警报。在实际的宇宙信号监听中你可以对持续流入的数据流进行滑动窗计算实时监控熵值的变化从而实现自动化检测。4. 奇异谱熵挖掘信号内部的空间结构功率谱熵在频域上很强大但对于一些特征不明显体现在频率上的信号或者数据量很少、噪声很强的短序列我们就需要另一种工具——奇异谱熵。它的思路非常巧妙把一维的时间序列重新组织成一个矩阵然后通过分析这个矩阵的“骨架”奇异值来揭示信号内在的复杂性。你可以把信号想象成一串珍珠项链。奇异谱分析SSA的工作是1把这串项链按一定规则盘绕成一个二维的“珍珠矩阵”2对这个矩阵进行一种叫做“奇异值分解”SVD的数学操作这相当于把矩阵拆解成几个具有不同“重要性”的组成部分3这些组成部分的“重要性”大小就是奇异值。如果信号主要是由少数几个强大的规律比如一个正弦波和一个趋势项构成的那么只有前几个奇异值会很大其他的都很小。如果信号是完全随机的噪声那么奇异值的大小会比较平均。奇异谱熵就是计算这些奇异值分布的均匀程度。分布越不均匀少数奇异值主导熵值越低说明信号结构越简单、越有规律分布越均匀熵值越高说明信号越接近随机噪声。4.1 奇异谱熵的计算步骤与代码实现第一步相空间重构构建轨迹矩阵。这是将一维信号升维到二维的关键。我们设定一个窗口长度L通常取信号长度的1/3到1/2然后像滑动窗口一样截取信号每一列就是窗口内的一段信号。function trajectory_matrix build_trajectory_matrix(signal, L) % 构建轨迹矩阵相空间重构 % 输入 % signal: 一维时间序列 (1 x N) % L: 窗口长度嵌入维度 % 输出 % trajectory_matrix: 轨迹矩阵 (L x K) N length(signal); K N - L 1; trajectory_matrix zeros(L, K); for i 1:K trajectory_matrix(:, i) signal(i:iL-1); end end第二步对轨迹矩阵进行奇异值分解SVD。使用MATLAB内置的svd函数我们可以轻松得到奇异值。奇异值通常按从大到小排列。function singular_values compute_svd_of_trajectory(trajectory_matrix) % 对轨迹矩阵进行奇异值分解并返回奇异值 % 输入trajectory_matrix - 轨迹矩阵 % 输出singular_values - 奇异值向量降序排列 [~, S, ~] svd(trajectory_matrix, econ); singular_values diag(S); end第三步计算奇异谱熵。将奇异值归一化为概率分布每个奇异值占所有奇异值总和的比例然后代入香农熵公式。function svd_entropy compute_svd_entropy(signal, L) % 计算信号的奇异谱熵 % 输入 % signal: 输入时域信号 % L: 窗口长度嵌入维度如果未指定则使用启发式规则 % 输出 % svd_entropy: 奇异谱熵值 if nargin 2 || isempty(L) % 一个常用的启发式规则L取信号长度的1/3 L floor(length(signal) / 3); L max(L, 2); % 确保L至少为2 end % 1. 构建轨迹矩阵 X build_trajectory_matrix(signal, L); % 2. 奇异值分解 s compute_svd_of_trajectory(X); % 3. 归一化奇异值为概率分布 s_total sum(s); if s_total 0 svd_entropy 0; return; end p s / s_total; % 4. 计算香农熵 svd_entropy shannon_entropy(p); end4.2 窗口长度L的选择一个实践中的小坑这里有个非常重要的细节窗口长度L的选择会显著影响奇异谱熵的结果。L太小重构的矩阵无法捕捉信号的长程相关性L太大矩阵会变得冗余且计算量增加对噪声更敏感。在我的项目中我通常遵循以下原则经验法则L取信号长度N的1/3到1/2是一个不错的起点。多次尝试对于关键应用我会计算L在不同取值下的熵值观察其稳定性。如果在一个范围内熵值变化平缓则说明结果相对可靠。领域知识如果你对信号潜在的周期有先验知识可以将L设为该周期的整数倍这样有助于分离出周期分量。让我们用奇异谱熵来测试之前的信号并与功率谱熵对比。% 计算奇异谱熵 L 50; % 尝试不同的L值例如50, 100, 150观察结果 svdE_noise compute_svd_entropy(signal_pure_noise, L); svdE_with_sine compute_svd_entropy(signal_with_info, L); fprintf(窗口长度 L %d\n, L); fprintf(纯噪声信号的奇异谱熵: %.4f\n, svdE_noise); fprintf(含正弦波信号的奇异谱熵: %.4f\n, svdE_with_sine); % 可视化奇异值分布 figure; X_noise build_trajectory_matrix(signal_pure_noise, L); s_noise compute_svd_of_trajectory(X_noise); p_noise s_noise / sum(s_noise); X_sine build_trajectory_matrix(signal_with_info, L); s_sine compute_svd_of_trajectory(X_sine); p_sine s_sine / sum(s_sine); subplot(2,1,1); stem(p_noise, filled); title([纯噪声信号奇异值分布 (SVD Entropy , num2str(svdE_noise, %.2f), )]); xlabel(奇异值序号); ylabel(归一化奇异值概率); grid on; subplot(2,1,2); stem(p_sine, filled); title([含正弦波信号奇异值分布 (SVD Entropy , num2str(svdE_with_sine, %.2f), )]); xlabel(奇异值序号); ylabel(归一化奇异值概率); grid on;你会发现对于包含正弦波的信号前几个奇异值的概率高度会远大于后面的分布很不均匀因此熵值较低。而对于白噪声奇异值的分布则相对平缓熵值较高。奇异谱熵特别擅长处理非平稳、短数据、含噪声的信号在天文信号处理中当我们接收到的数据段很短或者信号特性随时间变化时它往往能提供比功率谱熵更稳健的特征。5. 能量熵融合现代信号分解技术的灵活武器能量熵是一个更广义、更灵活的概念。它的核心思想是先用某种方法把原始信号“拆开”成若干个分量然后看能量在这些分量之间是如何分配的最后计算这种分配的不确定性。这个“拆开”的方法可以是经验模态分解EMD、变分模态分解VMD、小波分解Wavelet等等。能量熵的强大之处在于它能继承所用分解方法的优点。例如EMD能自适应地分解出信号的本征模态函数IMF每个IMF代表一个特定时间尺度的振荡。如果某个信号状态比如出现了一个瞬态脉冲导致能量在某个IMF上突然集中那么能量熵就会下降。小波分解则同时具有时域和频域定位能力能量熵可以反映信号能量在不同时频“格子”里的集中程度。5.1 以小波能量熵为例的实战我们以小波分解为例因为它有成熟且易用的工具箱。我们将使用wavedec函数进行多级小波分解。第一步选择小波基和分解层数进行分解。小波基如db4和层数需要根据信号特点选择。层数决定了分解的精细程度。function [c, l] wavelet_decomposition(signal, wavelet_name, level) % 使用小波进行多尺度分解 % 输入 % signal: 输入信号 % wavelet_name: 小波基名称如 db4, sym5 % level: 分解层数 % 输出 % c: 分解系数向量 % l: 记录各层系数长度的向量 [c, l] wavedec(signal, level, wavelet_name); end第二步计算各层细节系数和近似系数的能量。小波分解后我们得到各层的高频细节系数D1, D2, ...和最后一层的低频近似系数A。能量通常定义为系数的平方和。function energy_vec compute_wavelet_energy(c, l, level) % 计算小波分解后各分量的能量 % 输入 % c, l: wavedec函数的输出 % level: 分解层数 % 输出 % energy_vec: 能量向量 [E_A, E_Dlevel, ..., E_D1] energy_vec zeros(1, level 1); % 提取近似系数 (A) 的能量 approx_coef appcoef(c, l, db4, level); % 注意这里需要指定小波基 energy_vec(1) sum(approx_coef.^2); % 提取各层细节系数 (D) 的能量 for i 1:level detail_coef detcoef(c, l, i); energy_vec(i 1) sum(detail_coef.^2); end end第三步计算能量熵。将各分量能量归一化为概率分布计算香农熵。function we_entropy compute_wavelet_energy_entropy(signal, wavelet_name, level) % 计算信号的小波能量熵 % 输入 % signal: 输入信号 % wavelet_name: 小波基名称 % level: 分解层数 % 输出 % we_entropy: 小波能量熵值 % 1. 小波分解 [c, l] wavelet_decomposition(signal, wavelet_name, level); % 2. 计算各分量能量 energy compute_wavelet_energy(c, l, level); % 3. 归一化为概率 total_energy sum(energy); if total_energy 0 we_entropy 0; return; end p energy / total_energy; % 4. 计算香农熵 we_entropy shannon_entropy(p); end5.2 测试与对比能量熵的表现% 计算小波能量熵 wavelet_name db4; level 5; % 分解5层 weE_noise compute_wavelet_energy_entropy(signal_pure_noise, wavelet_name, level); weE_with_sine compute_wavelet_energy_entropy(signal_with_info, wavelet_name, level); fprintf(小波基: %s, 分解层数: %d\n, wavelet_name, level); fprintf(纯噪声信号的小波能量熵: %.4f\n, weE_noise); fprintf(含正弦波信号的小波能量熵: %.4f\n, weE_with_sine); % 可视化能量分布 [~, l_noise] wavelet_decomposition(signal_pure_noise, wavelet_name, level); energy_noise compute_wavelet_energy(c_noise, l_noise, level); p_noise energy_noise / sum(energy_noise); [~, l_sine] wavelet_decomposition(signal_with_info, wavelet_name, level); energy_sine compute_wavelet_energy(c_sine, l_sine, level); p_sine energy_sine / sum(energy_sine); figure; subplot(2,1,1); bar(p_noise); title([纯噪声信号小波能量分布 (Wavelet Energy Entropy , num2str(weE_noise, %.2f), )]); xlabel(分量 (A5, D5, D4, D3, D2, D1)); ylabel(能量占比); grid on; xticklabels({A5, D5, D4, D3, D2, D1}); subplot(2,1,2); bar(p_sine); title([含正弦波信号小波能量分布 (Wavelet Energy Entropy , num2str(weE_with_sine, %.2f), )]); xlabel(分量 (A5, D5, D4, D3, D2, D1)); ylabel(能量占比); grid on; xticklabels({A5, D5, D4, D3, D2, D1});对于我们的例子10Hz的正弦波能量可能会主要集中在小波分解的某个特定尺度比如D3或D4层导致该层的能量占比异常高从而使能量分布变得不均匀熵值降低。而白噪声的能量则会相对均匀地分布在各个尺度上熵值较高。能量熵的灵活性在于你可以根据信号特点选择最合适的分解方法。比如处理非平稳信号用EMD追求数学严谨性用VMD需要时频分析则用小波。它更像一个框架让你能融合最先进的信号处理工具来提取特征。6. 综合实战构建一个简易的宇宙信号监听原型理论和方法都清楚了现在我们来模拟一个真实的场景假设我们有一个实时数据流模拟从射电望远镜接收到的宇宙信号。数据流中大部分时间是背景噪声但偶尔会嵌入一段微弱的、可能是人工起源的信号。我们的任务是编写一个实时检测程序当“可疑信号”出现时发出警报。我们将采用滑动窗口的方法并同时计算功率谱熵和奇异谱熵作为双保险。能量熵计算稍慢这里暂不用于实时检测但可用于事后对报警片段的深入分析。6.1 系统设计思路与参数设置数据流模拟我们生成一段长时间序列在其中随机位置插入几段“智能信号”例如频率漂移的正弦波、线性调频信号等。滑动窗口定义一个固定长度的窗口例如包含1000个数据点以一定的步长例如100点在数据流上滑动。特征提取对每个窗口内的数据计算其功率谱熵和奇异谱熵。阈值判断根据历史噪声窗口的熵值统计确定一个报警阈值例如均值减去3倍标准差。当某个窗口的熵值低于阈值时触发报警。报警与记录记录报警时间、窗口数据和对应的熵值。6.2 MATLAB实现代码%% 模拟宇宙信号数据流 fs 1000; % 采样率 1kHz total_time 60; % 总时长60秒 t_total 0:1/fs:total_time-1/fs; N_total length(t_total); % 生成基础背景噪声可模拟不同强度的噪声 background_noise 1.0 * randn(1, N_total); % 在随机位置插入“智能信号” num_anomalies 5; anomaly_duration 0.5; % 每次异常信号持续0.5秒 anomaly_samples round(anomaly_duration * fs); signal_stream background_noise; % 初始化数据流为纯噪声 for i 1:num_anomalies % 随机选择插入起始点避免在开头和结尾 start_idx randi([fs*10, N_total - anomaly_samples - fs*10]); end_idx start_idx anomaly_samples - 1; % 生成一种“智能信号”这里用频率缓变的chirp信号为例 t_anomaly t_total(start_idx:end_idx); % 生成一个频率从5Hz线性变化到15Hz的线性调频信号 f0 5; f1 15; chirp_signal 0.7 * chirp(t_anomaly, f0, anomaly_duration, f1); % 将异常信号叠加到背景噪声上 signal_stream(start_idx:end_idx) signal_stream(start_idx:end_idx) chirp_signal; fprintf(植入异常信号段 %d: 时间 %.2fs 到 %.2fs\n, i, t_total(start_idx), t_total(end_idx)); end %% 实时监听系统参数设置 window_size 1000; % 滑动窗口大小1000个点1秒数据 step_size 100; % 滑动步长100个点0.1秒 num_windows floor((N_total - window_size) / step_size) 1; % 初始化存储数组 psd_entropy_vals zeros(1, num_windows); svd_entropy_vals zeros(1, num_windows); time_stamps zeros(1, num_windows); % 每个窗口中心点对应的时间 % 奇异谱熵窗口长度嵌入维度 L_svd 200; %% 滑动窗口计算熵值 fprintf(开始滑动窗口分析...\n); for win_idx 1:num_windows start_sample (win_idx - 1) * step_size 1; end_sample start_sample window_size - 1; current_window signal_stream(start_sample:end_sample); current_time t_total(round((start_sample end_sample)/2)); time_stamps(win_idx) current_time; % 计算功率谱熵 psd_entropy_vals(win_idx) compute_psd_entropy(current_window, fs); % 计算奇异谱熵 svd_entropy_vals(win_idx) compute_svd_entropy(current_window, L_svd); % 每处理100个窗口打印一次进度 if mod(win_idx, 100) 0 fprintf(已处理 %d/%d 个窗口...\n, win_idx, num_windows); end end fprintf(分析完成\n); %% 确定报警阈值基于前20秒的“纯噪声”基线 baseline_duration 20; % 前20秒假设为纯噪声 baseline_samples baseline_duration * fs; baseline_window_idx find(time_stamps baseline_duration); psd_baseline psd_entropy_vals(baseline_window_idx); svd_baseline svd_entropy_vals(baseline_window_idx); % 阈值设定均值 - 3倍标准差假设熵值降低为异常 psd_threshold mean(psd_baseline) - 3 * std(psd_baseline); svd_threshold mean(svd_baseline) - 3 * std(svd_baseline); fprintf(功率谱熵报警阈值: %.4f\n, psd_threshold); fprintf(奇异谱熵报警阈值: %.4f\n, svd_threshold); %% 检测报警 alarm_psd psd_entropy_vals psd_threshold; alarm_svd svd_entropy_vals svd_threshold; alarm_combined alarm_psd alarm_svd; % 双指标同时报警更可靠 alarm_times time_stamps(alarm_combined); if ~isempty(alarm_times) fprintf(检测到异常报警报警时间点秒: \n); disp(alarm_times); else fprintf(未检测到异常。\n); end %% 可视化结果 figure(Position, [100, 100, 1200, 800]); % 子图1原始数据流标记异常植入区域 subplot(4,1,1); plot(t_total, signal_stream); title(模拟宇宙信号数据流含植入的“智能信号”); xlabel(时间 (s)); ylabel(幅值); grid on; hold on; % 高亮显示植入信号的区域仅用于验证实际中未知 for i 1:num_anomalies start_idx find(t_total (i*10 5), 1); % 简化假设在之前植入循环中记录了位置 % 实际应用中应使用之前植入时记录的位置这里为演示简化 ylims ylim; fill([t_total(start_idx), t_total(start_idxanomaly_samples), t_total(start_idxanomaly_samples), t_total(start_idx)], ... [ylims(1), ylims(1), ylims(2), ylims(2)], y, FaceAlpha, 0.3, EdgeColor, none); end legend(数据流, 植入信号区域仅演示); % 子图2功率谱熵随时间变化 subplot(4,1,2); plot(time_stamps, psd_entropy_vals, b-, LineWidth, 1); hold on; yline(psd_threshold, r--, LineWidth, 1.5, Label, PSD熵报警阈值); title(滑动窗口功率谱熵); xlabel(时间 (s)); ylabel(功率谱熵值); grid on; legend(功率谱熵, 阈值); % 子图3奇异谱熵随时间变化 subplot(4,1,3); plot(time_stamps, svd_entropy_vals, g-, LineWidth, 1); hold on; yline(svd_threshold, r--, LineWidth, 1.5, Label, SVD熵报警阈值); title(滑动窗口奇异谱熵); xlabel(时间 (s)); ylabel(奇异谱熵值); grid on; legend(奇异谱熵, 阈值); % 子图4综合报警信号 subplot(4,1,4); stem(time_stamps, alarm_combined, r^, MarkerSize, 8, LineWidth, 1.5); title(综合异常报警PSD熵 SVD熵均低于阈值); xlabel(时间 (s)); ylabel(报警 (1是)); ylim([-0.1, 1.5]); grid on;运行这个模拟系统你会看到随着滑动窗口的移动两个熵值曲线在大部分时间保持高位波动对应背景噪声。而在我们植入“智能信号”的时间段两条熵值曲线会几乎同步地出现明显的“凹陷”并跌破阈值线从而触发报警。这个原型系统清晰地展示了如何将信息熵理论应用于一个近乎真实的自动化监听场景。在实际部署中你还需要考虑更复杂的噪声模型、阈值自适应更新、报警聚合与验证等环节但核心的骨架已经在这里了。