企业官网建设流程全解析

简介：高光谱图像异常检测中广泛采用的RX算法，如今有了可直接运行的MATLAB实现。这份压缩包面向遥感、环境监测、资源勘探等领域的开发者与研究人员，解决高光谱数据中快速识别异常像素的需求。包内共3个文件，均为.m脚本，其中主程序负责核心检测流程，另有脚本分别完成三维高光谱立方体与二维矩阵间的格式转换，便于统计计算与结果重构。整个压缩包仅1KB，轻量易用。已有276人浏览学习，适合希望掌握RX算法原理并快速上手实验的读者。通过阅读和运行这些脚本，可以理解数据预处理、统计量计算、阈值判定及异常标记的完整链路，还能为结合PCA、SVD等改进方法提供基础代码参考，对教学演示与二次开发都很友好。

1. 高光谱 RX 异常检测为何值得再读一遍

遥感行业里，高光谱 Transformer 已经铺天盖地，但很多真实项目在第一次筛异常时，仍然选择 1987 年提出的 RX 检测器。hyperRX.zip里的三个 MATLAB 文件（hyperGRX.m、hyperConvert2D.m、hyperConvert3D.m）组成了一个极简流水线：先把 M×N×B 的立方体压成二维矩阵，再用均值向量和协方差矩阵计算每个像素的马氏距离。模型参数只有 B 维均值向量和 B×B 协方差矩阵，与图像面积无关，这让它在环境评估、灾害监测乃至工业异常检测算法中仍然被当作不可绕过的基线。适合不熟悉 RX 细节的工程师快速上手，也适合研究者拿来做对照实验。只要正确理解“异常分数不是概率”这一点，它的应用场景立刻能扩到无监督缺陷检测和数据质量筛查。

2. hyperConvert2D 与 hyperConvert3D：高光谱立方体与二维统计矩阵之间的形态转换

2.1 RX 为什么必须先解决数据布局问题

高光谱数据通常以 M×N×B 的三维数组进入算法，M、N 是空间行列数，B 是波段数。统计模型以“像素样本”为基本单位：一个像素就是一串长度为 B 的光谱向量。RX 需要计算所有像素光谱的平均值和协方差，这对三维数组来说并不直观。直接从三维数组循环做双重求和，既容易在空间维和波段维之间来回搞错索引，也没有办法直接调用 MATLAB 的 mean、cov 这类优化过的矩阵函数。

常见做法是先通过 reshape 将空间维展平，把高光谱立方体变成二维矩阵，这就是 hyperConvert2D.m 的职责。变换只改变数据的组织方式，不改变光谱向量内部的数值顺序，所以完全不损失信息。检测结束后，又要把同一份二维矩阵还原成三维图时，再由 hyperConvert3D.m 做逆变换。可以把这个过程理解成流水线上的一进一出，中间夹住的才是真正的统计计算。

如果跳过这一步，直接手写三层循环去估计协方差，多数情况下会得到一个正确的均值和一个错误的协方差，因为不同波段的空间位置错位会让协方差矩阵变得毫无物理意义。因此在打开 hyperGRX.m 之前，先确认图像数据已经被转换成二维矩阵，是排错的第一步。

2.2 两段 reshape 代码还原转换过程

hyperConvert2D.m 最常见的等价实现如下：

function [X, rows, cols, bands] = hyperConvert2D(img) % 输入：M x N x B 的高光谱立方体 % 输出：像素数 M*N 行、波段数 B 列的二维矩阵 [rows, cols, bands] = size(img); X = reshape(img, rows * cols, bands); end

对应的 hyperConvert3D.m 写成：

function img = hyperConvert3D(X, rows, cols, bands) % 输入：pixels x bands 的二维矩阵和原始空间尺寸 % 输出：还原后的 M x N x B 高光谱立方体 img = reshape(X, rows, cols, bands); end

这里的核心是 MATLAB 的 reshape 按列优先方式填充。比如原始矩阵的大小是 [rows, cols, bands]，第一波段的一部分会被优先排进 X 的前几列。若想让波段维在前，原始代码也可能写成X = permute(reshape(img, [], bands), [3, 1, 2])之后再做转置，结果得到 bands×(M*N) 的布局。

两种布局在统计上等价，但后续代码方向完全不同。建议拿到 hyperGRX.m 后用size(X)实测一下，看清楚它期望的是像素在后还是波段在后。否则计算出来的协方差矩阵可能变成了像素空间的相关结构，检测结果看似有图，实际已经失真。

2.3 布局方向、内存开销与数据类型

下表列出三种典型尺寸下的内存与协方差量级：

可以看到，协方差矩阵只取决于波段数，这就是 RX 方法能被用到较大场景的原因。内存压力基本来自读入的原始高光谱立方体本身。如果内存紧，可以在读入后立即执行img = single(img)，内存占用会降一半。单精度仅带来约 1e-6 的浮点误差，对均值和协方差的估计影响几乎可以忽略。

还有一条容易被忽略的经验：uint16原始图像不要直接参与协方差计算。uint16 的类型会让部分 MATLAB 函数结果变成双精度，但数值范围可能达到 65535，导致协方差矩阵的元素非常大。更合理的是先把数据转成 double 并缩放到 [0,1]。转换矩阵时，同一批像素在空间上的排列顺序应当保持一致，不能一会按行光栅扫描，一会按列光栅扫描，否则会有空间排列误差。

验证两个转换函数是否版本匹配，可以跑下面这个最小回归测试：

[M, N, B] = deal(4, 5, 6); img = randn(M, N, B); X = hyperConvert2D(img); imgBack = hyperConvert3D(X, M, N, B); assert(isequal(img, imgBack));

这段代码虽然只有四行，却能绕开绝大多数布局错误。当后面把 hyperGRX 的结果 reshape 回 M×N 时，如果这次断言不过，那就没有继续讨论检测结果的意义。

3. hyperGRX.m 的核心算法：均值、协方差与马氏距离

3.1 全局 RX 检测器的统计本质

hyperGRX.m 中执行的全局 RX 公式为：

δ(x) = (x − μ)^T C^{-1} (x − μ)

这里 x 是一个 B 维光谱向量，μ 是整幅图像所有像素光谱的均值向量，C 是背景协方差矩阵。把 C^{-1} 视为马氏距离的度量矩阵，RX 实际上是在寻找那些偏离背景分布中心最远的像素。偏离方向不受重视：当背景方差大时，该维度的贡献被压低；当背景方差小时，该维度的贡献被抬高。这是它比纯欧氏距离更适合高光谱异常检测的原因。

从几何上看，δ(x) 相当于先把样本 x 映射到白化空间，再计算平方模长。如果背景各波段方差差异超过两个数量级，不做白化的话，检测分数会被强方差波段垄断，弱信号波段彻底失去作用。hyperGRX.m 里的协方差逆矩阵就是承担这个白化动作的核心参数。

需要单独强调的是，RX 输出的马氏距离不是概率，也不是 0/1 标签。如果注释里写了“大于阈值就判异常”，那阈值大概率不在模型内部，而在调用脚本中。理解这一点可以避免把检测分数直接当成置信度来投入业务系统。

这里还有个自相矛盾：RX 要求背景是正常的，可一开始并不知道什么是正常。当异常像素数量较多且光谱偏移大时，μ 和 C 会被异常样本拉偏，部分弱异常反而被压低。如果检测图上只有少数几个强异常被标出，却没有中间过渡像元，这通常是掩膜效应。常见做法是先跑一遍 RX，把第一次分数高于某个粗阈值的像素剔除掉，再用剩下的像素重新估计 μ 和 C，最后跑第二遍检测。这种做法常被称为 two-pass RX，实现成本很低，却能明显改善弱异常面积较大的场景。

3.2 一个可替换的 MATLAB 实现

下面的函数是 hyperGRX.m 的逻辑等价实现：

function [scoreMap, scoreVec] = hyperGRX(img) % 输入: M x N x B 高光谱立方体 % 输出: scoreMap 为原始尺寸的检测图，scoreVec 为展平的分数 [M, N, B] = size(img); X = hyperConvert2D(img); % M*N 行, B 列 mu = mean(X, 1); % 全图均值光谱 C = cov(X); % B x B 协方差矩阵 % 添加小正则项避免奇异 reg = 1e-6 * trace(C) / B; C = C + reg * eye(B); Xc = X - mu; invC = inv(C); d2 = sum((Xc * invC) .* Xc, 2); % 逐像素马氏距离 scoreMap = reshape(d2, M, N); scoreVec = d2; end

代码中的三个关键点：mu是 1×B 的行向量，mean(X,1)的方向必须与 X 的布局一致；C = cov(X)要求 X 的每一行是样本，列是特征，这和 hyperConvert2D 输出的 pixels×bands 布局配套；inv(C)要求矩阵满秩。

如果 hyperConvert2D 返回的是 bands×pixels，那么mean(X,2)是均值，cov(X')才能得到 B×B 协方差矩阵。常见的报错是矩阵维度不匹配，根源就是布局不统一。

reg正则项用trace(C)/B归一化，再乘以 1e-6，这样可以保证在背景波段方差差异很大时，弱方差波段不会被强方差波段的扰动淹没。部分工程版本直接用C = C + 0.001*eye(B)，在反射率数据上也成立，但波长范围跨度大的传感器最好按相对值加。如果数据本身已经做了方差归一化，正则项甚至可以设成 1e-10，只解决数值舍入问题。

3.3 阈值怎么取：统计阈值、百分位阈值与自由度

卡方阈值看起来优雅但容易误用。马氏距离在样本服从多元正态分布时服从卡方分布，自由度为 B。真实高光谱反射率往往带有噪声和异常边缘，不服从理想正态，因此卡方阈值只能作为一个初始位置。工业无损检测中更多是按“虚警率不能超过缺陷比例的某个倍数”来反查阈值。比如在生产线上要求一百万像素中虚警不超过 100 个，就取百分位prctile(score, 99.99)，再人工扫描确认是否为真实缺陷。

百分位阈值的好处是不受协方差奇异的影响，因为直接对分数排序取点。卡方阈值适合数据质量高的反射率影像，两者可以交叉验证。如果同一批数据上两者差异超过一个数量级，说明背景统计假设已经崩了，优先检查数据中是否包含大量未过滤的暗行、坏像元或饱和像元。

4. 完整实验：从模拟高光谱数据到 RX 异常检测图

4.1 构造带注入异常的高光谱测试立方体

为了验证 hyperGRX.m 是否有效，建议先在一幅完全可控制的模拟高光谱数据上运行。使用以下的 MATLAB 脚本：

rng(42); rows = 80; cols = 60; bands = 30; bgMean = 0.3 + 0.01 * (1:bands); img = zeros(rows, cols, bands); for b = 1:bands img(:, :, b) = bgMean(b) + 0.04 * randn(rows, cols); end truth = false(rows, cols); for k = 1:20 r = randi([5, rows-4]); c = randi([5, cols-4]); truth(r, c) = true; img(r, c, :) = 1.2 + 0.1 * randn(1, bands); end

这里地面背景使用了分段常数均值：靠后波段的均值略高，各波段独立噪声。异常像素则在整条光谱上被拉高到 1.2 附近，和背景均值相差几十倍标准差，RX 应当能把它们全部找出来。随机位置用 randi 限制在图像内边界处，是为了避免边界附近出现半个窗口统计导致评估失真，虽然这条规则对全局 RX 影响不大，但后续换成局部 RX 时仍然适用。

4.2 运行 hyperGRX 并输出检测结果

把函数放在 MATLAB 路径下后，执行下面的脚本：

scoreMap = hyperGRX(img); [rows, cols, bands] = size(img); th = chi2inv(0.999, bands); anomalyMask = scoreMap > th; figure; subplot(1,2,1); imagesc(scoreMap); colorbar; title('RX detection score'); subplot(1,2,2); imagesc(anomalyMask); colorbar; title('Thresholded anomaly mask'); labels = double(truth(:)); scores = scoreMap(:); [~, ~, ~, AUC] = perfcurve(labels, scores, 1); fprintf('AUC = %.4f\n', AUC);

perfcurve是 MATLAB Statistics and Machine Learning Toolbox 的函数，返回的多个参数中第四个才是 AUC。第一个和第二个分别是假正率、真正率，第三个是阈值向量，第四个是曲线下面积。真实业务中如果缺少这个工具箱，也可以用 for 循环遍历 1000 个阈值，统计真阳性率与虚警率到二维数组，再求梯形积分。

正常运行后，AUC 应当非常接近 1，scoreMap 上异常位置会呈现明显亮点。若 AUC 低于 0.98，优先检查函数边界布局是否一致、协方差正则项是否过大。正则项过大会把所有分数压缩到一个小范围，使异常与背景的对比变弱；通常表现为整幅 scoreMap 偏白。

如果 AUC 始终是 1.0000，需要把异常注入强度调弱，比如把异常光谱改成1.2 + 0.5*randn(1, bands)，或把平均偏移从 1.2 降到 0.8。太强的异常只检验得出“明显异常”这种简单情况，检验不出算法边界。

4.3 高光谱如何转反射率：DN 值对 RX 检测效果的影响

前面模拟数据使用的量纲范围是 0~1，真实高光谱数据则多以 DN 或辐射亮度存储。直接把 DN 丢进 hyperGRX 的风险在于，不同波段的传感器增益可能相差数十倍。增益较大的波段会主导协方差矩阵，导致马氏距离主要衡量该波段的噪声水平，而真正的异常光谱特征反而被挤到次要位置。

高光谱如何转反射率没有一个通用的公式，常见流程是：

1. 原始 DN 值通过辐射定标参数转换为辐亮度；
2. 去除大气影响，把辐亮度转换为地表反射率；
3. 对反射率做归一化或标准化，确保各波段处在近似量纲下。

如果暂时没有大气参数，也可以在 MATLAB 里做最小一致性处理：

img = double(rawImg) / 65535; % 假设原始为 16bit img = img / max(img(:)); % 全局归一化

这段代码不是严格反射率，它能解决数量级差异问题，却无法消除气溶胶和散射造成的波段间误差。这里仍建议查看影像头文件中的反射率刻度因子，如果没有就用归一化作为退路。在比较不同时段、不同太阳高度角的图像时，只使用全局归一化很可能把光照方向的一部分误判为异常像素，所以环境监测任务里最好还是做准反射率。

4.4 异常检测中的常见失败模式

运行高光谱立方体时，第一次失败的常见原因有三个。

第一，图像中含有大量暗行、坏像元和未校准像元。它们的光谱向量是 0 或饱和值，会在协方差矩阵中形成一个强轴，把正常目标推开。常见做法是先根据每个像元的辐射值做一次 mask，把 0 和饱和值排除在统计估计之外，然后再对剩余像元求均值和协方差。

第二，X 的行/列排列顺序不一致。尤其在把多个分段拼接到同一矩阵时，忘了按原始空间光栅顺序排列，scoreMap reshape 后会出现空间错位，明明异常点出现在真实位置附近但总是偏几行。

第三，混淆了批量统计和全图统计。全图 RX 估计出的协方差本来应该占满 B×B=200×200 的矩阵；若误把二维矩阵的转置当作输入，cov 操作返回的可能是 MN×MN 的矩阵，内存直接爆炸。遇到内存不足时，先检查cov(X)的矩阵尺寸是不是 B×B。若是，内存不足的根源只可能是原始像素矩阵太大，这时把大图切成若干 patch，每个 patch 单独检测，再合并分数图即可。

5. 进阶：局部窗口、低秩背景与更可靠的验证技巧

5.1 局部 RX 代替全局 RX

全局 RX 在背景光照不均匀时不稳定。比如森林、城区、水体混合在一起，整幅图的均值和协方差被多种地物拉扯，异常分数会呈现大范围的纹理响应。更好的做法是给每个像素套一个局部窗口，只用窗口内的像元估计 μ 和 C，再计算中心像素的马氏距离。局部 RX 虽好，但窗口太大会丢失局部性，窗口太小协方差矩阵又会严重奇异，常见经验是把窗口设为 11×11 或 25×25，具体取决于地物尺寸。

简单实现时不需要直接对每个像素重新计算整窗口的协方差，可以先取整个场景中很少一部分有代表性的像元估计全局协方差，再用局部均值替换全局均值。虽然从统计上讲不够严格，但许多工程场景下它比全局 RX 明显更稳。若追求完整滑动窗口，可以用im2col把邻域像元取成列，再逐列计算马氏距离；这种方式内存开销大，适合小图算法验证。

5.2 SVD 降维后运行 RX，解决协方差奇异性

如果波段数 B=200，可用像素却只有几万个，协方差矩阵不一定奇异，但背景强相关会让弱信号被数值噪声掩盖。常见做法是先对背景做 SVD/PCA 降维，再在降维后的特征空间执行 RX。

X = hyperConvert2D(img); mu = mean(X, 1); [~, ~, V] = svd(X - mu, 'econ'); k = 20; % 保留主成分数 Xp = (X - mu) * V(:, 1:k); % 用 Xp 估计 mu_k 和 Ck，再按 RX 公式计算马氏距离

降维后的 B 从 200 变成 20，可以让协方差矩阵更稳定。需要注意的是，RX 的灵敏度会因低秩投影剔除掉一部分非常微弱的异常，所以异常信号如果只存在于很小波段子集而不体现在主成分上，降维可能把它的信号推入残差。验收到这一步时，必须保留一小块不参与降维的原始维度异常注入样本，用来对比降维前后的召回率。

5.3 一种不会被虚警骗过去的验证技巧

最后提供一个实战验证技巧：不要逐像素对齐去计算正确率。真实异常边缘通常没有精确标签，逐像素比对会把偏离一个像素的检测视为错误，导致召回率看起来奇低。做法是先在真值掩膜上做一次 3×3 膨胀，再统计每个异常簇的响应。检测分数图经阈值分割后，凡是落入膨胀区域的连通域就算命中；命中数量除以真实异常簇数量才是召回率。每一条命中物再按 100×100 的局部窗口提取均值，评价信号强度。这个方法比逐像素指标更贴近工业异常检测算法在实际产线上的验收逻辑。具体做法是：把真值掩膜以 3×3 做一次膨胀，命中膨胀区域即算检出；在检出的异常簇里数簇数，而不是逐像素打分。

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