手里攒着一批多波段 TIFF想把它们合成一张综合评价图最卡人的地方往往不是代码而是权重——可见光波段给 0.3 还是 0.25植被指数和建筑指数谁更重要凭经验拍脑袋定完评审一问依据是什么就哑了。熵权法处理 TIFF 图像这条路我用了几年核心价值就在于把权重从主观判断变成由数据自身离散程度推出来的客观量某个波段在所有像元上差异越大它携带的区分信息越多权重就越高。这篇内容讲的是一套完整可复现的流程——从多波段 TIFF 的读取清洗到熵权法的逐像元改造再到把权重图写回带坐标的 GeoTIFF中间那些让结果悄悄跑偏的坑我也会一个个拆开。适合手上有栅格数据、做过基础 numpy 操作、想把叠图分析做得更讲道理的人看也适合刚接触多指标评价、还在纠结权重怎么定的朋友。1. 多波段 TIFF 的赋权困境与熵权法的切入点1.1 叠图分析里最难交代的那一步多波段 TIFF 最常见的用法是做加权叠置把若干指标波段各自归一化乘上权重相加成一张综合得分栅格。这套流程本身不难rasterio 读进来、numpy 算一算、再写出去几十行就能跑通。真正让人反复返工的是权重表。我见过太多项目的权重是专家打分 微调凑结果先按经验给一组数跑出来发现某个区域得分偏低不合理回头把对应波段权重调高一点再跑一遍直到图看起来顺眼。这种做法当然能出图但它有个致命问题——结论是被结果反向塑造的你无法证明这张图反映的是数据本身的结构还是你心里的预期。更麻烦的是波段数量一多主观定权就彻底失控。五六个波段还能凭直觉排个序十几个波段比如多光谱加几个纹理特征、几个地形因子时人脑根本没法同时权衡这么多维度的重要性最后往往退化成平均分配或者挑几个重要的给大权重剩下的象征性给一点。平均分配听着公平实际上等于宣布所有波段信息量相同这几乎不可能成立——一个几乎全图恒定的波段和一个剧烈起伏的波段凭什么同等对待。熵权法切入的正是这个位置。它不问你这个指标重不重要只问这个指标在样本间的分布有多不均匀。分布越集中说明它提供不了区分度权重自然压低分布越分散说明它能把样本拉开权重自然抬高。整个推导过程是纯数学的任何人拿同样的数据都能算出同样的权重可复现、可审计。这一点在需要写报告、需要向别人解释为什么是这个权重的场合价值远大于精度上的那点差异。1.2 熵权法在多波段影像上的独特适配性熵权法原本是评价领域的通用方法用在影像上之所以顺手是因为两者结构天然同构。标准熵权法要求一个 n×m 的决策矩阵n 是样本数、m 是指标数一张 m 个波段的 TIFF本质上就是一个像元数 × 波段数的巨大矩阵。每个像元是一个待评价对象每个波段是一个评价指标这套映射几乎不需要额外解释。由此衍生出两种用法差别很大选错了结果会完全不同。用法样本粒度权重形态适用场景计算量全局熵权全图所有有效像元一组标量权重m 个整幅影像统一评价、区域横向对比中等可抽样估计局部熵权每个滑动窗口内的像元一张权重栅格m 层空间异质性强、需要自适应高随窗口增大线性上升全局熵权的逻辑是整幅图共用一套价值观好处是不同区域、不同时相的结果可以直接比较坏处是它假设指标的重要性在空间上不变。局部熵权则让每个窗口自己算一套权重能捕捉城市区域建筑指数更重要、山区植被指数更重要这类空间分异代价是计算量暴涨而且窗口内样本少时权重会剧烈抖动。我的经验是做区域间横向对比、做时间序列分析用全局熵权把权重固定死做单幅图的内部结构刻画、做异常区识别用局部熵权。两者甚至可以叠起来用——先用全局权重出一版基准图再用局部权重看哪些地方和基准偏离最大偏离本身就是一个很有价值的信息层。1.3 这个方法的能力边界要先说清楚熵权法不是万能的用之前得认清它几个固有短板否则容易过度解读。第一它只反映数据内部的离散程度不反映指标的实际意义。假设某个波段里混入了传感器噪声噪声让数值剧烈波动熵权法反而会给它很高的权重——因为离散度大。这是最典型的翻车场景所以前置的噪声抑制和异常值处理绝对不能省。第二它对指标间的相关性完全不敏感。三个高度相关的波段会被各给一份权重等于这个维度被重复计权三次。如果波段间相关系数普遍在 0.9 以上要么先做主成分降维要么手动合并不能直接扔进去。第三权重是相对的、依赖样本集的。换一批像元、换一个研究区权重就会变。所以权重表必须和数据集绑定说明不能当成放之四海皆准的常数。把这三点想在前面后面很多争论根本不会发生。2. 拆开熵权法的计算链条2.1 指标方向统一先解决越大越好还是越小越好熵权法第一步是归一化而归一化的公式取决于指标方向。这一步搞错后面全废——我在第 5 节会专门讲这个坑有多隐蔽。多波段影像里三类方向的波段都可能出现正向指标越大越好植被指数、水体指数、亮度、纹理能量。归一化用(x - min) / (max - min)。负向指标越小越好地表温度做热岛评价时、坡度做建设适宜性时、建筑密度做生态评价时。归一化用(max - x) / (max - min)。区间型指标落在某个区间最好比如适宜坡度 0~8 度太陡不行、太平也可能排水不畅。这类要用分段公式先设定理想区间[a, b]落在区间内给满分超出后按距离衰减。实际操作里我习惯在脚本开头维护一张波段-方向映射表用字典写死避免中途忘掉。这张表同时充当文档别人接手一眼能看懂每个波段的处理逻辑。# 波段索引 - (方向, 参数) BAND_RULES { 0: (positive, None), # NDVI 类 1: (positive, None), # 纹理能量 2: (negative, None), # 地表温度 3: (interval, (0.0, 8.0)), # 坡度理想区间 0~8 度 }归一化之后还有个细节会出现 0。熵权法下一步要取对数ln(0)直接炸。标准做法是整体平移一个极小量业界常用x 0.999 * x 0.001或者x x 1e-6。前者把数据压缩到 [0.001, 1] 区间平移量固定不同批次之间可比性更好后者更轻但极小值的影响在样本量大时会被放大。我倾向用前者理由是它的扰动是有界的、可预期的。2.2 比重矩阵、信息熵、差异系数三步推导归一化完成后进入核心三步。假设有 n 个有效像元、m 个波段矩阵记作Xn×m元素x_ij。第一步算比重。对每一列每个波段单独归一化到占比p_ij x_ij / Σ_i x_ij这一步的含义是在这个波段的总量里第 i 个像元贡献了多少份额。所有像元的份额加起来等于 1。第二步算信息熵。经典香农熵公式e_j -k * Σ_i (p_ij * ln(p_ij))其中 k 1 / ln(n)系数k的作用是把熵值标准化到 [0, 1] 区间这样不同样本量下的熵可以横向比较。如果所有像元在该波段上取值完全相同p_ij全为1/n代入后e_j 1熵取最大值说明这个波段毫无区分能力。第三步算差异系数和权重。d_j 1 - e_j w_j d_j / Σ_j d_jd_j是差异系数也叫信息效用值熵越小它越大。最后按列求和归一化得到一组和为 1 的权重。整个链条只有加、乘、对数三种运算没有任何迭代和优化所以结果完全确定——这也是它比某些机器学习赋权方法更适合做报告的原因你能把公式贴在附录里任何人都能复算。2.3 手算一遍三个像元两个波段公式看着抽象用最小规模过一遍就清楚了。设两个波段、三个像元已完成归一化且已加 0.001 平移像元1: [0.801, 0.101] 像元2: [0.501, 0.501] 像元3: [0.101, 0.901]波段 1 的列和为 1.403比重为[0.571, 0.357, 0.072]。代入熵公式ln(3) 1.0986k 0.9102Σ p*ln(p) 0.571*(-0.560) 0.357*(-1.030) 0.072*(-2.632) -0.320 - 0.368 - 0.190 -0.878 e_1 -0.9102 * (-0.878) 0.799 d_1 1 - 0.799 0.201波段 2 的列和为 1.503比重为[0.067, 0.333, 0.600]对称地算出来e_2 0.799, d_2 0.201两个波段差异系数相同权重各 0.5。这个例子故意造成对称说明一件事熵权法衡量的是分布形状不是数值大小。波段 1 的绝对值普遍偏小、波段 2 普遍偏大但两者的相对分布结构一样权重就一样。很多人第一次用会误以为数值大的波段权重大这是典型误解必须先纠正。再换一组数据让波段 2 变成几乎恒定[0.501, 0.502, 0.503]。它的比重接近[0.333, 0.334, 0.333]熵逼近 1差异系数趋近 0权重会被压到几乎为零。这就是熵权法自动淘汰无用指标的直观体现。2.4 从标量权重推广到逐像元权重全局权重算完只是半程。如果要做局部熵权需要把上面的流程搬进每个滑动窗口。窗口大小k×k内样本数是k²扣除无效值后可能更少指标数还是波段数 m。每个窗口算出一组权重把它赋给窗口中心像元滑完一遍就得到 m 层权重栅格。这里有个隐含约束常被忽略窗口内的有效样本数必须显著大于波段数。熵公式里的k 1/ln(n)当 n 很小时熵的估计极不稳定。n2 时ln(2)0.693两个点就能算出熵但完全没有统计意义。我的经验阈值是n ≥ 5m即窗口内有效像元至少是波段数的五倍若波段有 6 个窗口至少 30 个有效像元对应 7×7 窗口去掉边缘和无效值后大概能剩 30 左右。窗口再小权重图会变成一片噪声看着花哨实际不能用。3. 把 TIFF 读进来并清洗干净3.1 三条读取路线怎么选Python 生态里读 TIFF 主要有三条路各有适用边界选错了后期会很别扭。库优势局限适用tifffile纯读写速度快支持 BigTIFF、多页不解析地理坐标纯图像处理、无空间参考rasterio完整地理元数据、分块读写、重投影依赖 GDAL安装稍重遥感栅格、需要坐标PIL/imageio轻量适合 8 位 RGB多波段支持弱精度受限普通图片做综合评价一般都要保留空间参考因为最终结果要叠到地图上用所以rasterio是默认选择。它的open()返回的对象带profile里面存着坐标系、仿射变换、数据类型、分块大小最后写结果时直接复用这份 profile能省掉一大堆坐标对齐的麻烦。import rasterio import numpy as np with rasterio.open(input_stack.tif) as src: profile src.profile.copy() arr src.read() # 形状 (bands, rows, cols) nodata src.nodata print(arr.shape, arr.dtype, nodata)有个细节src.read()会把整幅图读进内存。一张 10000×10000、6 波段、float32 的影像约 2.4 GB勉强能塞进内存但要是 30000×30000直接就是 21 GB机器扛不住。这种情况必须走分块。3.2 NoData 与异常值最容易被跳过的一步原始 TIFF 里几乎必然存在无效区域影像边缘的黑边、云遮挡、水体掩膜、重采样产生的空洞。它们的值可能是-9999、0、nan也可能是某个奇怪的哨兵值。这些值如果不处理后果有两个层面。表层是数值污染——一个-9999塞进归一化min直接被拉到-9999整幅图的动态范围被压缩到几乎为零所有像元挤在一起熵权法算出来的权重毫无意义。深层是统计污染——即便你用了稳健的 min/max无效像元进入比重矩阵后会把某一列的总和撑大其他像元的份额被稀释权重计算整体失真。处理流程我一般这么走按src.nodata生成第一层掩膜叠加物理合理范围掩膜比如反射率应在[-0.1, 1.2]温度应在[-50, 70]超出的一律判无效叠加np.isfinite过滤 NaN 和 Inf三个掩膜取交集得到最终的valid_mask。valid np.ones(arr.shape[1:], dtypebool) if nodata is not None: valid ~np.all(arr nodata, axis0) valid ~np.any(arr nodata, axis0) valid np.all(np.isfinite(arr), axis0) valid np.all((arr -0.1) (arr 1.2), axis0) # 按实际物理量调整注意这里用的是np.any和np.all的区别只要任何一个波段无效这个像元就整体剔除。理由很实在——综合评价要求每个像元在各维度上都有值缺一个维度就没法算加权和与其插补引入人为误差不如直接标为无效最后在成果图上留白反倒更诚实。3.3 大文件的分块策略rasterio的block_windows()会按文件内部的分块组织tiling返回窗口这个分块在写入时就是按原始布局生成的读取效率最高。对只读需求也可以自己按固定尺寸切from rasterio.windows import Window BLOCK 2048 with rasterio.open(input_stack.tif) as src: for row_off in range(0, src.height, BLOCK): for col_off in range(0, src.width, BLOCK): win Window(col_off, row_off, min(BLOCK, src.width - col_off), min(BLOCK, src.height - row_off)) block src.read(windowwin) # ... 处理 ...分块有两个坑必须处理。第一个是边界效应如果做局部熵权块边缘的窗口需要邻块数据否则窗口不完整权重会在块与块之间出现明显断层。解决办法是读取时外扩half_win像素的重叠区算完再裁掉。第二个是全局统计量依赖全局熵权需要全图的min/max和列和这些必须先在完整数据上统计出来不能边分块边累加除非你严格按可加性设计比如只累加列和但 min/max 需要单独一遍扫描。所以标准做法是两遍扫描第一遍只读元信息算统计量第二遍才正式处理。3.4 数据类型与精度取舍原始数据可能是uint16如 0-10000 的反射率放大值而归一化后必然是浮点。中间过程的 dtype 选择会影响结果精度也影响内存。我的习惯是读取后立刻转float32全流程用float32只有熵值求和这种累积量用float64累加器。理由很实际——float32只有约 7 位有效数字样本量上百万时Σ p*ln(p)的累加误差可能到 1e-5 量级虽然对权重影响不大但如果直接在整个大数组上做np.sum误差会累积得更明显。numpy 的np.sum对float32输入默认用成对求和pairwise summation精度已经不错但遇到关键结果我还是会显式.astype(np.float64).sum()。写回时反过来如果成果只是给人看的得分图float32足够如果要参与后续定量分析用float32存 GeoTIFF 并设置compressdeflate体积能压到原始的三分之一左右。4. 逐像元熵权的向量化写法4.1 双层循环为什么必须干掉最直白的实现是两层 for 循环遍历每个像元把该像元所有波段的值取出来当成一个样本算熵权——但这是错的。熵权法的样本维度是参与比较的对象集合单个像元只有一行数据n1ln(1)0公式直接退化。所以逐像元的正确理解是窗口内所有像元构成样本集窗口中心像元获得这套权重下的综合得分。按这个理解写朴素实现就是外层遍历行、内层遍历列每个位置切片取窗口算一遍完整熵权流程。一张 2000×2000 的图就是 400 万次窗口运算每次内部还有若干次 numpy 小数组操作Python 解释器开销会把它拖到几十分钟甚至几小时。这个量级必须向量化。4.2 用滑动窗口视图替代切片numpy.lib.stride_tricks.sliding_window_view是关键工具。它不复制数据只通过调整 stride 构造一个视图形状从(H, W, C)变成(H-k1, W-k1, k, k, C)。之后所有运算都是整块数组操作走的是底层 C 循环。from numpy.lib.stride_tricks import sliding_window_view def entropy_weights_local(arr, valid, win7, eps1e-3): arr: (C, H, W) float32已归一化到 [0,1] valid: (H, W) bool win: 窗口边长奇数 return:(C, H, W) 权重图 (H, W) 综合得分 C, H, W arr.shape pad win // 2 # 边缘用 0 填充后续用 valid 掩膜排除 a np.pad(arr, ((0, 0), (pad, pad), (pad, pad)), modeconstant) v np.pad(valid, ((pad, pad), (pad, pad)), modeconstant) # (H, W, k, k, C) win_a sliding_window_view(a, (win, win), axis(1, 2)) win_v sliding_window_view(v, (win, win), axis(0, 1)) win_a np.moveaxis(win_a, -1, -2) # (H, W, C, k, k) win_v win_v[..., None, :, :] # (H, W, 1, k, k) n win_v.sum(axis(-2, -1)) # (H, W, 1) n_safe np.maximum(n, 2) x np.where(win_v, win_a, 0.0) col_sum x.sum(axis(-2, -1), keepdimsTrue) # (H, W, C, 1, 1) col_sum np.where(col_sum 1e-12, 1e-12, col_sum) p np.where(win_v, x / col_sum, 0.0) plogp np.where(p 0, p * np.log(np.maximum(p, 1e-12)), 0.0) e -plogp.sum(axis(-2, -1)) / np.log(n_safe) # (H, W, C) d 1.0 - e d np.where(d 0, 0.0, d) # 数值误差兜底 d_sum d.sum(axis-1, keepdimsTrue) d_sum np.where(d_sum 1e-12, 1e-12, d_sum) w d / d_sum # (H, W, C) # 综合得分中心像元值 × 权重对波段求和 center np.moveaxis(win_a, -1, -2)[..., win // 2, win // 2] # (H, W, C) center np.where(win_v[:, :, :, win // 2, win // 2], center, 0.0) score (center * w).sum(axis-1) # (H, W) w_out np.moveaxis(w, -1, 0) # (C, H, W) return w_out.astype(np.float32), score.astype(np.float32)这段代码有几个刻意的设计值得说明。np.where(win_v, win_a, 0.0)把所有无效像元在运算前清零这样列求和只会统计窗口内的有效样本不需要在公式里额外加分支。n_safe np.maximum(n, 2)是为了防止窗口内一个有效像元都没有时除以ln(0)或ln(1)这种情况下算出来的权重其实是垃圾值最后会由valid掩膜统一剔除所以没必要在计算中做精细判断只要不产生 inf/nan 污染邻居就行。d np.where(d 0, 0.0, d)这行看着多余其实很有必要。理论上e_j ∈ [0,1]但浮点误差可能让e_j略微超过 1d变成 -1e-8之后归一化时这个微小负值会被放大成奇怪的权重分布。加个下限钳位就能彻底杜绝。center那段的取值逻辑中心像元在窗口切片里的位置是(win//2, win//2)前面moveaxis之后维度顺序是(H, W, C, k, k)所以切片是[..., pad, pad]。这里容易写错索引调试时建议先用一个 5×5 全 1 数组验证一下中心取值是否对得上。4.3 内存与速度的实测取舍上面这段在 2000×2000、6 波段、7×7 窗口下内存峰值约 3~4 GB——因为sliding_window_view虽然是视图但后面np.where、除法、log都会生成同尺寸的新数组每个中间结果都是H×W×C×k×k的float32。7×7 就是 49 倍膨胀600 万像元 × 6 波段 × 49 × 4 字节 ≈ 7 GB非常吃内存。三个优化手段按性价比排序降窗口尺寸。从 15×15 降到 7×7内存直接降 4.6 倍。5×5 往往就够用前提是波段数不超过 3~4 个且满足n ≥ 5m。分块 重叠。按 512×512 分块每块外扩pad像素算完裁剪。内存峰值只跟块大小有关与整图无关。减少中间变量。用np.einsum替代显式乘加或者把p*log(p)写成原地操作能省一两个临时数组但收益有限且代码可读性下降除非确实卡内存否则不建议。速度方面纯 numpy 版本在 8 核机器上跑 2000×2000 大约 20~40 秒可接受。如果追求更快numba的njit(parallelTrue)能再快 5~10 倍但要把循环写回去调试成本上升。我的建议是先跑通 numpy 版本验证结果正确确认逻辑无误后再考虑用 numba 加速顺序别反。4.4 写回 GeoTIFF 时容易丢的东西结果写回最省事的写法是复用输入 profileout_profile profile.copy() out_profile.update( count1, dtypefloat32, nodata-9999.0, compressdeflate, predictor3, tiledTrue, blockxsize256, blockysize256, ) score np.where(valid, score, -9999.0).astype(float32) with rasterio.open(entropy_score.tif, driverGTiff, **{k: v for k, v in out_profile.items() if k ! driver}) as dst: dst.write(score, 1)几个必须检查的点坐标不能丢。如果只复制了crs和transform却忘了height/width或者分块处理时没做坐标偏移成果叠到地图上会整体偏移。分块写回时要用dst.write(block, 1, windowwin)指定窗口不能让 rasterio 自己去猜。predictor3只对浮点有效。predictor2是给整型的混用会让压缩失效甚至数据错乱。NoData 值要一致。输出声明nodata-9999实际数组里也要真的填-9999不能只在元数据里写、数据里留着 NaN——很多 GIS 软件会因此把整幅图判为无效。权重图建议单独存。m 层权重图对分析很有价值能看出哪些区域哪个波段占主导别只存综合得分把中间过程扔了。5. 五个让结果白算的坑5.1 方向搞反结果看起来很有道理但完全颠倒这是最隐蔽的坑因为出图不会报错甚至配色出来还挺好看。我在一个热环境评价的项目里踩过地表温度波段本该是负向指标写规则表时手滑写成正向最终成果图显示中心城区环境质量最优——因为温度高的地方被当成高分。图看着平滑、对比度也不错差点就这么交出去了是同事顺口问了句市中心怎么是最优区才发现。排查链路其实很简单但要在出图前主动做打印每个波段的归一化结果直方图正向指标和负向指标的直方图应该互为镜像抽 10 个已知类别的样本点比如确定的城区点、林地五点看它们的归一化值是否与预期方向一致计算每个波段与综合得分的相关系数如果某个波段相关系数为负而它被标为正向立刻报警。第 3 条我后来做成了脚本里的断言跑完自动检查省得靠肉眼。for j in range(C): r np.corrcoef(norm[j][valid], score[valid])[0, 1] direction BAND_RULES[j][0] if direction positive and r -0.05: print(f[警告] 波段{j} 相关系数 {r:.3f}方向可能标反) if direction negative and r 0.05: print(f[警告] 波段{j} 相关系数 {r:.3f}方向可能标反)5.2 常数列与 log(0)两个都会让权重凭空消失常数列的典型来源是掩膜波段、全图同一值的分类波段、或者某个波段在裁剪后只剩下一个取值。它的熵为 1、差异系数为 0本身处理起来没问题权重为 0 是合理的但麻烦在于它会让Σ d_j变小其他波段的权重被按比例放大——如果只有一个波段有效而其他全是常数列那个波段权重直接等于 1综合得分就等于它自己看着合理其实毫无信息。log(0)的来源是归一化后出现精确的 0。虽然前面加了 0.001 平移但如果某些像元的值恰好等于列最小值x会是 0.001 而非 0通常安全。真正的风险在别处col_sum为 0。如果某个窗口内该波段全为 0比如有效像元都被掩膜掉了窗口内只剩几个有效点而它们的值恰好为 0p会变成 0/0 nannan 会扩散。防御方式是三层col_sum做下限钳位、np.log前用np.maximum(p, 1e-12)、算完统一np.nan_to_num。这三层看着冗余但真跑大图时能救命——我曾因为漏了第三层nan 从一个窗口扩散到整张权重图最后输出全黑排查了两个小时。5.3 归一化被离群值带偏熵权法对极值极其敏感因为归一化的分母是max - min。一个像元上的异常高值传感器坏点、镜面反射、火点就能把整个波段的动态范围拉大导致其余 99.9% 的像元被压缩到很窄的区间熵值偏高、权重偏低这个波段等于被废掉。处理手段有三种按侵入性递增方案做法优点风险分位数截断用 2% 和 98% 分位代替 min/max简单效果明显超出部分被钳到边界损失极端信息稳健标准化用中位数和 IQR 替代均值方差抗离群会改变分布形态熵值含义略变先滤波再算中值滤波或多时相合成去异常从源头解决引入额外处理步骤可能损失细节我一般用分位数截断2%~98% 或 1%~99% 视数据质量定。判据是截断后重新看直方图如果分布从一根尖峰 长尾变成近似对称的钟形说明异常值确实被压住了如果截断后分布明显被削平说明分位点选得太激进要往回收。有个容易忽略的细节分位数必须在全图有效像元上统计不能分块各自统计否则不同块的归一化基准不一致拼起来会有明显的块状接缝。这个现象我第一次遇到时以为是分块代码写错了检查了半天才发现是统计口径的问题。5.4 窗口尺度和波段数量不匹配局部熵权的窗口大小不是随便定的它受两个约束同时限制。下限来自样本量窗口内有效像元数要满足n ≥ 5m。6 波段的最小窗口约 7×749 个像元扣除无效后通常剩 35~45。如果强行用 3×3只有 9 个样本算出来的权重会在相邻像元之间剧烈跳变权重图看起来像随机噪声。上限来自空间分辨率窗口不能大到抹平你要观察的空间结构。做城市内部的功能区识别几十米分辨率下窗口超过 30×30约 1 公里就会把街区尺度差异全抹平权重图变成一片均匀的色块失去局部自适应的意义。中间的取舍只能靠试。我的做法是先固定波段数跑 5×5、9×9、15×15 三组看综合得分的空间变异系数标准差/均值随窗口的变化曲线曲线从陡降转为平缓的那个拐点通常就是合适尺度。这个方法不算严谨但比拍脑袋强。还有个实操细节窗口边长用奇数。偶数边长时中心像元没有唯一定义不同库的默认处理不一致容易出现半像素偏移。5.5 权重图没人看得懂呈现方式的问题算完 m 层权重图如果直接扔给别人多半得到的反馈是这花花绿绿的是什么。权重图的解读需要引导。有效做法是生成一张主导波段图对每个像元取权重最大的波段索引用分类配色显示。这样一眼能看出东部这片全是第 2 波段主导西部是第 4 波段主导直接对应到空间格局上。再配一张权重堆叠面积图随机抽 5000 个点按行堆叠展示各波段权重的整体占比和波动范围。综合得分的呈现也有讲究。原始得分范围通常在 [0.1, 0.9]直接线性拉伸会显得灰蒙蒙。我的处理是先按 2%~98% 分位做拉伸再用自然断点法Jenks分 5 级配色这样既保留了对比度分级又符合数据的实际聚集结构。分级数不要超过 7 级人眼分辨不了那么多。配色上尽量避开红绿同时高饱和度出现的方案色觉障碍人群看起来会混淆。用蓝-黄、紫-橙这类对立色对更稳妥。6. 结果可信度的交叉验证6.1 权重合理性看三个指标熵权算完先别急着出图花十分钟做三项自检。第一权重排序是否符合常识。如果一个明显充满噪声的波段权重排第一说明噪声没压住如果一个信息量很足的波段权重接近 0说明它可能被归一化压扁了。这里的常识不是主观偏好而是基于数据质量的判断——你能从直方图、异常值比例上找到依据。第二权重分布的集中度。计算权重的变异系数。如果所有波段权重都在1/m附近变异系数 0.1说明各波段离散程度接近熵权法没提供多少信息量此时用等权法反而更稳、更好解释。如果某个波段权重超过 0.5要特别警惕检查是不是其他波段的归一化出了问题。第三敏感性。把权重上下扰动 ±10%重新算综合得分看排序前 10% 的像元集合变化多大。变化很小说明结果稳健变化剧烈说明结果对权重高度敏感此时权重的小误差会放大成结论的大误差需要慎重表述。6.2 和等权法、PCA 做对照单一方法的结论总是让人不放心我习惯同时算三个版本做对照。等权法是最简单的baselinew_j 1/m。它假设所有波段同等重要虽然粗暴但胜在无偏。主成分法取第一主成分的载荷绝对值归一化作为权重它捕捉的是方差最大的方向跟熵权法思路不同但目标相近两者结果接近时可信度大幅提升。相关系数法用各波段与综合得分的相关系数归一化作为权重属于结果导向的赋权适合做验证而不适合做定权。三者的排序一致性用 Spearman 秩相关衡量。我实测的经验值是熵权与 PCA 的秩相关在 0.85 以上说明数据的主结构清晰结论可靠低于 0.6 就要回头查数据通常是波段间共线性太强或者方向标错了。from scipy.stats import spearmanr r1, _ spearmanr(score_entropy[valid], score_equal[valid]) r2, _ spearmanr(score_entropy[valid], score_pca[valid]) print(f熵权 vs 等权: {r1:.3f}, 熵权 vs PCA: {r2:.3f})注意这个比较要在大样本上做抽 10 万个有效像元就够别在全图上跑spearmanr那个函数对百万级样本很慢。6.3 用已知参照点做落地检验纯统计验证之外最有说服力的是拿真实参照点检验。做法是找一批已知属性的位置——比如几处确认的高质量绿地、几处确认的重度扰动区、几处确认的水体——看它们的综合得分是否落在预期区间。这个检验不需要很多点二三十个就足以暴露系统性问题。如果高质量区的得分中位数低于扰动区无论统计指标多漂亮都得推翻重来。我在一个项目里就是靠这个检验发现了负向指标漏标的问题统计层面一切正常熵值分布合理、权重和为 1只有落到具体点位上才看出方向反了。顺带说一个经验参照点要挑在影像中心区域、远离边缘和云影的地方。边缘像元受重采样和掩膜影响大值本身就不稳拿它做检验容易误判。最后再分享一个我在实际项目里养成的小习惯每次算完熵权把权重表、样本量、有效像元比例、各波段归一化前后的 min/max/分位数全部写进一个 JSON 文件和成果图一起归档。当时觉得多余但半年后回头复现同一套流程、或者要向别人解释为什么这个波段权重是 0.08时这份记录能省掉大量的重跑和回忆。数据处理的信任成本往往就藏在这些不起眼的元数据里。 
长河君 · 2026/10/26 20:31:17 · 阅读 1.2k 多波段TIFF熵权法赋权与GeoTIFF综合评价实践
手里攒着一批多波段 TIFF想把它们合成一张综合评价图最卡人的地方往往不是代码而是权重——可见光波段给 0.3 还是 0.25植被指数和建筑指数谁更重要凭经验拍脑袋定完评审一问依据是什么就哑了。熵权法处理 TIFF 图像这条路我用了几年核心价值就在于把权重从主观判断变成由数据自身离散程度推出来的客观量某个波段在所有像元上差异越大它携带的区分信息越多权重就越高。这篇内容讲的是一套完整可复现的流程——从多波段 TIFF 的读取清洗到熵权法的逐像元改造再到把权重图写回带坐标的 GeoTIFF中间那些让结果悄悄跑偏的坑我也会一个个拆开。适合手上有栅格数据、做过基础 numpy 操作、想把叠图分析做得更讲道理的人看也适合刚接触多指标评价、还在纠结权重怎么定的朋友。1. 多波段 TIFF 的赋权困境与熵权法的切入点1.1 叠图分析里最难交代的那一步多波段 TIFF 最常见的用法是做加权叠置把若干指标波段各自归一化乘上权重相加成一张综合得分栅格。这套流程本身不难rasterio 读进来、numpy 算一算、再写出去几十行就能跑通。真正让人反复返工的是权重表。我见过太多项目的权重是专家打分 微调凑结果先按经验给一组数跑出来发现某个区域得分偏低不合理回头把对应波段权重调高一点再跑一遍直到图看起来顺眼。这种做法当然能出图但它有个致命问题——结论是被结果反向塑造的你无法证明这张图反映的是数据本身的结构还是你心里的预期。更麻烦的是波段数量一多主观定权就彻底失控。五六个波段还能凭直觉排个序十几个波段比如多光谱加几个纹理特征、几个地形因子时人脑根本没法同时权衡这么多维度的重要性最后往往退化成平均分配或者挑几个重要的给大权重剩下的象征性给一点。平均分配听着公平实际上等于宣布所有波段信息量相同这几乎不可能成立——一个几乎全图恒定的波段和一个剧烈起伏的波段凭什么同等对待。熵权法切入的正是这个位置。它不问你这个指标重不重要只问这个指标在样本间的分布有多不均匀。分布越集中说明它提供不了区分度权重自然压低分布越分散说明它能把样本拉开权重自然抬高。整个推导过程是纯数学的任何人拿同样的数据都能算出同样的权重可复现、可审计。这一点在需要写报告、需要向别人解释为什么是这个权重的场合价值远大于精度上的那点差异。1.2 熵权法在多波段影像上的独特适配性熵权法原本是评价领域的通用方法用在影像上之所以顺手是因为两者结构天然同构。标准熵权法要求一个 n×m 的决策矩阵n 是样本数、m 是指标数一张 m 个波段的 TIFF本质上就是一个像元数 × 波段数的巨大矩阵。每个像元是一个待评价对象每个波段是一个评价指标这套映射几乎不需要额外解释。由此衍生出两种用法差别很大选错了结果会完全不同。用法样本粒度权重形态适用场景计算量全局熵权全图所有有效像元一组标量权重m 个整幅影像统一评价、区域横向对比中等可抽样估计局部熵权每个滑动窗口内的像元一张权重栅格m 层空间异质性强、需要自适应高随窗口增大线性上升全局熵权的逻辑是整幅图共用一套价值观好处是不同区域、不同时相的结果可以直接比较坏处是它假设指标的重要性在空间上不变。局部熵权则让每个窗口自己算一套权重能捕捉城市区域建筑指数更重要、山区植被指数更重要这类空间分异代价是计算量暴涨而且窗口内样本少时权重会剧烈抖动。我的经验是做区域间横向对比、做时间序列分析用全局熵权把权重固定死做单幅图的内部结构刻画、做异常区识别用局部熵权。两者甚至可以叠起来用——先用全局权重出一版基准图再用局部权重看哪些地方和基准偏离最大偏离本身就是一个很有价值的信息层。1.3 这个方法的能力边界要先说清楚熵权法不是万能的用之前得认清它几个固有短板否则容易过度解读。第一它只反映数据内部的离散程度不反映指标的实际意义。假设某个波段里混入了传感器噪声噪声让数值剧烈波动熵权法反而会给它很高的权重——因为离散度大。这是最典型的翻车场景所以前置的噪声抑制和异常值处理绝对不能省。第二它对指标间的相关性完全不敏感。三个高度相关的波段会被各给一份权重等于这个维度被重复计权三次。如果波段间相关系数普遍在 0.9 以上要么先做主成分降维要么手动合并不能直接扔进去。第三权重是相对的、依赖样本集的。换一批像元、换一个研究区权重就会变。所以权重表必须和数据集绑定说明不能当成放之四海皆准的常数。把这三点想在前面后面很多争论根本不会发生。2. 拆开熵权法的计算链条2.1 指标方向统一先解决越大越好还是越小越好熵权法第一步是归一化而归一化的公式取决于指标方向。这一步搞错后面全废——我在第 5 节会专门讲这个坑有多隐蔽。多波段影像里三类方向的波段都可能出现正向指标越大越好植被指数、水体指数、亮度、纹理能量。归一化用(x - min) / (max - min)。负向指标越小越好地表温度做热岛评价时、坡度做建设适宜性时、建筑密度做生态评价时。归一化用(max - x) / (max - min)。区间型指标落在某个区间最好比如适宜坡度 0~8 度太陡不行、太平也可能排水不畅。这类要用分段公式先设定理想区间[a, b]落在区间内给满分超出后按距离衰减。实际操作里我习惯在脚本开头维护一张波段-方向映射表用字典写死避免中途忘掉。这张表同时充当文档别人接手一眼能看懂每个波段的处理逻辑。# 波段索引 - (方向, 参数) BAND_RULES { 0: (positive, None), # NDVI 类 1: (positive, None), # 纹理能量 2: (negative, None), # 地表温度 3: (interval, (0.0, 8.0)), # 坡度理想区间 0~8 度 }归一化之后还有个细节会出现 0。熵权法下一步要取对数ln(0)直接炸。标准做法是整体平移一个极小量业界常用x 0.999 * x 0.001或者x x 1e-6。前者把数据压缩到 [0.001, 1] 区间平移量固定不同批次之间可比性更好后者更轻但极小值的影响在样本量大时会被放大。我倾向用前者理由是它的扰动是有界的、可预期的。2.2 比重矩阵、信息熵、差异系数三步推导归一化完成后进入核心三步。假设有 n 个有效像元、m 个波段矩阵记作Xn×m元素x_ij。第一步算比重。对每一列每个波段单独归一化到占比p_ij x_ij / Σ_i x_ij这一步的含义是在这个波段的总量里第 i 个像元贡献了多少份额。所有像元的份额加起来等于 1。第二步算信息熵。经典香农熵公式e_j -k * Σ_i (p_ij * ln(p_ij))其中 k 1 / ln(n)系数k的作用是把熵值标准化到 [0, 1] 区间这样不同样本量下的熵可以横向比较。如果所有像元在该波段上取值完全相同p_ij全为1/n代入后e_j 1熵取最大值说明这个波段毫无区分能力。第三步算差异系数和权重。d_j 1 - e_j w_j d_j / Σ_j d_jd_j是差异系数也叫信息效用值熵越小它越大。最后按列求和归一化得到一组和为 1 的权重。整个链条只有加、乘、对数三种运算没有任何迭代和优化所以结果完全确定——这也是它比某些机器学习赋权方法更适合做报告的原因你能把公式贴在附录里任何人都能复算。2.3 手算一遍三个像元两个波段公式看着抽象用最小规模过一遍就清楚了。设两个波段、三个像元已完成归一化且已加 0.001 平移像元1: [0.801, 0.101] 像元2: [0.501, 0.501] 像元3: [0.101, 0.901]波段 1 的列和为 1.403比重为[0.571, 0.357, 0.072]。代入熵公式ln(3) 1.0986k 0.9102Σ p*ln(p) 0.571*(-0.560) 0.357*(-1.030) 0.072*(-2.632) -0.320 - 0.368 - 0.190 -0.878 e_1 -0.9102 * (-0.878) 0.799 d_1 1 - 0.799 0.201波段 2 的列和为 1.503比重为[0.067, 0.333, 0.600]对称地算出来e_2 0.799, d_2 0.201两个波段差异系数相同权重各 0.5。这个例子故意造成对称说明一件事熵权法衡量的是分布形状不是数值大小。波段 1 的绝对值普遍偏小、波段 2 普遍偏大但两者的相对分布结构一样权重就一样。很多人第一次用会误以为数值大的波段权重大这是典型误解必须先纠正。再换一组数据让波段 2 变成几乎恒定[0.501, 0.502, 0.503]。它的比重接近[0.333, 0.334, 0.333]熵逼近 1差异系数趋近 0权重会被压到几乎为零。这就是熵权法自动淘汰无用指标的直观体现。2.4 从标量权重推广到逐像元权重全局权重算完只是半程。如果要做局部熵权需要把上面的流程搬进每个滑动窗口。窗口大小k×k内样本数是k²扣除无效值后可能更少指标数还是波段数 m。每个窗口算出一组权重把它赋给窗口中心像元滑完一遍就得到 m 层权重栅格。这里有个隐含约束常被忽略窗口内的有效样本数必须显著大于波段数。熵公式里的k 1/ln(n)当 n 很小时熵的估计极不稳定。n2 时ln(2)0.693两个点就能算出熵但完全没有统计意义。我的经验阈值是n ≥ 5m即窗口内有效像元至少是波段数的五倍若波段有 6 个窗口至少 30 个有效像元对应 7×7 窗口去掉边缘和无效值后大概能剩 30 左右。窗口再小权重图会变成一片噪声看着花哨实际不能用。3. 把 TIFF 读进来并清洗干净3.1 三条读取路线怎么选Python 生态里读 TIFF 主要有三条路各有适用边界选错了后期会很别扭。库优势局限适用tifffile纯读写速度快支持 BigTIFF、多页不解析地理坐标纯图像处理、无空间参考rasterio完整地理元数据、分块读写、重投影依赖 GDAL安装稍重遥感栅格、需要坐标PIL/imageio轻量适合 8 位 RGB多波段支持弱精度受限普通图片做综合评价一般都要保留空间参考因为最终结果要叠到地图上用所以rasterio是默认选择。它的open()返回的对象带profile里面存着坐标系、仿射变换、数据类型、分块大小最后写结果时直接复用这份 profile能省掉一大堆坐标对齐的麻烦。import rasterio import numpy as np with rasterio.open(input_stack.tif) as src: profile src.profile.copy() arr src.read() # 形状 (bands, rows, cols) nodata src.nodata print(arr.shape, arr.dtype, nodata)有个细节src.read()会把整幅图读进内存。一张 10000×10000、6 波段、float32 的影像约 2.4 GB勉强能塞进内存但要是 30000×30000直接就是 21 GB机器扛不住。这种情况必须走分块。3.2 NoData 与异常值最容易被跳过的一步原始 TIFF 里几乎必然存在无效区域影像边缘的黑边、云遮挡、水体掩膜、重采样产生的空洞。它们的值可能是-9999、0、nan也可能是某个奇怪的哨兵值。这些值如果不处理后果有两个层面。表层是数值污染——一个-9999塞进归一化min直接被拉到-9999整幅图的动态范围被压缩到几乎为零所有像元挤在一起熵权法算出来的权重毫无意义。深层是统计污染——即便你用了稳健的 min/max无效像元进入比重矩阵后会把某一列的总和撑大其他像元的份额被稀释权重计算整体失真。处理流程我一般这么走按src.nodata生成第一层掩膜叠加物理合理范围掩膜比如反射率应在[-0.1, 1.2]温度应在[-50, 70]超出的一律判无效叠加np.isfinite过滤 NaN 和 Inf三个掩膜取交集得到最终的valid_mask。valid np.ones(arr.shape[1:], dtypebool) if nodata is not None: valid ~np.all(arr nodata, axis0) valid ~np.any(arr nodata, axis0) valid np.all(np.isfinite(arr), axis0) valid np.all((arr -0.1) (arr 1.2), axis0) # 按实际物理量调整注意这里用的是np.any和np.all的区别只要任何一个波段无效这个像元就整体剔除。理由很实在——综合评价要求每个像元在各维度上都有值缺一个维度就没法算加权和与其插补引入人为误差不如直接标为无效最后在成果图上留白反倒更诚实。3.3 大文件的分块策略rasterio的block_windows()会按文件内部的分块组织tiling返回窗口这个分块在写入时就是按原始布局生成的读取效率最高。对只读需求也可以自己按固定尺寸切from rasterio.windows import Window BLOCK 2048 with rasterio.open(input_stack.tif) as src: for row_off in range(0, src.height, BLOCK): for col_off in range(0, src.width, BLOCK): win Window(col_off, row_off, min(BLOCK, src.width - col_off), min(BLOCK, src.height - row_off)) block src.read(windowwin) # ... 处理 ...分块有两个坑必须处理。第一个是边界效应如果做局部熵权块边缘的窗口需要邻块数据否则窗口不完整权重会在块与块之间出现明显断层。解决办法是读取时外扩half_win像素的重叠区算完再裁掉。第二个是全局统计量依赖全局熵权需要全图的min/max和列和这些必须先在完整数据上统计出来不能边分块边累加除非你严格按可加性设计比如只累加列和但 min/max 需要单独一遍扫描。所以标准做法是两遍扫描第一遍只读元信息算统计量第二遍才正式处理。3.4 数据类型与精度取舍原始数据可能是uint16如 0-10000 的反射率放大值而归一化后必然是浮点。中间过程的 dtype 选择会影响结果精度也影响内存。我的习惯是读取后立刻转float32全流程用float32只有熵值求和这种累积量用float64累加器。理由很实际——float32只有约 7 位有效数字样本量上百万时Σ p*ln(p)的累加误差可能到 1e-5 量级虽然对权重影响不大但如果直接在整个大数组上做np.sum误差会累积得更明显。numpy 的np.sum对float32输入默认用成对求和pairwise summation精度已经不错但遇到关键结果我还是会显式.astype(np.float64).sum()。写回时反过来如果成果只是给人看的得分图float32足够如果要参与后续定量分析用float32存 GeoTIFF 并设置compressdeflate体积能压到原始的三分之一左右。4. 逐像元熵权的向量化写法4.1 双层循环为什么必须干掉最直白的实现是两层 for 循环遍历每个像元把该像元所有波段的值取出来当成一个样本算熵权——但这是错的。熵权法的样本维度是参与比较的对象集合单个像元只有一行数据n1ln(1)0公式直接退化。所以逐像元的正确理解是窗口内所有像元构成样本集窗口中心像元获得这套权重下的综合得分。按这个理解写朴素实现就是外层遍历行、内层遍历列每个位置切片取窗口算一遍完整熵权流程。一张 2000×2000 的图就是 400 万次窗口运算每次内部还有若干次 numpy 小数组操作Python 解释器开销会把它拖到几十分钟甚至几小时。这个量级必须向量化。4.2 用滑动窗口视图替代切片numpy.lib.stride_tricks.sliding_window_view是关键工具。它不复制数据只通过调整 stride 构造一个视图形状从(H, W, C)变成(H-k1, W-k1, k, k, C)。之后所有运算都是整块数组操作走的是底层 C 循环。from numpy.lib.stride_tricks import sliding_window_view def entropy_weights_local(arr, valid, win7, eps1e-3): arr: (C, H, W) float32已归一化到 [0,1] valid: (H, W) bool win: 窗口边长奇数 return:(C, H, W) 权重图 (H, W) 综合得分 C, H, W arr.shape pad win // 2 # 边缘用 0 填充后续用 valid 掩膜排除 a np.pad(arr, ((0, 0), (pad, pad), (pad, pad)), modeconstant) v np.pad(valid, ((pad, pad), (pad, pad)), modeconstant) # (H, W, k, k, C) win_a sliding_window_view(a, (win, win), axis(1, 2)) win_v sliding_window_view(v, (win, win), axis(0, 1)) win_a np.moveaxis(win_a, -1, -2) # (H, W, C, k, k) win_v win_v[..., None, :, :] # (H, W, 1, k, k) n win_v.sum(axis(-2, -1)) # (H, W, 1) n_safe np.maximum(n, 2) x np.where(win_v, win_a, 0.0) col_sum x.sum(axis(-2, -1), keepdimsTrue) # (H, W, C, 1, 1) col_sum np.where(col_sum 1e-12, 1e-12, col_sum) p np.where(win_v, x / col_sum, 0.0) plogp np.where(p 0, p * np.log(np.maximum(p, 1e-12)), 0.0) e -plogp.sum(axis(-2, -1)) / np.log(n_safe) # (H, W, C) d 1.0 - e d np.where(d 0, 0.0, d) # 数值误差兜底 d_sum d.sum(axis-1, keepdimsTrue) d_sum np.where(d_sum 1e-12, 1e-12, d_sum) w d / d_sum # (H, W, C) # 综合得分中心像元值 × 权重对波段求和 center np.moveaxis(win_a, -1, -2)[..., win // 2, win // 2] # (H, W, C) center np.where(win_v[:, :, :, win // 2, win // 2], center, 0.0) score (center * w).sum(axis-1) # (H, W) w_out np.moveaxis(w, -1, 0) # (C, H, W) return w_out.astype(np.float32), score.astype(np.float32)这段代码有几个刻意的设计值得说明。np.where(win_v, win_a, 0.0)把所有无效像元在运算前清零这样列求和只会统计窗口内的有效样本不需要在公式里额外加分支。n_safe np.maximum(n, 2)是为了防止窗口内一个有效像元都没有时除以ln(0)或ln(1)这种情况下算出来的权重其实是垃圾值最后会由valid掩膜统一剔除所以没必要在计算中做精细判断只要不产生 inf/nan 污染邻居就行。d np.where(d 0, 0.0, d)这行看着多余其实很有必要。理论上e_j ∈ [0,1]但浮点误差可能让e_j略微超过 1d变成 -1e-8之后归一化时这个微小负值会被放大成奇怪的权重分布。加个下限钳位就能彻底杜绝。center那段的取值逻辑中心像元在窗口切片里的位置是(win//2, win//2)前面moveaxis之后维度顺序是(H, W, C, k, k)所以切片是[..., pad, pad]。这里容易写错索引调试时建议先用一个 5×5 全 1 数组验证一下中心取值是否对得上。4.3 内存与速度的实测取舍上面这段在 2000×2000、6 波段、7×7 窗口下内存峰值约 3~4 GB——因为sliding_window_view虽然是视图但后面np.where、除法、log都会生成同尺寸的新数组每个中间结果都是H×W×C×k×k的float32。7×7 就是 49 倍膨胀600 万像元 × 6 波段 × 49 × 4 字节 ≈ 7 GB非常吃内存。三个优化手段按性价比排序降窗口尺寸。从 15×15 降到 7×7内存直接降 4.6 倍。5×5 往往就够用前提是波段数不超过 3~4 个且满足n ≥ 5m。分块 重叠。按 512×512 分块每块外扩pad像素算完裁剪。内存峰值只跟块大小有关与整图无关。减少中间变量。用np.einsum替代显式乘加或者把p*log(p)写成原地操作能省一两个临时数组但收益有限且代码可读性下降除非确实卡内存否则不建议。速度方面纯 numpy 版本在 8 核机器上跑 2000×2000 大约 20~40 秒可接受。如果追求更快numba的njit(parallelTrue)能再快 5~10 倍但要把循环写回去调试成本上升。我的建议是先跑通 numpy 版本验证结果正确确认逻辑无误后再考虑用 numba 加速顺序别反。4.4 写回 GeoTIFF 时容易丢的东西结果写回最省事的写法是复用输入 profileout_profile profile.copy() out_profile.update( count1, dtypefloat32, nodata-9999.0, compressdeflate, predictor3, tiledTrue, blockxsize256, blockysize256, ) score np.where(valid, score, -9999.0).astype(float32) with rasterio.open(entropy_score.tif, driverGTiff, **{k: v for k, v in out_profile.items() if k ! driver}) as dst: dst.write(score, 1)几个必须检查的点坐标不能丢。如果只复制了crs和transform却忘了height/width或者分块处理时没做坐标偏移成果叠到地图上会整体偏移。分块写回时要用dst.write(block, 1, windowwin)指定窗口不能让 rasterio 自己去猜。predictor3只对浮点有效。predictor2是给整型的混用会让压缩失效甚至数据错乱。NoData 值要一致。输出声明nodata-9999实际数组里也要真的填-9999不能只在元数据里写、数据里留着 NaN——很多 GIS 软件会因此把整幅图判为无效。权重图建议单独存。m 层权重图对分析很有价值能看出哪些区域哪个波段占主导别只存综合得分把中间过程扔了。5. 五个让结果白算的坑5.1 方向搞反结果看起来很有道理但完全颠倒这是最隐蔽的坑因为出图不会报错甚至配色出来还挺好看。我在一个热环境评价的项目里踩过地表温度波段本该是负向指标写规则表时手滑写成正向最终成果图显示中心城区环境质量最优——因为温度高的地方被当成高分。图看着平滑、对比度也不错差点就这么交出去了是同事顺口问了句市中心怎么是最优区才发现。排查链路其实很简单但要在出图前主动做打印每个波段的归一化结果直方图正向指标和负向指标的直方图应该互为镜像抽 10 个已知类别的样本点比如确定的城区点、林地五点看它们的归一化值是否与预期方向一致计算每个波段与综合得分的相关系数如果某个波段相关系数为负而它被标为正向立刻报警。第 3 条我后来做成了脚本里的断言跑完自动检查省得靠肉眼。for j in range(C): r np.corrcoef(norm[j][valid], score[valid])[0, 1] direction BAND_RULES[j][0] if direction positive and r -0.05: print(f[警告] 波段{j} 相关系数 {r:.3f}方向可能标反) if direction negative and r 0.05: print(f[警告] 波段{j} 相关系数 {r:.3f}方向可能标反)5.2 常数列与 log(0)两个都会让权重凭空消失常数列的典型来源是掩膜波段、全图同一值的分类波段、或者某个波段在裁剪后只剩下一个取值。它的熵为 1、差异系数为 0本身处理起来没问题权重为 0 是合理的但麻烦在于它会让Σ d_j变小其他波段的权重被按比例放大——如果只有一个波段有效而其他全是常数列那个波段权重直接等于 1综合得分就等于它自己看着合理其实毫无信息。log(0)的来源是归一化后出现精确的 0。虽然前面加了 0.001 平移但如果某些像元的值恰好等于列最小值x会是 0.001 而非 0通常安全。真正的风险在别处col_sum为 0。如果某个窗口内该波段全为 0比如有效像元都被掩膜掉了窗口内只剩几个有效点而它们的值恰好为 0p会变成 0/0 nannan 会扩散。防御方式是三层col_sum做下限钳位、np.log前用np.maximum(p, 1e-12)、算完统一np.nan_to_num。这三层看着冗余但真跑大图时能救命——我曾因为漏了第三层nan 从一个窗口扩散到整张权重图最后输出全黑排查了两个小时。5.3 归一化被离群值带偏熵权法对极值极其敏感因为归一化的分母是max - min。一个像元上的异常高值传感器坏点、镜面反射、火点就能把整个波段的动态范围拉大导致其余 99.9% 的像元被压缩到很窄的区间熵值偏高、权重偏低这个波段等于被废掉。处理手段有三种按侵入性递增方案做法优点风险分位数截断用 2% 和 98% 分位代替 min/max简单效果明显超出部分被钳到边界损失极端信息稳健标准化用中位数和 IQR 替代均值方差抗离群会改变分布形态熵值含义略变先滤波再算中值滤波或多时相合成去异常从源头解决引入额外处理步骤可能损失细节我一般用分位数截断2%~98% 或 1%~99% 视数据质量定。判据是截断后重新看直方图如果分布从一根尖峰 长尾变成近似对称的钟形说明异常值确实被压住了如果截断后分布明显被削平说明分位点选得太激进要往回收。有个容易忽略的细节分位数必须在全图有效像元上统计不能分块各自统计否则不同块的归一化基准不一致拼起来会有明显的块状接缝。这个现象我第一次遇到时以为是分块代码写错了检查了半天才发现是统计口径的问题。5.4 窗口尺度和波段数量不匹配局部熵权的窗口大小不是随便定的它受两个约束同时限制。下限来自样本量窗口内有效像元数要满足n ≥ 5m。6 波段的最小窗口约 7×749 个像元扣除无效后通常剩 35~45。如果强行用 3×3只有 9 个样本算出来的权重会在相邻像元之间剧烈跳变权重图看起来像随机噪声。上限来自空间分辨率窗口不能大到抹平你要观察的空间结构。做城市内部的功能区识别几十米分辨率下窗口超过 30×30约 1 公里就会把街区尺度差异全抹平权重图变成一片均匀的色块失去局部自适应的意义。中间的取舍只能靠试。我的做法是先固定波段数跑 5×5、9×9、15×15 三组看综合得分的空间变异系数标准差/均值随窗口的变化曲线曲线从陡降转为平缓的那个拐点通常就是合适尺度。这个方法不算严谨但比拍脑袋强。还有个实操细节窗口边长用奇数。偶数边长时中心像元没有唯一定义不同库的默认处理不一致容易出现半像素偏移。5.5 权重图没人看得懂呈现方式的问题算完 m 层权重图如果直接扔给别人多半得到的反馈是这花花绿绿的是什么。权重图的解读需要引导。有效做法是生成一张主导波段图对每个像元取权重最大的波段索引用分类配色显示。这样一眼能看出东部这片全是第 2 波段主导西部是第 4 波段主导直接对应到空间格局上。再配一张权重堆叠面积图随机抽 5000 个点按行堆叠展示各波段权重的整体占比和波动范围。综合得分的呈现也有讲究。原始得分范围通常在 [0.1, 0.9]直接线性拉伸会显得灰蒙蒙。我的处理是先按 2%~98% 分位做拉伸再用自然断点法Jenks分 5 级配色这样既保留了对比度分级又符合数据的实际聚集结构。分级数不要超过 7 级人眼分辨不了那么多。配色上尽量避开红绿同时高饱和度出现的方案色觉障碍人群看起来会混淆。用蓝-黄、紫-橙这类对立色对更稳妥。6. 结果可信度的交叉验证6.1 权重合理性看三个指标熵权算完先别急着出图花十分钟做三项自检。第一权重排序是否符合常识。如果一个明显充满噪声的波段权重排第一说明噪声没压住如果一个信息量很足的波段权重接近 0说明它可能被归一化压扁了。这里的常识不是主观偏好而是基于数据质量的判断——你能从直方图、异常值比例上找到依据。第二权重分布的集中度。计算权重的变异系数。如果所有波段权重都在1/m附近变异系数 0.1说明各波段离散程度接近熵权法没提供多少信息量此时用等权法反而更稳、更好解释。如果某个波段权重超过 0.5要特别警惕检查是不是其他波段的归一化出了问题。第三敏感性。把权重上下扰动 ±10%重新算综合得分看排序前 10% 的像元集合变化多大。变化很小说明结果稳健变化剧烈说明结果对权重高度敏感此时权重的小误差会放大成结论的大误差需要慎重表述。6.2 和等权法、PCA 做对照单一方法的结论总是让人不放心我习惯同时算三个版本做对照。等权法是最简单的baselinew_j 1/m。它假设所有波段同等重要虽然粗暴但胜在无偏。主成分法取第一主成分的载荷绝对值归一化作为权重它捕捉的是方差最大的方向跟熵权法思路不同但目标相近两者结果接近时可信度大幅提升。相关系数法用各波段与综合得分的相关系数归一化作为权重属于结果导向的赋权适合做验证而不适合做定权。三者的排序一致性用 Spearman 秩相关衡量。我实测的经验值是熵权与 PCA 的秩相关在 0.85 以上说明数据的主结构清晰结论可靠低于 0.6 就要回头查数据通常是波段间共线性太强或者方向标错了。from scipy.stats import spearmanr r1, _ spearmanr(score_entropy[valid], score_equal[valid]) r2, _ spearmanr(score_entropy[valid], score_pca[valid]) print(f熵权 vs 等权: {r1:.3f}, 熵权 vs PCA: {r2:.3f})注意这个比较要在大样本上做抽 10 万个有效像元就够别在全图上跑spearmanr那个函数对百万级样本很慢。6.3 用已知参照点做落地检验纯统计验证之外最有说服力的是拿真实参照点检验。做法是找一批已知属性的位置——比如几处确认的高质量绿地、几处确认的重度扰动区、几处确认的水体——看它们的综合得分是否落在预期区间。这个检验不需要很多点二三十个就足以暴露系统性问题。如果高质量区的得分中位数低于扰动区无论统计指标多漂亮都得推翻重来。我在一个项目里就是靠这个检验发现了负向指标漏标的问题统计层面一切正常熵值分布合理、权重和为 1只有落到具体点位上才看出方向反了。顺带说一个经验参照点要挑在影像中心区域、远离边缘和云影的地方。边缘像元受重采样和掩膜影响大值本身就不稳拿它做检验容易误判。最后再分享一个我在实际项目里养成的小习惯每次算完熵权把权重表、样本量、有效像元比例、各波段归一化前后的 min/max/分位数全部写进一个 JSON 文件和成果图一起归档。当时觉得多余但半年后回头复现同一套流程、或者要向别人解释为什么这个波段权重是 0.08时这份记录能省掉大量的重跑和回忆。数据处理的信任成本往往就藏在这些不起眼的元数据里。