企业官网建设流程全解析

简介：这份资源是面向统计信号处理学习者与科研人员的广义最大似然比检验（GLRT）MATLAB仿真实现，聚焦弱信号检测与噪声环境下的异常判别问题，适合已具备假设检验基础、希望动手复现检测算法的高年级本科生或研究生。压缩包共7个文件，全部为m脚本，整体约3KB，涵盖数据生成、似然函数计算、临界值确定、决策规则与性能评估等环节，可配合连续波信号模型分析检测性能。已有452人学习下载，说明该方向具备一定关注度。读者可借助脚本理解似然比统计量的构造逻辑，观察不同显著性水平与信噪比下虚警率与检测概率的变化，并掌握蒙特卡洛仿真评估检测器性能的完整流程，为雷达、通信等场景中的信号检测研究提供可复用的代码框架与排错参考。

1. 从一段“跑出来全是1”的检测代码说起：GLRT 与似然比检验到底在解决什么

如果你写过信号检测或者异常检测的代码，大概率遇到过这种场景：仿真数据里明明有目标，检测统计量画出来却和噪声分不开，判决门限怎么调都是虚警。这时候很多人第一反应是换模型、加特征，但真正的问题往往出在检验统计量本身没构造对。MLE.rar_GLRT_Statistical test_spellcw5_似然比_似然比检验这个标题拆开看，核心就是三件事：广义似然比检验（GLRT）、统计检验框架、以及似然比这个统计量怎么算。它解决的是在未知参数存在时，如何构造一个仍然具备最优检测性能的判决规则。适合做雷达、通信、声呐、工业异常检测的工程师，也适合任何需要从含噪观测里判断“有没有信号”的人。这一章先把 GLRT 和似然比检验的关系讲清楚，后面再落到代码和参数上。

2. 似然比检验的统计底子：从 Neyman-Pearson 到 GLRT 的推导链

2.1 简单假设下的似然比检验为什么是最优的

统计检验的基本问题是：观测到一组数据 $x$，要判断它来自 $H_0$（只有噪声）还是 $H_1$（信号加噪声）。Neyman-Pearson 引理给出了一个非常干净的结论——在给定虚警概率 $P_{fa}$ 的上限约束下，使检测概率 $P_d$ 最大的判决规则，一定可以写成似然比的形式：

$$ \Lambda(x) = \frac{p(x|H_1)}{p(x|H_0)} \underset{H_0}{\overset{H_1}{\gtrless}} \eta $$

其中 $\eta$ 是由 $P_{fa}$ 确定的门限。这个结论不依赖任何具体的分布形式，只要两个假设下的概率密度函数已知，它就是最优的。实际工程里，我们通常对 $\Lambda(x)$ 取对数，得到对数似然比 $\ln\Lambda(x)$，这样可以把乘积变成求和，数值上更稳定。

但这里有个前提：$p(x|H_1)$ 和 $p(x|H_0)$ 必须是完全已知的。现实中这个前提经常不成立——信号幅度未知、相位未知、噪声功率未知，这些未知参数让似然函数没法直接写出来。

2.2 未知参数出现后 GLRT 是怎么接手的

当假设下的概率密度含有未知参数时，Neyman-Pearson 引理不能直接用了。GLRT 的思路很直接：既然参数未知，那就先用最大似然估计（MLE）把参数估出来，再把估计值代回似然函数，构造广义似然比：

$$ \Lambda_G(x) = \frac{\max_{\theta_1} p(x|\theta_1, H_1)}{\max_{\theta_0} p(x|\theta_0, H_0)} $$

分子是在 $H_1$ 下对未知参数 $\theta_1$ 做最大似然，分母是在 $H_0$ 下对 $\theta_0$ 做最大似然。这个统计量不再保证最优，但在大样本下具有渐进最优性，而且实际表现通常足够好。标题里的MLE指的就是这一步——最大似然估计是 GLRT 的内核。

提示：GLRT 不是万能的。当未知参数维度很高、样本量又不够时，MLE 的估计误差会直接恶化检测性能。这时候要考虑降维或者引入先验。

2.3 从统计量到判决：门限怎么定

构造出 $\Lambda_G(x)$ 之后，判决规则是 $\Lambda_G(x) > \eta$ 则判 $H_1$。门限 $\eta$ 的确定方式取决于你用什么准则：

工程上最常用的是 Neyman-Pearson 准则，因为虚警率通常有明确的指标要求。但 $\Lambda_G$ 的分布往往没有闭式解，门限得靠蒙特卡洛仿真或者渐近分布来定。

3. 用 Python 把 GLRT 跑通：从数据生成到判决的完整链路

3.1 生成仿真数据并定义假设模型

先搭一个最小可复现的仿真环境。假设接收信号模型为 $x[n] = A s[n] + w[n]$，$H_0$ 时 $A=0$，$H_1$ 时 $A>0$ 但具体值未知，$w[n]$ 是高斯白噪声，功率 $\sigma^2$ 也未知。这就是经典的未知幅度、未知噪声功率场景。

import numpy as np def generate_data(N, A_true, sigma_true, seed=42): """ 生成 GLRT 仿真数据 N: 样本数 A_true: 真实信号幅度，H0 时为 0 sigma_true: 真实噪声标准差 """ rng = np.random.default_rng(seed) s = np.ones(N) # 确定性信号，实际可替换为任意已知波形 w = rng.normal(0, sigma_true, N) x = A_true * s + w return x, s # 生成 H0 和 H1 各一组数据 x0, s = generate_data(N=64, A_true=0.0, sigma_true=1.0, seed=1) x1, _ = generate_data(N=64, A_true=0.8, sigma_true=1.0, seed=2)

这段代码定义了信号模型和噪声模型。s是已知波形，这里用全 1 向量简化，实际中可以是匹配滤波器的参考信号。A_true和sigma_true是仿真用的真值，GLRT 算法本身不知道它们。

3.2 推导并实现 GLRT 统计量

对于上面的模型，$H_1$ 下 $x \sim \mathcal{N}(A s, \sigma^2 I)$，$H_0$ 下 $x \sim \mathcal{N}(0, \sigma^2 I)$。未知参数是 $A$ 和 $\sigma^2$。在 $H_1$ 下对 $A$ 和 $\sigma^2$ 做 MLE：

$$ \hat{A} = \frac{s^T x}{s^T s}, \quad \hat{\sigma}_1^2 = \frac{1}{N}(x - \hat{A}s)^T(x - \hat{A}s) $$

在 $H_0$ 下 $\sigma^2$ 的 MLE 是 $\hat{\sigma}_0^2 = \frac{1}{N}x^T x$。代入广义似然比并取对数，化简后得到：

$$ \ln\Lambda_G(x) = \frac{N}{2}\ln\frac{\hat{\sigma}_0^2}{\hat{\sigma}_1^2} $$

def glrt_statistic(x, s): """ 计算 GLRT 统计量 x: 观测向量 s: 已知信号波形 """ N = len(x) # H1 下的 MLE A_hat = np.dot(s, x) / np.dot(s, s) residual_h1 = x - A_hat * s sigma1_sq = np.dot(residual_h1, residual_h1) / N # H0 下的 MLE sigma0_sq = np.dot(x, x) / N # 对数广义似然比 stat = (N / 2) * np.log(sigma0_sq / sigma1_sq) return stat stat_h0 = glrt_statistic(x0, s) stat_h1 = glrt_statistic(x1, s) print(f"H0 统计量: {stat_h0:.4f}") print(f"H1 统计量: {stat_h1:.4f}")

这里的关键是sigma0_sq / sigma1_sq这个比值。当信号存在时，$H_1$ 下的残差功率明显小于 $H_0$ 下的总功率，比值大于 1，对数值为正。信号越强，统计量越大。N/2这个系数来自高斯分布的对数似然推导，实际使用时可以吸收到门限里，不影响检测性能。

3.3 用蒙特卡洛仿真确定门限

门限不能拍脑袋定。标准做法是在 $H_0$ 下跑大量仿真，得到统计量的经验分布，然后取对应 $P_{fa}$ 的分位数。

def estimate_threshold(N, sigma, pfa, n_trials=100000, seed=0): """ 蒙特卡洛估计门限 N: 样本数 sigma: 噪声标准差 pfa: 目标虚警概率 n_trials: 仿真次数 """ rng = np.random.default_rng(seed) s = np.ones(N) stats = np.zeros(n_trials) for i in range(n_trials): x = rng.normal(0, sigma, N) # H0 下只有噪声 stats[i] = glrt_statistic(x, s) threshold = np.quantile(stats, 1 - pfa) return threshold thr = estimate_threshold(N=64, sigma=1.0, pfa=0.01) print(f"Pfa=0.01 对应的门限: {thr:.4f}")

n_trials取 100000 是为了让 0.01 分位数估计足够稳定。如果目标 $P_{fa}$ 更低，比如 $10^{-4}$，仿真次数至少要 100 万以上，否则分位数估计的方差会很大。这是很多人在做低虚警率检测时容易翻车的地方——仿真次数不够，门限偏了，实测虚警率完全失控。

3.4 检测性能评估与 ROC 曲线

有了门限之后，在 $H_1$ 下跑仿真统计检测概率，再扫一遍门限画 ROC 曲线。

def evaluate_detection(N, A, sigma, threshold, n_trials=10000, seed=10): """评估给定门限下的检测概率""" rng = np.random.default_rng(seed) s = np.ones(N) detections = 0 for _ in range(n_trials): x = A * s + rng.normal(0, sigma, N) if glrt_statistic(x, s) > threshold: detections += 1 return detections / n_trials pd = evaluate_detection(N=64, A=0.8, sigma=1.0, threshold=thr) print(f"Pfa=0.01, A=0.8 时 Pd={pd:.4f}")

参数说明：A是信号幅度，sigma是噪声标准差，两者之比就是信噪比。A=0.8, sigma=1.0对应 SNR 约 $-1.9$ dB（按 $20\log_{10}(A/\sigma)$ 算）。这个信噪比下 GLRT 还能保持不错的检测概率，说明统计检验框架在低信噪比下确实有效。

4. 避坑与排查：GLRT 落地时最容易翻车的五个地方

4.1 统计量算出来全是 NaN 或者 Inf

现象：跑仿真时glrt_statistic返回nan或inf，后续判决全部失效。

原因：sigma1_sq或sigma0_sq出现零或负数。常见于信号波形s全零、观测x全零、或者数值精度导致残差平方和算成负数。

解决：在计算前加保护，sigma0_sq = max(np.dot(x, x) / N, 1e-12)，sigma1_sq同理。另外检查s是否归一化，np.dot(s, s)为零会导致A_hat除零。

4.2 门限仿真次数不够导致虚警率失控

现象：仿真设定 $P_{fa}=0.01$，实测虚警率却是 0.05 甚至更高。

原因：蒙特卡洛次数太少，高分位数估计不准。100000 次仿真对 0.01 分位数的标准差大约是 $\sqrt{0.01 \times 0.99 / 100000} \approx 0.0003$，看起来还行，但如果统计量分布尾部较厚，实际误差会更大。

解决：目标 $P_{fa}$ 每降低一个数量级，仿真次数至少增加 10 倍。或者改用渐近分布——GLRT 统计量在 $H_0$ 下通常近似服从卡方分布，可以用解析式定门限，再用少量仿真验证。

4.3 未知参数维度太高导致 MLE 不稳定

现象：检测性能远低于理论预期，ROC 曲线几乎贴着对角线。

原因：未知参数太多，MLE 估计方差大，代回似然函数后统计量的区分度被稀释。比如同时估计幅度、相位、频率、噪声功率，四个参数一起估，样本量又只有几十个，MLE 本身就不靠谱。

解决：能确定参数就确定，能降维就降维。相位未知可以用非相干积累，频率未知可以先做粗估计再代入。实在要估，增加样本量或者引入正则化。

4.4 把 GLRT 当成最优检验用

现象：明明用了 GLRT，性能却不如一个简单的能量检测器。

原因：GLRT 只是渐进最优，小样本下不一定比简单方法好。能量检测器在未知幅度、已知噪声功率的场景下就是最优的，GLRT 反而因为要估噪声功率而损失性能。

解决：先想清楚哪些参数真的未知。如果噪声功率已知，就不要估它，直接用已知值构造似然比。GLRT 的“广义”是有代价的。

4.5 信号波形与实际不匹配

现象：仿真里性能很好，换到实测数据就崩了。

原因：GLRT 里的s是假设已知的参考信号，实际中信号波形可能有畸变、时延、多普勒。波形失配会让 MLE 估出的幅度偏小，统计量整体下移。

解决：用匹配滤波的思路，对可能的时延和多普勒做搜索，取最大统计量作为最终判决量。或者改用广义匹配滤波器组，覆盖参数不确定范围。

5. 进阶技巧：用渐近分布替代蒙特卡洛定门限

蒙特卡洛定门限虽然通用，但低虚警率场景下计算量太大。一个更高效的做法是利用 GLRT 统计量的渐近分布。在正则条件下，$2\ln\Lambda_G$ 在 $H_0$ 下近似服从自由度为“$H_1$ 比 $H_0$ 多出的自由参数个数”的卡方分布。上面例子中，$H_1$ 比 $H_0$ 多了一个未知幅度 $A$，所以自由度是 1。

from scipy.stats import chi2 def threshold_from_chi2(pfa, df=1): """用卡方分布解析计算门限""" return chi2.ppf(1 - pfa, df) / 2 # 除以2是因为统计量里有个N/2因子 thr_analytic = threshold_from_chi2(pfa=0.01, df=1) print(f"解析门限: {thr_analytic:.4f}") print(f"仿真门限: {thr:.4f}")

两个门限应该接近但不完全相同，因为卡方近似在大样本下才准确。N=64时已经比较接近了。如果差距大，说明样本量不够，或者模型不满足正则条件。

我一般会先用解析门限跑一遍，再用蒙特卡洛验证。如果两者差距在 10% 以内，就直接用解析门限，省掉大量仿真时间。差距大就老老实实跑仿真，同时检查模型假设是否成立。

另一个实用技巧是给统计量做归一化。不同样本数 $N$ 下统计量的分布不同，门限没法通用。把统计量除以 $N$，得到单位样本的统计量，门限就可以跨 $N$ 复用了。代价是低 $N$ 时归一化统计量的方差更大，检测性能会略降。

def normalized_glrt(x, s): """归一化 GLRT 统计量，便于跨样本数复用门限""" N = len(x) return glrt_statistic(x, s) / N

这个习惯是我踩过坑之后养成的。早期做变长帧检测时，每换一个帧长就要重新仿真门限，后来统一归一化，门限一次定好，所有帧长通用，省了大量时间。希望帮到你。

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