拿到一份 LED 光谱数据，怎么秒算色坐标、色温和显色指数？

> 关键词：LED、光谱分析、色度学、显色指数、Python、CIE1931

        做照明、做光源的人对下面这个场景一定不陌生：分光光度计导出来一份 `波长(nm) – 辐照度` 的数据，老板/客户问一句"这个灯色温多少？显色指数够不够？主波长在哪？"——如果你还在一个个公式手算，这篇文章能帮你把整套流程跑通，文末还放了一个我自研的小工具。

—

## 一、先搞懂：光谱功率分布（SPD）才是"源头"

        所有颜色指标，都是从一份 **光谱功率分布（Spectral Power Distribution, SPD）** 算出来的。SPD 就是"每个波长上有多少光能量"：

– 横轴：波长 λ，可见光 380–780 nm；
– 纵轴：相对辐照度 / 相对功率 S(λ)。

        比如一颗白光 LED，典型的 SPD 长这样：450 nm 附近有个尖峰（蓝光芯片），再叠一段宽宽的黄绿光（荧光粉转换）——蓝+黄绿混出"白"。

        后面所有指标（色坐标、色温、显色指数）都是对这条曲线做积分/运算得到的。**所以只要你有了 SPD，颜色就完全确定了。**

—

## 二、从 SPD 到色坐标（CIE 1931）

        人眼有三种视锥细胞，对应 CIE 1931 标准观察者的三条配色函数 x̄(λ)、ȳ(λ)、z̄(λ)。把 SPD 和它们分别相乘再积分，得到三刺激值：

```
X = ∫ S(λ)·x̄(λ) dλ
Y = ∫ S(λ)·ȳ(λ) dλ      ← Y 同时代表人眼亮度
Z = ∫ S(λ)·z̄(λ) dλ
```

归一化得到色坐标：

```
x = X / (X + Y + Z)
y = Y / (X + Y + Z)
```

        把 (x, y) 画到 **CIE 1931 色度图**上，就能直观看到这个光源"偏什么颜色"。下面是示例白光 LED 在色度图上的位置，以及黑体（普朗克）轨迹和 D65/D50 标准白点：


—

## 三、色温 CCT 与普朗克轨迹

        黑体被加热到不同温度会发出不同颜色的光，这条"黑体发光轨迹"就是**普朗克轨迹**（图2 橙色线）。一个光源的**相关色温 CCT**，就是"色度图上离它最近的黑体温度"。

        严格地讲，CCT 要在**黑体（普朗克）辐射谱**上求解。把普朗克定律对每个波长积分、再投影到 CIE 1931，就得到温度 T 对应的黑体色坐标 (x, y)——也就是普朗克轨迹上一点：

```
S(λ,T) = 1 / [ λ⁵ · (exp(c₂/(λT)) − 1) ]      c₂ = 0.0143877696 m·K，λ 取米
x(T), y(T) = ∫ S·x̄ dλ / Σ ,  y(T) = ∫ S·ȳ dλ / Σ    （Σ = ∫S·(x̄+ȳ+z̄)dλ）
```

        求解 CCT 就是在普朗克轨迹上找一点，使它与目标色坐标 (x, y) 的"距离"最小。**本文工具（也是第八节 demo）的真正做法**：在 **CIE 1960 UCS (u, v) 空间**（普朗克轨迹在这里近似为直线，比 (x,y) 更稳）里，用**牛顿法**在上面的精确 Planck 轨迹上直接迭代收敛到真实 CCT（中心差分取一阶/二阶数值导数，≤15 步，容差 1e-4 K）。为了加速收敛，先用 **McCamy 闭式公式**给一个毫秒级初值——它**只作"种子"，不参与最终精度**：

```
n = (x−0.3320)/(y−0.1858);  T_seed = −449·n³ + 3525·n² + 6823.3·n + 5520.33
```

> 注意：McCamy 单独用误差通常在几个百分点，只能当起点；本文工具最终 CCT 由牛顿法在**精确 Planck 积分**轨迹上收敛得到，并非 McCamy 直出。第八节完整 Python 按此复刻，可直接跑。

—

## 四、主波长与色纯度（饱和度）

– **主波长（Dominant Wavelength）**：从该光源颜色向光谱轨迹作延长线，交点对应的单色光波长。它描述颜色的"色调"（红/绿/蓝倾向）。
– **色纯度（Purity）**：`|白点→样品点| / |白点→光谱轨迹交点|`，描述颜色的"饱和程度"。

—

## 五、峰值波长、半峰宽（LED 最关心）

        对 LED 工程师来说，这两个比色温还重要：

– **峰值波长 λ_peak**：SPD 曲线最高点对应的波长（直接 `argmax` 即可）。
– **半峰宽 FWHM（Full Width at Half Maximum）**：曲线峰值一半处的宽度，反映光谱"窄不窄"——窄带 LED 颜色纯，宽带 LED 混光更白。

—

## 六、显色指数 CRI（重头戏）

        显色指数衡量"这个灯照在物体上，颜色还原准不准"。

– **CRI（CIE 13.3-1995）**：用 15 块标准色板（R1–R14，另加 R15 亚洲人肤色），分别算光源下与参考光源下的色差，得到 R1–R15，取前 8 块的平均就是**一般显色指数 Ra**。


—

## 七、光通量与光视效能

– **光通量 Φ（lm）**：`Φ = K_m · ∫ S(λ)·ȳ(λ) dλ`，`K_m = 683 lm/W` 是最大光谱光视效能。
– **光视效能（lm/W）**：光通量 ÷ 电功率，衡量"一度电能换多少可见光"，是节能的核心指标。

> **关于"绝对定标"——这点很容易踩坑**：上式能直接算出"绝对流明"的前提是 `S(λ)` 本身是**绝对辐射量**（如 W·nm⁻¹·sr⁻¹ 或经校准的 W·nm⁻¹）。但绝大多数分光光度计导出的 SPD 是**相对值**（常归一化到峰值=1 或某参考），此时 `∫S·ȳ` 只是相对数，算出的 Φ 是**相对光通量**，不能当"多少流明"读。要得到真正的绝对流明：
> 1. 用**已校准**的光谱辐射计（溯源到标准灯 / 积分球系统）测量，使 S(λ) 落在绝对单位上；或
> 2. 用一支**已知光通量的参考源**标定比例系数 `k`，再算 `Φ_abs = K_m · k · ∫S_rel·ȳ dλ`。
> 本文工具按输入 SPD 的单位直接套用该公式——喂相对谱就得相对值；要出绝对流明，请先把输入数据做好绝对定标。

—

## 八、把上面整套流程写成可运行 Python

        下面这段自包含代码，输入任意 SPD，就能算出色坐标、CCT、主波长、色纯度、峰值波长、半峰宽。配色函数用了 Wyman 2013 的解析近似，不用额外下载数据文件，**复制即可跑**：

import numpy as np

WL = np.arange(380, 781, 1, dtype=float)          # 380-780 nm, 1nm 步进

# —- CIE 1931 2度观察者配色函数（Wyman 2013 解析近似）—-
def _g(x, mu, s1, s2):
    s = np.where(x < mu, s1, s2)
    t = (x – mu) / s
    return np.exp(-0.5 * t * t)

def xbar(x): return 1.056*_g(x,599.8,37.9,31.0) + 0.362*_g(x,442,16,26.7) – 0.065*_g(x,501.1,20.4,26.2)
def ybar(x): return 0.821*_g(x,568.8,46.9,40.5) + 0.286*_g(x,530.9,16.3,31.1)
def zbar(x): return 1.217*_g(x,437,11.8,36) + 0.681*_g(x,459,26,13.8)

XB, YB, ZB = xbar(WL), ybar(WL), zbar(WL)

def spd_to_xyz(spd):
    X = np.trapezoid(spd*XB, WL); Y = np.trapezoid(spd*YB, WL); Z = np.trapezoid(spd*ZB, WL)
    s = X+Y+Z
    return X/s, Y/s, Z/s

# —- 示例：白光 LED（蓝光芯片 + 荧光粉宽带）—-
spd = (np.exp(-((WL-450)**2)/(2*18**2))
       + 0.9*np.exp(-((WL-565)**2)/(2*80**2))
       + 0.35*np.exp(-((WL-610)**2)/(2*70**2)))
spd = spd / spd.max()

x, y, _ = spd_to_xyz(spd)

# —- 精确 Planck 积分 + 牛顿法求 CCT（与工具内 calc_cct 实现一致）—-
def planckian_xy(T):
    c2 = 0.0143877696                     # 第二辐射常数 h·c/k (m·K)
    wl_m = WL * 1e-9                      # 波长 nm -> m
    S = 1.0 / (wl_m**5 * (np.exp(c2/(wl_m*T)) – 1.0))   # Planck 谱（相对值）
    X = np.trapezoid(S*XB, WL); Y = np.trapezoid(S*YB, WL); Z = np.trapezoid(S*ZB, WL)
    s = X+Y+Z
    return X/s, Y/s

# CIE 1960 UCS 坐标：普朗克轨迹在此空间近似直线，CCT 在此空间求解
def _uv(x, y):
    d = -2*x + 12*y + 3
    return 4*x/d, 6*y/d

# McCamy 给初值，牛顿法在精确 Planck 轨迹上精修
def calc_cct(x, y, max_iter=15, eps=1e-3, tol=1e-4):
    ut, vt = _uv(x, y)
    n = (x-0.3320)/(y-0.1858)
    T = float(np.clip(-449*n**3 + 3525*n**2 + 6823.3*n + 5520.33, 1000.0, 20000.0))
    for _ in range(max_iter):
        xx, yy = planckian_xy(T)
        u, v = _uv(xx, yy)
        du, dv = u-ut, v-vt
        f = du*du + dv*dv
        if f < 1e-12:
            break
        # 数值一阶/二阶导数（中心差分）
        xh, yh = planckian_xy(T+eps); uh, vh = _uv(xh, yh)
        xl, yl = planckian_xy(T-eps); ul, vl = _uv(xl, yl)
        up, vp   = (uh-ul)/(2*eps), (vh-vl)/(2*eps)
        upp, vpp = (uh-2*u+ul)/eps**2, (vh-2*v+vl)/eps**2
        fp  = 2*(du*up + dv*vp)
        fpp = 2*(up*up + vp*vp + du*upp + dv*vpp)
        if abs(fpp) < 1e-15:
            break
        T_new = float(np.clip(T – fp/fpp, 1000.0, 20000.0))
        if abs(T_new – T) < tol:
            T = T_new; break
        T = T_new
    return T

CCT = calc_cct(x, y)

# 峰值波长 & 半峰宽（白光 LED 通常指蓝光芯片峰的半峰宽）
peak = WL[np.argmax(spd)]
blue = np.exp(-((WL-450)**2)/(2*18**2))      # 蓝光芯片峰（高斯）
hb = blue.max()/2
ib = np.where(blue >= hb)[0]
fwhm = WL[ib[-1]] – WL[ib[0]]

# 主波长 & 色纯度（以等能白点 E 为参考）
We = np.array([1/3, 1/3])
# 光谱轨迹 + 紫线，构造边界多边形
lx = []; ly = []
for w in WL:
    m = (WL==w).astype(float)
    X = np.trapezoid(m*XB, WL); Y = np.trapezoid(m*YB, WL); Z = np.trapezoid(m*ZB, WL)
    sx, sy = X/(X+Y+Z), Y/(X+Y+Z); lx.append(sx); ly.append(sy)
lx, ly = np.array(lx), np.array(ly)
boundary = np.vstack([np.column_stack([lx, ly]), [lx[-1], ly[-1]], [lx[0], ly[0]]])  # 闭合紫线

P = np.array([x, y])
best = None
for i in range(len(boundary)-1):
    A, B = boundary[i], boundary[i+1]
    # 射线 We->P 与线段 A->B 求交
    v1 = P – We; v2 = B – A
    denom = np.cross(v1, v2)
    if abs(denom) < 1e-12: continue
    t = np.cross(A-We, v2)/denom
    u = np.cross(A-We, v1)/denom
    if t >= 0 and 0 <= u <= 1:
        inter = We + t*v1
        if best is None or t < best[0]:
            best = (t, inter)
dom_wl = None
if best:
    ip = best[1]
    # 交点对应的波长：找最近的轨迹点（紫线部分人为指定）
    d = np.hypot(lx-ip[0], ly-ip[1])
    k = int(np.argmin(d))
    dom_wl = WL[k]
    purity = np.linalg.norm(P-We)/np.linalg.norm(ip-We)

print(f"色坐标 x={x:.4f}, y={y:.4f}")
print(f"相关色温 CCT ≈ {CCT:.0f} K")
print(f"峰值波长 ≈ {peak:.0f} nm, 蓝光芯片半峰宽 FWHM ≈ {fwhm:.0f} nm")
if dom_wl: print(f"主波长 ≈ {dom_wl:.0f} nm, 色纯度 ≈ {purity*100:.1f}%")

        跑这段，你会得到类似：`x=0.3424, y=0.3380，CCT≈5057K，峰值≈452nm，FWHM≈42nm，主波长≈587nm，色纯度≈4.1%`——一个典型的中性白光 LED。

> 显色指数 CRI（R1–R15）的算法涉及 15 块标准色板和 CIE 13.3 配色函数，篇幅所限完整实现不展开，工程实现可参考 CIE 13.3-1995 官方文档。

—

## 九、懒人方案：我做的「LED 光谱性能分析工具」

        把上面这一整套（导入→计算→出图→导出）做成了一个**桌面小工具**，免装 Python、点几下就出结果：

– **导入**：支持 Excel（.xls/.xlsx）和 CSV，CSV 表头 `波长(nm),辐照度S(λ)`，波长范围 380–780 nm；
– **内置预设谱**：CIE D65（6504K）、CIE D50（5003K）、标准光源 A、以及 LED-R(625nm)/G(525nm)/B(455nm) 单色模型，可一键加载对照；
– **一键计算**：色坐标、CCT、主波长、色纯度、峰值波长、半峰宽 FWHM、光通量、光视效能、**显色指数 CRI（Ra + R1–R15 明细）** 全列出来；
– **双图联动**：左图 SPD 光谱曲线，右图 CIE 1931 色度图（含光谱轨迹、普朗克轨迹、当前点），算完自动刷新；
– **导出**：结果可导出 CSV，并支持设置导出波长间隔。

> 图为界面示意，实际功能以软件为准。

—

## 十、这工具能帮你省掉哪些重复活

        原理和代码上面都给了，代码部分完全可自用、不收费。真正耗人的往往不是"算一次"，而是"天天算、批量算、还要出报告"：

– **批量光谱计算**：产线 / 来料检验每天几十上百份分光数据，手动一份份算色温、显色指数不现实。工具支持一次导入多份 SPD，批量出全部指标，几分钟扫完一天的量。
– **自动导出检测报告**：算完直接导出 CSV，色坐标 / CCT / CRI（Ra + R1–R15）/ 光通量全字段带出，还能设导出波长间隔，省掉手抄和排版。
– **对接产线判色**：把"合格区间"做成对照，快速筛出偏色、色温漂移、显色指数不达标批次，比人眼稳得多。

        如果你正好卡在"批量 + 出报告 + 产线落地"这一步，想省掉搭环境、对付 Excel 的功夫，欢迎私信聊聊，或看看我的主页简介 😊。我也在持续更新「自研光电 & 办公小工具实战」系列，感兴趣可以关注一下。