企业官网建设流程全解析

蚀变数据筛选这件事，说穿了就是在一堆红点里判断哪个点值得背起包上山。刚入行那会儿我以为提取算法才是难点，后来发现提取只是把门打开，真正的功夫全花在“删”上——一张工作区的羟基异常图摊开，几千个斑块从比例尺上看几乎铺满，谁都知道那不可能是矿，但要从里面挑出二三十个靠谱的，靠的不是软件默认参数。这篇东西写给正在做岩矿遥感、手里已经有一叠异常图却不知道怎么收敛的人，也写给刚接手蚀变数据处理、还在纠结阈值取几倍标准差的同行。我不会讲太多教科书，主要讲我实际项目里反复用的那套筛选逻辑、参数怎么定、以及在哪些地方容易翻车，蚀变数据筛选的完整链路会从头到尾过一遍，能直接拿去改参数复现。

1. 蚀变数据筛选到底在解决什么问题

1.1 从光谱到异常图：一条每一步都在积累误差的链路

先要把这条链讲明白，不然后面的筛选全是凭感觉。原始影像进来，第一步辐射定标把 DN 转成辐亮度，第二步大气校正转成地表反射率，第三步才是蚀变信息提取，常见的是波段比值、主成分分析、光谱角匹配、匹配滤波这几类。每一步都在往结果里注入误差，大气校正里残留的水汽和气溶胶会在 2.2 微米附近制造假吸收，正射校正做不干净会让异常斑块和地质界线错开两三百米，多景影像拼接处的辐射差异会直接在异常图上画出一条笔直的边。我见过最典型的一个案例，工作区横跨两景不同季节的影像，拼接缝在蚀变异常图上呈现为一条南北向的高值带，长度将近二十公里，看着比任何真实异常都“漂亮”，最后用直方图匹配把两景拉到同一辐射基准才消掉。

所以筛选不是独立的最后一步，它其实是对整条链路的一次反向质检。你在筛异常的同时，实际上也在检查辐射校正做没做到位、掩膜做没做干净、阈值定得合不合理。理解这一点很关键，因为它决定了你遇到问题时的排查方向——很多“异常筛选结果不对”，根子根本不在筛选环节，而在三个月前那次大气校正。

1.2 假异常的四大来源，决定了筛选的整体思路

为什么不能直接信提取结果？因为异常像元里混进了大量与矿化毫无关系的东西，我习惯分成四类。

第一类是光谱相似地物。铁氧化物异常最爱被红层砂岩、砖瓦、裸土、煤系地层干扰，这些地物的三价铁吸收特征和褐铁矿几乎重合；羟基异常则常被盐碱地、石膏、明矾石化风化壳、甚至含高岭石的耕地土壤污染。这一类的麻烦在于它们的光谱特征是真的，不是噪声，任何靠光谱阈值的方法都删不掉，只能靠空间和地质信息去压。

第二类是环境干扰。云、云影、积雪、水体、浓密植被、地形阴影，这些东西本身没有蚀变特征，但因为整体反射率畸低或畸高，在比值运算里会被放大成假高值。云影是最阴险的一种，它在可见光波段暗得像水体，在短波红外又保留了一部分地表信号，做比值法时经常在云影边缘拉出一圈异常环。

第三类是处理引入的伪异常。除了前面说的拼接缝，还有条带噪声、DEM 高程误差导致的阴影误判、掩膜边界处的卷积效应。这类异常有个共同特征：形态极其规整，要么是笔直的线，要么是紧贴影像边缘的框，肉眼一看就不像自然产物。

第四类是阈值设定不当。同一套算法，在不同景观区用同一个阈值，背景均值高的区域会大面积出异常，背景均值低的区域一个都出不来。我早期的毛病就是全工作区一个 k 值走天下，后来发现山区和平原的直方图形态差得太远，必须分区统计。

把这四类想清楚，筛选策略基本就定了：环境干扰靠掩膜删，处理伪异常靠形态和位置删，光谱相似地物靠多方法交集压，阈值问题靠分区统计和分级解决。

1.3 筛选的三层判据体系

我把整个筛选设计成三层判据，层层收敛。

光谱层管“这个像元的光谱形态对不对”，判据是吸收深度、吸收位置、比值大小、光谱角。这一层能删掉噪声和大部分非蚀变地物，但删不掉红层。

空间层管“这个斑块长得像不像地质体”，判据是面积、形状指数、长宽比、聚集密度、边界平滑度。道路是长条形，阴影是不规则大片，矿化蚀变通常呈现沿构造线分布的串珠状或团块状。这一层对删道路、删条带特别有效。

地质层管“这个位置有没有成矿条件”，判据是地层接触带、断裂缓冲距离、岩体接触晕、已知矿化点、化探异常、物探梯度带。这一层是最终收敛的关键，也是最耗时间的，因为它需要你把遥感结果和一堆非遥感资料叠起来看。

三层下来，一个典型工作区的收敛比例大致是这样：提取得到十万量级异常像元，掩膜加阈值后剩几万，斑块过滤后剩几千，聚类成几百个异常集中区，多方法交集后剩一百多，地质约束打分后剩二三十，野外验证后真正值得布钻的可能就三五个。这个漏斗每往下走一步都要留中间成果，别一次删干净，不然出了问题没法回溯。

2. 筛选前的数据底座准备

2.1 辐射与几何处理的底线要求

筛选精度有一个天花板，天花板由预处理决定。对蚀变数据来说，我给自己定的底线是三条。

反射率必须做过大气校正，且校正后同一稳定地物在不同景之间的反射率差异控制在 5% 以内。检查方法很土但有效：挑一片大面积裸岩或干河床，在两景里各取一千个像元算均值，差异超了就得回炉。没做大气的相对归一化也能凑合用，但只能做单景内部的相对判断，跨景比较会出问题。经验上，短波红外波段的校正难度比可见光高，2.2 微米附近如果有明显的水汽残留，羟基异常的背景值会整体抬高，阈值就跟着水涨船高。

正射校正的相对误差要压到一到两个像元以内。蚀变异常最终要和地质图、断裂、化探点位叠，错位三个像元以上，缓冲距离分析就全是假的。山区地形起伏大的地方尤其要注意，用低分辨率 DEM 做正射会在山脊两侧留下系统性偏移。

同一工作区尽量用同一季节、相近太阳高度角的影像。太阳高度角差 10 度，阴影面积能差出一大截，而阴影边缘恰恰是假异常的高发区。如果实在避不开跨季节，就分幅独立统计阈值，别强行统一。

2.2 干扰地物掩膜：一张一张数着做

掩膜这一步最枯燥，但它决定了后面所有工作量的下限。我的做法是把干扰地物拆成独立图层，一张一张单独生成，最后取并集，这样出了问题能立刻定位是哪张掩膜漏了。

表格里我特意没写死具体数值，因为这些东西真的没有通用阈值。植被阈值在湿润区和干旱区能差一倍，积雪阈值在高原和平原完全不同。我的做法是每次先做一张掩膜叠加在真彩色影像上，目视检查边界贴不贴合，不贴合就调，调完再看。这个过程大概要来回三四次，别嫌烦，比后面人工删假异常省事得多。

提示：掩膜边界要做 1 到 2 个像元的膨胀处理。因为植被和水体的边缘是混合像元，直接按硬阈值切会在边界留下一圈半植被半裸土的像元，这些像元在比值运算里最容易出假高值。

2.3 阈值分割方法的选型对比

异常图到二值图这一步，方法选型直接影响后面筛什么。我把常用的几种摆在一起对比。

我个人的默认组合是分区均值加 k 倍标准差做初筛，再用直方图拐点法做交叉验证。两个方法给出的阈值如果差得不多，说明数据分布比较规矩，可以直接用；如果差得远，说明分布有明显长尾或者多峰，这时候就得回去看是不是掩膜没做干净，或者工作区里混进了两种不同背景的景观单元。

主成分特征向量法这里多说一句，因为很多人用错。判读规则不是看哪个分量方差大，而是看目标矿物吸收波段对应的载荷符号是否与其他波段相反，且绝对值足够大。举个具体例子，做羟基异常时用一种短波红外为主的四波段组合做主成分分析，如果某个分量里 2.2 微米波段的载荷为负、其余三个为正，那这个分量就代表羟基吸收的增强方向，取它的负值或者取反后作为异常分量。载荷的绝对值越大，说明该波段在这个分量里的权重越高，物理意义越纯粹。

3. 三级筛选流程的完整实现

3.1 一级筛选：统计阈值加斑块形态过滤

一级筛选解决的是“像元级”的问题，目标是快速把明显不是异常的像元删掉。这一步我一般用代码跑，因为量大，人工看不过来。

import numpy as np import rasterio from scipy import ndimage # 读入蚀变分量栅格与有效区掩膜 with rasterio.open("pc_hydroxyl.tif") as src: pc = src.read(1).astype("float32") profile = src.profile nodata = src.nodata with rasterio.open("mask_valid.tif") as msrc: mask = msrc.read(1).astype(bool) # True 表示有效像元 valid = mask & (pc != nodata) mu, sigma = pc[valid].mean(), pc[valid].std() k = 2.5 thr = mu + k * sigma print(f"mean={mu:.3f} std={sigma:.3f} thr={thr:.3f}") anom = valid & (pc > thr) # 去除小于 min_px 像元的碎斑 min_px = 50 lab, n = ndimage.label(anom, structure=np.ones((3, 3))) sizes = ndimage.sum(anom, lab, index=np.arange(1, n + 1)) keep = np.zeros(n + 1, dtype=bool) keep[1:][sizes >= min_px] = True anom_clean = keep[lab] profile.update(dtype="uint8", count=1, nodata=0, compress="lzw") with rasterio.open("anom_level1.tif", "w", **profile) as dst: dst.write(anom_clean.astype("uint8"), 1)

代码不复杂，但有三处细节值得说。ndimage.label的structure参数一定要传八邻域，默认是四邻域，会把对角相连的斑块拆成两个，导致本来够面积的斑块被误删。min_px的取值和影像分辨率强相关，30 米分辨率下 50 个像元大约是 4.5 万平方米，差不多是一个小矿化蚀变带的规模；换成 15 米分辨率，同样面积就要 200 个像元，不能照搬。第三，计算mu和sigma时一定要限定在有效区内，如果把掩膜掉的云、水体也纳入统计，标准差会被拉得很大，阈值跟着虚高，异常就全没了。

形态过滤在这一步也要做。我最常用的是形状指数，定义为4πA/P²，A 是面积，P 是周长。圆形接近 1，细长条接近 0。经验取值是把形状指数低于 0.15 的斑块删掉，这一刀能砍掉大量道路、条带噪声和拼接缝。另一个指标是长宽比，用最小外接矩形算，长宽比大于 8 的基本可以判定是线性地物。这两个指标配合用，比单纯看面积有效得多。

注意：形态过滤不要一刀切太狠。真实的蚀变带如果沿断裂发育，本身也可能呈长条状。稳妥的做法是把长条状斑块单独导出来，叠加到影像上目视确认一遍，确认是道路再删，别直接按指标批量干掉。

3.2 二级筛选：多方法交集与空间聚类

一级筛完之后，异常像元数量通常会降到原来的百分之几，但斑块数量可能还有几千个。这一级要做两件事：压掉光谱相似地物的干扰，把散点聚成有地质意义的异常区。

交集策略上，我不追求三种方法全交，因为全交的保留率往往不到 10%，容易把真异常也交掉。我的经验是三种方法取两种以上，保留率大约在 30% 到 50%，这个比例既能压掉大部分单方法特有的假异常，又不至于过度收缩。具体哪三种，得看数据源：多光谱数据常用比值法、主成分分析法、光谱角匹配三选二；高光谱数据可以加上匹配滤波和矿物丰度图。

交集之后跑空间聚类，用密度聚类比较合适，因为它不需要预先指定簇的数量，还能识别噪声点。

import numpy as np from sklearn.cluster import DBSCAN pix = 30.0 # 像元大小，单位米 eps = 4 * pix # 邻域半径，约 120 米 min_samples = 8 # 邻域内最少像元数 rr, cc = np.nonzero(anom_clean) xy = np.column_stack([cc * pix, -rr * pix]) # 转为米制，y 轴翻转 db = DBSCAN(eps=eps, min_samples=min_samples).fit(xy) labels = db.labels_ print("簇数量:", len(set(labels)) - (1 if -1 in labels else 0))

这里有两个参数要定，eps和min_samples，具体怎么算我在第 4 章细说。跑完聚类，每个簇就是一片异常集中区，接下来按簇内像元总数分级：像元数最多的前 10% 定为一级异常，10% 到 30% 定为二级，其余为三级。这个分级只是工作量排序，不代表找矿价值排序，别混淆。

聚类这一步有个容易被忽略的用处：它能暴露条带噪声。如果聚类结果里出现一长串沿某一方向均匀排列、簇内像元数都差不多的异常区，十有八九是传感器条带或者地形阴影造成的。这种规律性用单个斑块很难看出来，一聚类就一目了然。

3.3 三级筛选：地质、物探、化探的交叉约束

到这一步，剩下的异常区可能还有一百多个，光靠遥感已经压不下去了，必须引入非遥感资料。我用的是一个加权打分的思路，把每个异常区的各项指标量化，最后按总分排序。

地质约束这一块，我主要看四个东西。地层方面，中酸性侵入体与碳酸盐岩或碎屑岩的接触带、火山岩机构的环状断裂带、变质岩区的韧性剪切带，这几类是最有利的。构造方面，统计异常区到最近断裂的距离，通常 500 米以内加分，1 公里以外不加分，同时看是否位于断裂交汇部位，交汇部位的权重更高。岩体方面，看异常区到岩体接触带的距离，以及接触带外侧的矽卡岩化晕宽度。已知矿化点方面，2 公里缓冲范围内有已知矿点的直接加分，这个信号最强。

物探和化探的用法要克制。化探异常圈的范围往往很大，用它做硬约束会把大部分遥感异常都圈进去，失去筛选意义；我一般只取化探的高值浓集中心，还是用距离缓冲的方式叠加。物探方面，航磁的梯度带、重力资料的密度界面、激电的高极化率区，这些信息的分辨率通常比遥感低一个数量级，只适合做背景参考，不适合做点对点的匹配。

打分表我按下面这个思路设计，权重可以根据工作区的矿化类型调整。

总分 100，实际跑下来通常能筛出二三十个 60 分以上的异常区。这个数量级是可以承受野外工作的。分数只是排序工具，我见过总分不高但野外验证见矿的情况，所以别把打分当成判决书，它只是帮你决定先上哪座山。

3.4 野外验证与置信度回写

野外这一环决定了整套流程能不能形成闭环。我的标准动作是：到异常区中心用便携式光谱仪实测三到五个点，记录坐标和光谱曲线；在吸收特征最强的位置取样，送薄片鉴定和 X 射线衍射；同时记录地表出露情况、蚀变类型、蚀变强度和分带特征。

回填的置信度分四级。确认：野外见到对应蚀变矿物，光谱实测吸收特征吻合。基本确认：见到蚀变但矿物组合与提取目标不完全一致。存疑：地表覆盖严重，取不到新鲜样品。否定：野外证实是岩性、风化壳或人工地物引起的假异常。

这些标签回填到异常区数据库里，下一轮做同类工作的时候可以直接作为训练样本或者参考阈值。我手上有一套跑了三年的样本库，每次新工作区开张，先拿样本库里的光谱曲线做一遍匹配，误报率能比从零开始低三成左右。这个积累过程很慢，但复利效应明显。

4. 关键参数怎么定：三个可复现的算例

4.1 阈值算例：k 到底该取多少

这是问得最多的问题，我给一个完整的推演过程。

假设工作区影像 8000 乘 8000 像元，共 6400 万个像元，做完掩膜后有效像元占 70%，也就是 4480 万个。某羟基主分量在有效区内的分布近似正态，均值 μ 接近 0（主成分变换后的分量均值天然接近零），标准差 σ 为 120（单位无所谓，取决于数据量纲）。

在不同 k 值下，单尾超出阈值的像元比例如下。

看这组数就明白为什么不能拍脑袋取 k=2。取 2 的话有上百万异常像元，后面的流程根本跑不动；取 3.5 只剩一万，除非成矿条件极好，否则很可能把真异常也切掉了。我的经验区间是2.5 到 3.0，对应的异常像元占比 0.1% 到 0.6%，这个比例和全球范围内已知蚀变带的面积占比大致吻合，算是有物理依据的经验值。

再往下推。取 k=2.5，得到 27.8 万异常像元。假设这些像元成簇分布，平均每个斑块 30 个像元，能形成约 9200 个斑块。用 50 像元的最小面积过滤，能砍掉大约八成，剩 1800 个左右。走三条方法的二取二交集，保留率按 40% 算，剩 700 多个。聚类之后形成 200 个左右的异常集中区。这时候数量级就进入人工可处理的范围了。

这套推演的价值在于：如果你跑出来的结果和这个估算差了十倍以上，说明某个环节出了问题，可能是掩膜没做、可能是阈值取错、也可能是主分量选错了。先算预期再跑数据，这个习惯能省掉很多返工。

分区统计这件事必须强调。同一个工作区里，如果一半是山区一半是盆地，背景均值可能差两三倍。这时候要按地貌单元或者按影像分幅分别算 μ 和 σ，各用各的阈值。我一般用滑动窗口法做局部统计，窗口大小取 500 乘 500 像元，窗口内有效像元不足 30% 就跳过。这样出来的阈值面是渐变的，比全域统一阈值贴合实际得多。

4.2 聚类算例：eps 与 min_samples 怎么定

eps的物理含义是“多远的两个像元算作邻居”。它的取值要跟地质体的尺度挂钩，不能只看像元大小。

假设影像分辨率 30 米，目标蚀变带的典型宽度 200 到 500 米，那么eps取 120 米（4 个像元）比较合适，能把宽度几百米的蚀变带内部连成一片，又不至于把相隔一公里的两个独立异常合并。如果研究的是斑岩型蚀变，蚀变晕范围可能上公里，eps可以放到 300 米甚至更大；如果做的是脉状矿化，蚀变带只有几十米宽，eps就要收到 60 米左右。

min_samples控制的是密度下限，防止把零星散点聚成簇。常用的经验规则是取min_samples ≈ ln(N)，N 是异常像元总数。27.8 万个像元，ln 值约 12.5，取 10 到 15 比较合理。这个规则背后的逻辑是：在随机分布假设下，eps邻域内出现ln(N)个点的概率很低，超过这个数就说明是聚集而非随机。

参数敏感度要实测。我的做法是固定一个参数扫另一个，画一张簇数量与参数的关系曲线，找曲线趋于平缓的那一段。簇数量对eps特别敏感，eps从 90 米加到 150 米，簇数量可能从 500 掉到 100，这段区间里结果极不稳定，说明参数选在了悬崖边上。稳妥的做法是选曲线平台上的一段。

提示：DBSCAN 的朴素实现内存占用是 O(N²)，几十万个点还好，上百万点就得注意了。两个办法，一是先按图幅分块跑再合并边界处的簇，二是先做一遍规则网格降采样，把数据压到十万量级再聚类。降采样会损失一点精度，但对簇的边界影响很小。

4.3 斑块度量与吸收深度的实际算法

除了面积和形状，我还会算三个指标。

吸收深度是最能反映蚀变强度的。算法是先在吸收波段两侧的肩部做线性内插，得到连续统，再算实测值与连续统的比值。

举个具体算例。某像元在中心波长 1.65 微米处反射率 0.30，在 2.395 微米处反射率 0.26，目标吸收在 2.205 微米处实测反射率 0.21。先内插连续统：

Rc = 0.30 + (0.26 − 0.30) × (2.205 − 1.65) ÷ (2.395 − 1.65) = 0.30 − 0.04 × 0.555 ÷ 0.745 = 0.30 − 0.0298 = 0.2702

然后 D = 1 − Rb / Rc = 1 − 0.21 / 0.2702 = 0.223。

吸收深度 0.223 属于中等偏强的羟基吸收水平，在绢云母和明矾石发育的区域常见。这个值本身没有绝对意义，因为它受颗粒度、湿度、混合像元比例影响很大，只有同一景影像内部的相对比较才有参考价值。我一般把全区吸收深度做百分位归一化，取 80 百分位作为强度分级的参考线。

异常聚集度，定义为聚类簇内异常像元数与簇外接圆面积的比值。这个指标能区分“致密蚀变中心”和“稀疏晕”，前者聚集度高，找矿意义更大。

面积密度，统计 1 公里乘 1 公里网格内的异常像元占比。这个指标适合做靶区圈定，把密度等值线图叠在地质图上，高值区直接圈成找矿靶区。

5. 实操踩坑与排查手册

5.1 高频问题速查表

下面这些是我这些年实际遇到并记录下来的，按出现频率排序。

5.2 几条不写在文档里的经验

第一条，先做减法再做加法。新手容易犯的错是拼命调算法想把异常提得更全，结果越调越乱。正确顺序是先想办法把不该出现的都删掉，让保留下来的每一个异常都能解释得通，然后再回头考虑是不是漏了。宁可漏掉几个真异常，也不要让假异常淹没真异常，因为筛选成本主要在人工判读，假异常多一个就多一份工作量。

第二条，中间成果全部留档。每一级筛选后的栅格、每一版的阈值参数、每一次的野外记录，都按日期和版本号存好。我吃过这个亏，半年前的一个工作区被要求补充报告，结果发现当时的交集阈值没记录，只能从头推一遍。现在我的目录结构是固定的：原始数据、预处理中间件、提取结果、筛选各级、打分表、野外验证，层层分目录，文件名带日期前缀。

第三条，别信软件默认参数。无论哪个光谱匹配工具，默认阈值都是针对通用场景设的，放到你的工作区大概率不合适。我的习惯是拿到默认结果后先看分布，再自己定阈值，从来不用软件给的默认二值图。

第四条，用已知区标定。每个工作区如果能找到一片已经勘探验证过的区域，哪怕很小，都拿它来标定整个流程。具体做法是把已知区的一部分作为训练，剩下的作为检验，看你的阈值能不能把已知矿化区圈出来、圈出来的面积占全区多少。这一步花不了两天，但能让后面所有参数的确定变得有据可依。

注意：工作区里如果有大面积的耕地、果园、人工林，短波红外波段的信号会被耕作层和植被残体强烈影响。这类区域建议直接整片掩掉，别试图区分。我试过用多种指数想保留耕地里的真实异常，效果很差，投入产出完全不划算。

第五条，异常的分级要服务于工作节奏。一级异常先跑，跑完回来更新置信度，再决定要不要跑二级三级。不要一次性把几百个异常全部验证，那既不现实也没必要。野外验证的排序逻辑是：先看一级异常里打分最高的，再看交通条件，最后看季节窗口。冬天雪盖和夏天植被都会影响验证效果，干旱区反而是夏季植被稀少的季节更好。

5.3 结果交付前必须过的五道检查

我自己定了五条硬性检查，任何一条不过就不能出图。

第一道，异常像元占比检查。全区异常像元占比应在 0.05% 到 1% 之间。低于 0.05% 说明阈值过严，可能漏矿；高于 1% 说明掩膜或阈值有问题，结果没法用。

第二道，形态合理性检查。把所有异常斑块用最小外接矩形框出来，快速浏览一遍，看有没有明显笔直的、规整矩形的、紧贴影像边缘的。这三种形态基本可以判定是处理伪异常。

第三道，与已知信息的一致性检查。叠加已知矿化点、化探异常、断裂构造，看异常的空间分布是否落在合理位置。如果一个已知矿点周围三公里内一个异常都没有，要么是阈值问题，要么是提取方法不适合该区的蚀变类型，必须回头检查。

第四道，野外抽检覆盖率检查。一级异常至少要抽检 30%，覆盖不同的地貌单元和异常强度等级。抽检结果里如果否定率超过 50%，说明前面的筛选逻辑有系统性偏差，要整条链重新审。

第五道，可复现性检查。把参数、代码、中间数据交给另一个人，让他独立跑一遍，看能不能复现出基本一致的结果。这一条最难，但价值最大，因为它是检验你有没有把偶然当必然的唯一办法。我见过太多结果看着漂亮、换个人重跑就出不来的情况，那种结果不能写进报告。

蚀变数据筛选这活，做久了会发现它本质上是个不断收敛的过程，从一片红点收敛到几个坐标，每一步都在用不同的信息维度做减法。文献里讲得最多的是提取算法，实际项目里最费时间的永远是筛选和验证。我个人最深的体会是，一套稳定的筛选流程，比一个花哨的提取算法值钱得多——算法可以换，流程一旦跑顺了，换个工作区、换种矿化类型，改改参数就能接着用。