1. 为什么你的高光谱图像颜色总是不对？
我猜很多刚开始玩高光谱图像的朋友都踩过这个坑：辛辛苦苦拿到了一堆波段数据，想把它变成一张人眼能看的彩色图，结果合成出来的颜色要么灰蒙蒙的，要么绿得发慌，跟真实场景完全对不上。你可能试过直接挑三个波段（比如红、绿、蓝附近的波段）拼成RGB，也试过用一些图像处理软件自带的“伪彩色”功能，但效果总是不尽人意，感觉颜色“假假的”。
这背后的核心问题在于，我们人眼看到的“颜色”，并不是光谱信息的简单堆砌。我们的眼睛里有三种视锥细胞，分别对长波（偏红）、中波（偏绿）、短波（偏蓝）的光敏感。但它们的敏感曲线并不是方方正正的，而是有重叠、有起伏的复杂函数。一台高光谱相机能记录几十甚至几百个窄波段的光强，这远远超过了人眼三种细胞的“分辨率”。所以，把高光谱变成真彩色RGB，本质上是一个“降维”过程：如何把上百个波段的信息，精准地“压缩”成三个值（R, G, B），并且让这三个值刺激我们眼睛后产生的颜色感觉，和真实世界看到的感觉一致。
这就是CIE XYZ色彩空间出场的时候了。你可以把它想象成一个绝对客观的“颜色翻译官”。它不依赖于任何具体的设备（比如你的显示器或打印机），而是基于大量人类视觉实验数据建立的一套标准。它的核心工具就是色匹配函数。这套函数描述了对于每一个波长的光，需要多少“虚拟”的X、Y、Z三原色光来匹配，才能让人眼觉得“颜色一样”。通过CIE XYZ这个桥梁，我们就能把物理的光谱数据，翻译成符合人眼感知的颜色信号，然后再转换到我们显示器用的RGB空间。这个方法不是“大概齐”，而是一种科学、可复现的精准合成方法，在遥感、医学影像、艺术品数字化等领域，对颜色保真度要求高的场景里特别有用。
2. 理解核心：色匹配函数与CIE XYZ色彩空间
2.1 色匹配函数：人眼颜色的“密码本”
要理解色匹配函数，我们得从人眼的老底说起。咱们眼睛里负责看颜色的视锥细胞有三种，就像三款不同的光传感器。它们各有偏好：L型喜欢长波红光，M型喜欢中波绿光，S型喜欢短波蓝光。但它们的“喜好曲线”是平滑且相互重叠的，并不是只认单一波长。
上世纪二三十年代，科学家们做了一系列非常精巧的实验来量化这种感知。实验大概是这样：让被试者看一个分成两半的屏幕。左边打上一束特定波长的单色光，比如一道纯净的555nm黄光。右边则用红、绿、蓝三束可以调节强度的光同时照射。然后让被试者调节右边三色光的强度，直到他觉得左右两边的颜色看起来一模一样为止。记录下此时红、绿、蓝光的强度值（r, g, b）。这个（r, g, b）就是匹配这一波长单色光所需的“配方”。
科学家们把可见光范围内的波长从头到尾测了个遍，就得到了三组数据，这就是最初的CIE RGB色匹配函数。但这里出现了一个有趣的现象：匹配某些颜色时，比如某些蓝绿色，即使用尽右边所有的红光，颜色还是对不上。怎么办？实验者想了个办法：把一点红光加到左边的被测光里，然后再调右边。这在数学上就等价于在右边的配方里，给红光强度赋了一个负值。所以最初的CIE RGB色匹配函数曲线是包含负值区域的。
负值在计算上很麻烦。所以，科学家们通过一个数学上的线性变换，将整个空间变换到了一个所有分量都为正值的新空间。这个新空间就是CIE XYZ色彩空间，而变换后的三组数据，就是我们现在常用的 CIE 1931 XYZ标准色匹配函数，记作 x̄(λ), ȳ(λ), z̄(λ)。这个变换是纯数学的，XYZ是三个假想的、不存在于现实中的原色，但它们构建的空间非常完美，涵盖了所有人眼能见的颜色，且计算方便。
2.2 获取与使用色匹配函数数据
在实际编程中，我们不需要自己去做实验。色匹配函数有标准化的数据表。你可以从国际照明委员会（CIE）的官网、维基百科或一些权威的色彩科学网站（如cvrl.org）下载到这些数据。数据通常是波长（单位纳米）和对应的x̄, ȳ, z̄三个值。
对于编程来说，我们有两种使用方式：
- 查表法：如果你的高光谱波段波长是离散且与标准数据表波长对齐或接近，可以直接用插值（如线性插值）获取每个波段对应的色匹配函数值。这是最直接的方法。
- 解析近似法：有时查表不方便，或者需要连续函数。色匹配函数可以用一组高斯函数的和来高精度近似。维基百科上给出了1931年标准的近似公式。我在实际项目里用过这个近似方法，效果和查表几乎没差别，而且代码写起来更干净。
下面这个Python函数就是用来计算近似值的，你可以直接拿去用：
import numpy as np
def piecewise_gaussian(x, u, sigma1, sigma2):
"""分段高斯函数，用于近似色匹配函数"""
return np.where(x < u,
np.exp(-0.5 * (x - u) ** 2 / sigma1 ** 2),
np.exp(-0.5 * (x - u) ** 2 / sigma2 ** 2))
def cie_1931_xyz_cmf(wavelengths):
"""输入波长数组（纳米），返回近似的CIE 1931 XYZ色匹配函数值"""
# 使用维基百科提供的参数进行近似
x_bar = (1.056 * piecewise_gaussian(wavelengths, 599.8, 37.9, 31.0) +
0.362 * piecewise_gaussian(wavelengths, 442.0, 16.0, 26.7) -
0.065 * piecewise_gaussian(wavelengths, 501.1, 20.4, 26.2))
y_bar = (0.821 * piecewise_gaussian(wavelengths, 568.8, 46.9, 40.5) +
0.286 * piecewise_gaussian(wavelengths, 530.9, 16.3, 31.1))
z_bar = (1.217 * piecewise_gaussian(wavelengths, 437.0, 11.8, 36.0) +
0.681 * piecewise_gaussian(wavelengths, 459.0, 26.0, 13.8))
return x_bar, y_bar, z_bar
你可以把计算出的三条曲线画出来看看，它们非常直观地反映了人眼三种颜色通道对不同波长光的相对敏感度。其中 ȳ(λ) 函数本身就被定义为“标准光度函数”，它直接代表了人眼对亮度的感知，这也是为什么后续计算中，Y分量常常被当作明度值来使用。
3. 从光谱到XYZ：计算三刺激值
拿到了色匹配函数这把“尺子”，我们就可以开始度量高光谱图像中每个像素的颜色了。假设我们有一个像素，它的高光谱数据是一个向量，记录了在N个波段上的强度值 I₁, I₂, ..., Iₙ，对应的波长是 λ₁, λ₂, ..., λₙ。
这个像素的颜色在CIE XYZ空间中的坐标（X, Y, Z），称为三刺激值。计算原理就是“加权求和”：每个波段的强度，乘以该波段对人眼X、Y、Z通道的刺激贡献（即色匹配函数值），然后对所有波段求和。在连续光谱下，这是一个积分；在我们的离散波段数据下，就是一个点乘运算。
公式很简单：
X = Σ [I(λ_i) * x̄(λ_i)]
Y = Σ [I(λ_i) * ȳ(λ_i)]
Z = Σ [I(λ_i) * z̄(λ_i)]
这里的求和（Σ）是对所有波段 i 进行。
用Python和NumPy实现起来非常高效，可以利用数组广播机制一次性计算整个图像：
def hyper_to_xyz(hyper_cube, wavelengths):
"""
将高光谱数据立方体转换为XYZ图像。
hyper_cube: numpy数组，形状为 (高度, 宽度, 波段数)
wavelengths: 一维数组，长度等于波段数，单位纳米
返回: XYZ图像，形状为 (高度, 宽度, 3)
"""
# 1. 获取色匹配函数权重
x_bar, y_bar, z_bar = cie_1931_xyz_cmf(wavelengths)
# 2. 计算每个像素的XYZ值
# 利用np.tensordot进行张量点乘，效率很高
X = np.tensordot(hyper_cube, x_bar, axes=([2], [0]))
Y = np.tensordot(hyper_cube, y_bar, axes=([2], [0]))
Z = np.tensordot(hyper_cube, z_bar, axes=([2], [0]))
# 堆叠成XYZ图像
xyz_image = np.stack([X, Y, Z], axis=-1)
return xyz_image
这段代码的核心是 np.tensordot，它沿着高光谱数据的波段维度（axis=2）和色匹配函数的波长维度进行点乘求和，直接得到每个像素的X、Y、Z值。这一步之后，我们就得到了一个XYZ空间的图像，它包含了符合人眼感知的颜色信息，但还不能直接显示在普通屏幕上。
注意：这里有一个关键点，你的高光谱数据
I(λ_i)最好是经过辐射定标后的物理量（如辐射亮度），或者至少是线性响应的原始DN值。如果数据已经被非线性压缩（比如做了伽马校正），需要先转换回线性空间，否则颜色计算会出错。这是很多新手容易忽略的地方。
4. 关键转换：从XYZ空间到你的RGB屏幕
现在我们有了一张XYZ图像，但我们的显示器、手机屏幕使用的是RGB色彩空间。不同的设备（sRGB, Adobe RGB, Display P3等）有不同的色域，即它们能显示的颜色范围。同时，它们用来混合出颜色的“红”、“绿”、“蓝”三原色，在CIE XYZ空间中有不同的坐标。因此，我们需要一个转换矩阵，将设备无关的XYZ值，映射到特定设备的RGB值。
这个转换是一个线性变换，可以用一个3x3的矩阵M来表示：
[R_lin] [X]
[G_lin] = M * [Y]
[B_lin] [Z]
这里的 R_lin, G_lin, B_lin 是线性RGB值。
矩阵M从哪里来？它由目标RGB色彩空间的三原色和白色点在XYZ空间中的坐标决定。幸运的是，我们不需要自己推导。色彩科学家Bruce Lindbloom维护了一个非常全面的网站，列出了各种标准RGB空间的转换矩阵。我们以最常用的sRGB色彩空间（D65白点）为例，其逆转换矩阵（从XYZ到线性sRGB）如下：
| 色彩空间 | 白色点 | 矩阵 M (XYZ -> linear sRGB) |
|---|---|---|
| sRGB | D65 | [[ 3.2404542, -1.5371385, -0.4985314], [-0.9692660, 1.8760108, 0.0415560], [ 0.0556434, -0.2040259, 1.0572252]] |
在代码中，我们这样实现转换：
def xyz_to_linear_srgb(xyz_image):
"""
将XYZ图像转换到线性sRGB空间。
xyz_image: numpy数组，形状为 (H, W, 3)，代表XYZ
返回: 线性sRGB图像，形状为 (H, W, 3)
"""
# sRGB (D65) 的逆矩阵，从XYZ到线性RGB
M_inv = np.array([[ 3.2404542, -1.5371385, -0.4985314],
[-0.9692660, 1.8760108, 0.0415560],
[ 0.0556434, -0.2040259, 1.0572252]])
# 重塑以便进行矩阵乘法
orig_shape = xyz_image.shape
xyz_flat = xyz_image.reshape(-1, 3)
rgb_lin_flat = np.dot(xyz_flat, M_inv.T) # 注意是点乘矩阵的转置
rgb_lin_image = rgb_lin_flat.reshape(orig_shape)
return rgb_lin_image
这里有一个细节：矩阵M通常定义为 RGB = M * XYZ，所以我们在代码里用 np.dot(xyz, M_inv.T)，其中 M_inv 是M的逆矩阵。上面表格给出的就是可以直接使用的逆矩阵。完成这一步，我们就得到了线性sRGB图像。但工作还没完，直接显示这个图像，颜色还是会显得暗淡、不准确。
5. 最后一步：伽马校正——让颜色“活”过来
显示器（以及人眼）对光强的感知是非线性的。为了在有限的数字编码范围（如0-255）内更好地分配数据位深，并补偿显示器的非线性响应，RGB数据在存储和传输前会进行一次伽马编码。最常见的标准是sRGB伽马曲线。
因此，我们的线性RGB值 R_lin, G_lin, B_lin 必须经过伽马校正，才能得到最终用于显示或存储的 R_sRGB, G_sRGB, B_sRGB。
sRGB的伽马校正不是简单的幂律函数，而是一段线性函数和一段幂函数的组合：
如果 C_lin <= 0.0031308:
C_sRGB = 12.92 * C_lin
否则:
C_sRGB = 1.055 * (C_lin ** (1/2.4)) - 0.055
其中C代表R、G、B中的任何一个通道。
这个分段函数在低亮度区是线性的，主要是为了避免数值精度问题。在要求不极端严格的场合，比如快速预览，你也可以用一个简单的 γ=2.2 的幂律函数来近似：C_sRGB = C_lin ** (1/2.2)。但在追求精准的合成中，建议使用标准sRGB公式。
实现代码如下：
def linear_srgb_to_srgb(rgb_lin_image):
"""
将线性sRGB图像进行伽马校正，转换为标准的sRGB图像。
rgb_lin_image: numpy数组，值域理论上应在[0, +∞)，实际计算后需要裁剪。
返回: sRGB图像，值域[0, 1]
"""
# 首先，将负值裁剪为0（超出色域的颜色转换后可能出现负值）
rgb_lin_clipped = np.maximum(rgb_lin_image, 0)
# 应用sRGB伽马校正
mask = rgb_lin_clipped <= 0.0031308
rgb_srgb = np.where(mask,
12.92 * rgb_lin_clipped,
1.055 * np.power(rgb_lin_clipped, 1/2.4) - 0.055)
# 将值域限制在[0, 1]之间，超过1的可能是高光过曝区域，裁剪为1
rgb_srgb = np.clip(rgb_srgb, 0, 1)
return rgb_srgb
经过伽马校正后，图像的颜色对比度会变得更加自然，亮部和暗部的细节也得到更好的呈现。这时候的RGB图像，才是真正意义上符合人眼感知、可以在标准sRGB显示器上正确显示的真彩色图像。
6. 实战整合与效果调优
让我们把上面的所有步骤串起来，写一个完整的、健壮的合成类。这个类会处理从波长信息加载到最终RGB图像生成的全过程，并且我加入了一些我在实际项目中总结的调优技巧。
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
class HyperspectralToRGBConverter:
def __init__(self, sensor_wavelengths, target_space='sRGB_D65'):
"""
初始化转换器。
sensor_wavelengths: 一维数组，传感器各波段的中心波长（单位：纳米）。
target_space: 目标RGB色彩空间，支持 'sRGB_D65'（默认）, 'AdobeRGB_D65'。
"""
self.wavelengths = sensor_wavelengths
self._load_color_matching_functions()
self._set_transformation_matrix(target_space)
def _load_color_matching_functions(self):
"""加载CIE 1931标准色匹配函数，并插值到传感器波长上。"""
# 这里使用标准数据点（示例，实际应使用更完整的数据）
# 波长范围约380nm到780nm，间隔1nm或5nm
std_wave = np.arange(380, 781, 1) # 假设有1nm间隔的标准数据
# 假设我们已经有了 std_x, std_y, std_z 标准值
# 可以从文件加载，这里用近似函数生成示例
std_x, std_y, std_z = cie_1931_xyz_cmf(std_wave)
# 将标准函数插值到传感器波长上
self.cmf_x = np.interp(self.wavelengths, std_wave, std_x)
self.cmf_y = np.interp(self.wavelengths, std_wave, std_y)
self.cmf_z = np.interp(self.wavelengths, std_wave, std_z)
def _set_transformation_matrix(self, target_space):
"""根据目标色彩空间设置XYZ->线性RGB的转换矩阵。"""
matrices = {
'sRGB_D65': np.array([
[3.2404542, -1.5371385, -0.4985314],
[-0.9692660, 1.8760108, 0.0415560],
[0.0556434, -0.2040259, 1.0572252]
]),
'AdobeRGB_D65': np.array([
[2.0413690, -0.5649464, -0.3446944],
[-0.9692660, 1.8760108, 0.0415560],
[0.0134474, -0.1183897, 1.0154096]
])
}
self.matrix = matrices.get(target_space, matrices['sRGB_D65'])
def convert(self, hyper_cube, normalize=True, clip_percentile=(0.5, 99.5)):
"""
主转换函数。
hyper_cube: 高光谱数据立方体，(H, W, C)。
normalize: 是否对结果RGB进行归一化拉伸以增强视觉效果。
clip_percentile: 归一化时使用的百分比裁剪，避免极端值影响。
返回: sRGB图像，值域[0, 1]，形状(H, W, 3)。
"""
# 1. 计算XYZ三刺激值
X = np.tensordot(hyper_cube, self.cmf_x, axes=([2], [0]))
Y = np.tensordot(hyper_cube, self.cmf_y, axes=([2], [0]))
Z = np.tensordot(hyper_cube, self.cmf_z, axes=([2], [0]))
# 堆叠并重塑以进行矩阵乘法
h, w = X.shape
xyz_flat = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
# 2. 转换到线性RGB
rgb_lin_flat = np.dot(xyz_flat, self.matrix.T)
rgb_lin = rgb_lin_flat.reshape(h, w, 3)
# 3. 伽马校正
rgb = self._gamma_correct(rgb_lin)
# 4. 可选：归一化拉伸以优化显示
if normalize:
rgb = self._percentile_normalize(rgb, clip_percentile)
return np.clip(rgb, 0, 1)
def _gamma_correct(self, rgb_lin):
"""应用sRGB伽马校正。"""
rgb_lin_clipped = np.maximum(rgb_lin, 0)
mask = rgb_lin_clipped <= 0.0031308
rgb = np.where(mask,
12.92 * rgb_lin_clipped,
1.055 * np.power(rgb_lin_clipped, 1/2.4) - 0.055)
return rgb
def _percentile_normalize(self, image, percentiles):
"""按百分比裁剪并归一化到[0,1]。用于增强图像对比度。"""
low, high = np.percentile(image, percentiles[0]), np.percentile(image, percentiles[1])
image_clipped = np.clip(image, low, high)
image_normalized = (image_clipped - low) / (high - low + 1e-10)
return image_normalized
# 使用示例
if __name__ == "__main__":
# 假设你已经有了高光谱数据 hyper_data (H, W, 128) 和波长列表 wavelengths_128
# converter = HyperspectralToRGBConverter(wavelengths_128, target_space='sRGB_D65')
# rgb_result = converter.convert(hyper_data, normalize=True, clip_percentile=(0.5, 99.5))
#
# plt.imshow(rgb_result)
# plt.axis('off')
# plt.show()
pass
效果调优经验分享：
- 归一化技巧：直接转换出来的RGB值动态范围可能很窄，导致图像对比度低。我常用的方法是使用百分比裁剪归一化（如代码所示），比如裁剪掉0.5%的最低值和0.5%的最高值，然后再拉伸到0-1。这能有效去除噪声极值，显著提升视觉观感，且不损失主要信息。
- 色域外颜色：有些高光谱颜色（比如非常饱和的单色光）可能超出了sRGB的色域范围，转换后RGB值会出现负值或大于1。我的处理是先在线性RGB阶段将负值裁剪为0，伽马校正后再将大于1的值裁剪为1。这相当于将超色域颜色映射到色域边界，是通用做法。
- 波段响应修正：如果你的高光谱相机各波段的相对响应度不一致（通常如此），最理想的情况是使用经过辐射定标和大气校正后的地表反射率数据。如果只有原始DN值，合成颜色可能会整体偏色。这时可以尝试在计算XYZ前，对高光谱数据应用一个基于白板或已知参考物的简单增益校正。
- 尝试不同RGB空间：如果你感觉sRGB合成的颜色不够鲜艳，可以尝试转换到色域更广的Adobe RGB（1998）空间。代码中提供了矩阵选项。但注意，最终显示设备要支持该色彩空间，否则颜色管理不当反而会出错。
踩过几次坑之后，我发现基于CIE XYZ的转换方法最大的优点是稳定和可解释。只要你的输入数据是线性的，波长信息准确，那么合成结果在颜色关系上就是科学准确的。它可能不会每次都产生最“惊艳”的视觉效果，但它提供的是一种可靠的、与物理世界和人眼感知对齐的颜色基准。在这个基础上，你再去做风格化的调整，心里就有底了。
147
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
