用MATLAB计算普朗克公式：黑体辐射计算与单位换算全指南

简介：面向红外仿真与黑体辐射研究的MATLAB代码包，实现普朗克公式对辐射出射度的数值计算，适合需要分析不同温度与波长组合下黑体辐射特性的工程师、科研人员与相关专业学习者。普朗克公式是描述黑体辐射能量分布的经典定律，其数值计算在红外物理、遥感与热辐射分析中应用广泛。压缩包内共5个功能互补的.m文件，大小仅3KB，覆盖物理常数的定义、温度与波长范围的参数设置、双层循环遍历计算辐射度、基于slice函数的三维可视化，并附带曲线拟合与映射函数脚本，便于对辐射曲线作进一步处理和扩展。使用时只需修改温度与波长区间即可快速得到结果与三维辐射强度图，可直接支撑红外成像仿真、热辐射计算等场景，也可作为MATLAB数值仿真教学中的简明范例。目前已有4842人学习，对于需要掌握普朗克公式编程实现和MATLAB可视化表征的读者具有实用参考价值。 做热辐射计算的人，十有八九都写过“用MATLAB计算普朗克公式”这个小需求。无论你是做红外探测器标定、热成像系统仿真，还是做遥感大气校正、太阳能电池光谱响应分析，黑体辐射公式都是绕不开的起点。我最早接触这个需求是在做红外系统仿真的时候，为了估算不同温度目标在8~14μm波段的辐射亮度，需要先在MATLAB里把普朗克公式完整实现一遍。写完之后发现：公式本身并不难，真正让人反复折腾的是单位换算、积分区间的选择，以及如何用维恩位移定律、斯特藩-玻尔兹曼定律这类经典结论来验证代码正确性。这篇文章就围绕这几件事，把从零开始用MATLAB计算普朗克公式的完整过程，以及我在调试中踩过的坑整理出来，希望对正在做类似计算的朋友有帮助。

1. 一个公式，三种写法：先分清你算的是哪种普朗克公式

1.1 波长域、频率域、波数域，连换算规则都不同

做黑体辐射计算，第一步就容易被各种公式形式绕晕。同一套物理内容，波长域写一个形式，频率域写一个形式，波数域再换一套参数，网上随便一搜的公式经常长得完全不一样。实际上它们描述的是同一个物理对象：绝对黑体在热平衡状态下的光谱辐射亮度。区别只在于自变量选了波长还是频率。

我建议你直接锁定波长域形式，其他形式先别管。波长域普朗克公式的标准写法是：

B(λ,T) = (2hc²) / (λ⁵ × (e^{hc/(λkT)} − 1))

其中h是普朗克常数，k是玻尔兹曼常数，c是真空光速。这个形式的B是光谱辐射亮度，单位是W/(m²·sr·m)。很多初学者算出来的结果和参考值差好几个数量级，问题几乎都出在这个单位上：最后那个“/m”代表光谱密度，也就是说这个量是“每单位波长”的辐射亮度。如果你把光谱量当成总量来读，结果自然对不上。

如果换成频率域，公式长这样：

B(ν,T) = (2hν³) / (c² × (e^{hν/(kT)} − 1))

这里要特别提醒：波长域和频率域的曲线不能直接拿来对比。因为Bλ dλ = Bν dν才代表相同的辐射能量，所以两个密度函数之间必须乘一个坐标变换因子Bλ = Bν · c/λ²。换句话说，谱密度和自变量是绑定的，换了自变量，密度值必须跟着换算。这个细节在复现论文和外文资料时特别容易踩雷，我看到过不止一个人拿频率域的结果去对比波长域曲线，最后怎么都对不上。

1.2 工程常用“微米版”，两个常数怎么来的

实际工程计算里，更常用的是把常数合并好的“微米版”：

B(λ,T) = (1.191042×10⁸) / (λ⁵ × (e^{14387.77/(λT)} − 1))

这个形式里波长λ的单位是微米，温度T的单位是开尔文。那这两个看着很奇怪的常数到底怎么来的？其实很简单：把2hc²算出来是1.191×10⁻¹⁶ W·m²，再把波长从米换算成微米，也就是把λ⁵里的单位换掉，就得到1.191×10⁸；把hc/k算出来是1.4388×10⁻² m·K，换算成微米·K就是14387.77。这两个常数不是凭空蹦出来的，而是原始物理常数和单位换算系数合并后的结果。

我平时的习惯是：做理论验证时用SI单位版，保证和教材公式一一对应；写工程计算脚本时全部换成微米版，因为红外系统给出的谱段范围通常直接就是“8~14 μm”“3~5 μm”，没人愿意每次都把它们写成“8×10⁻⁶ ~ 14×10⁻⁶ m”再去参与运算。这里有个容易混淆的点：同一物理波长下，SI版输出的“每米”谱辐射亮度和微米版输出的“每微米”值之间差一个10⁻⁶系数。这不是程序bug，而是因为谱密度是“每单位波长间隔”的量，1μm = 10⁻⁶ m，坐标尺度变了，密度值自然等比例缩放。你可以把它理解为“每公里有几棵树”和“每米有几棵树”的差别——树的总量没变，单位密度数值完全不同。

2. 最小实现：把普朗克公式翻译成MATLAB函数

2.1 先写一个SI单位版本

写代码要遵循一个原则：对输入参数的单位做明确区分。否则三个月后回来看代码，大概率要猜这个变量到底是米还是微米。我一般这样封装：

function B = planck_lambda(lambda_m, T) % 普朗克公式，波长域，SI 单位输入 % lambda_m: 波长，单位 m % T: 温度，单位 K % B: 光谱辐射亮度，单位 W/(m^2 sr m) h = 6.62607015e-34; c = 2.99792458e8; k = 1.380649e-23; B = (2*h*c^2 ./ lambda_m.^5) ./ (exp(h*c ./ (k .* lambda_m .* T)) - 1); end

这个版本的好处是逻辑和物理定义一一对应，适合做正确性验证。但注意我写的是lambda_m.^5、lambda_m .* T，全部用点运算。MATLAB里^默认矩阵幂，如果你写成lambda_m^5而输入是行向量，会直接报维度错误。向量化是后面所有绘图和积分操作的基础，没有这一步什么都做不了。

2.2 工程计算用微米版

另一个我更常用的版本是微米版，用了合并后的常数：

function B = planck_lambda_um(lambda_um, T) % 普朗克公式，波长单位 μm，温度单位 K % 输出单位 W/(m^2 sr μm) c1L = 1.191042e8; c2L = 14387.77; B = (c1L ./ (lambda_um.^5)) ./ (exp(c2L ./ (lambda_um .* T)) - 1); end

封装完函数后，建议先做一个手算验证。以10μm、300K为例，把λ=10、T=300代入微米版：指数参数为14387.77/3000 ≈ 4.796，e的4.796次方约121，减1后约120；分子1.191042×10⁸除以10⁵是1190，除以120后约9.9 W/(m²·sr·μm)。你可以在命令行里跑一下，如果输出和这个量级明显不一致，说明函数或调用方式有问题。每次写完这类函数，我都建议先找一个已知点做这种“手算校核”，这是最高效的自检手段，比直接扔进大项目里出错后再回头查要省事得多。

3. 光画图不够，还得验证维恩位移定律

3.1 一次性画出多条温度曲线

有了函数，画图就是水到渠成的事。我的习惯是把不同温度的曲线画在同一张图里，直观比较峰值移动和曲线整体形状：

lambda = linspace(0.1e-6, 30e-6, 5000); T_list = [3000, 4000, 5000, 5800]; figure; hold on; for T = T_list B = planck_lambda(lambda, T); plot(lambda*1e6, B, 'LineWidth', 1.5); end hold off; xlabel('\lambda (\mum)'); ylabel('B_\lambda (W m^{-2} sr^{-1} m^{-1})'); legend(arrayfun(@(t) [num2str(t) ' K'], T_list, 'UniformOutput', false));

这里我把横轴换算成了微米，方便工程阅读。波长范围取0.1~30μm，能覆盖3000K到5800K的峰值区域。如果你要研究室温物体，比如300K，那曲线峰值在10μm附近，横轴取1~50μm更合适。画出来的曲线特征很明显：温度越高，整体辐射亮度越大，峰值波长越短。这是普朗克公式本身决定的规律，也是维恩位移定律的直观体现。

如果你想比较宽温域范围内的小信号和大信号，建议用semilogy画对数纵轴。线性图会让低温曲线几乎贴在零轴上，什么都看不出来；对数图能把200K到2000K每个温度下的曲线层次都拉开，这对于观察低温下的长波红外辐射特别有用。

3.2 从曲线中定位峰值，验证λ_max·T = b

维恩位移定律的常见写法是λ_max·T = 2898 μm·K。用数值方法找峰值很简单：

T = 5800; lambda = linspace(0.05e-6, 5e-6, 20000); B = planck_lambda(lambda, T); [~, idx] = max(B); lambda_peak = lambda(idx); fprintf('峰值波长: %.2f nm\n', lambda_peak*1e9); fprintf('λ_max*T = %.1f μm·K\n', lambda_peak*1e6*T);

运行之后你会发现λ_max·T大致在2898附近，但很少刚好等于2898。原因有两个：一是数值离散化，线性网格上最大值点不正好落在连续函数极值处；二是浮点精度。解决办法是“先粗扫、再加密”：第一次定位到大致位置后，在峰值附近重新用更密的网格搜索。比如第一次找到0.5μm附近，就重新linspace(0.45e-6, 0.55e-6, 50000)，精度可以轻松到0.1nm以内。

这个步骤千万别省。我实测过，如果只用100个点的粗网格，峰值波长甚至可能偏差到530nm而不是理论上的500nm。用来验证定律时这个偏差还能接受，但如果做探测器响应标定或者滤光片通带设计，几纳米的偏差就可能让设计完全跑偏。这也是数值计算和理论公式之间一个很有趣的差别：理论是连续的，数值计算永远离散，关键是你怎么把离散误差控制在可接受范围内。为了对数量级有感觉，可以对照下面这组理论峰值波长：

这组数据来自λ_max = 2898/T。当你的数值结果和这组数据明显不一致时，先别急着怀疑离散化，回去检查单位和函数实现，那才是大概率出问题的地方。

4. 波段积分：从曲线到工程上真正关心的数字

4.1 计算8~14μm波段辐射亮度

很多工程场景关心的不是一个波长点上的辐射亮度，而是某个谱段内的积分值。拿热红外系统举例，8~14μm是大气窗口，300K地面目标在这个窗口内的辐射亮度，直接决定了探测器的信号水平。画一条曲线只是定性认识，定量积分才是设计的输入。

实现核心就是一行积分，两种方式都可以：

T = 300; lambda = linspace(8e-6, 14e-6, 50000); B = planck_lambda(lambda, T); radiance_band = trapz(lambda, B); disp(radiance_band);

trapz的原理是把积分区间切成很多小段，每段近似为梯形累加。普朗克函数在8~14μm上没有奇点，形状平滑，5万点足够把误差压得很小。如果喜欢用自适应积分，可以写：

radiance_band = integral(@(lambda) planck_lambda(lambda, T), 8e-6, 14e-6);

integral会自动加密采样，对光滑函数可能更精准，但每次调用都有额外开销。我的建议：一次性的波段积分用trapz，因为网格和区间你都看得见摸得着，方便复核；如果要在循环里反复扫温度或扫波段，用integral更省心。

这里再强调一遍单位：对同一个8~14μm波段，如果用SI版函数，积分变量要从8e-6积到14e-6，结果单位是W/(m²·sr)；如果换成微米版函数，积分变量是从8积到14，虽然数值上同样都是对6个单位的波长宽度积分，但物理单位不同，结果量级也不同。把两个结果直接混着比较，是波段积分中最常见的错误来源。

4.2 用斯特藩-玻尔兹曼定律做全谱自检

拿到波段积分后，怎么确认这个数是对的？最可靠的办法是用斯特藩-玻尔兹曼定律验证全波谱积分：

∫₀∞ B_λ(λ,T) dλ = σT⁴ / π

其中σ=5.670374419×10⁻⁸ W/(m²·K⁴)。为什么除以π？因为辐射亮度是单位立体角上的量，半球出射度和亮度之间差一个π的几何系数。数值积分做不到真的从0积到∞，但可以把积分区间压缩到“远小于峰值、远大于峰值”的范围：

T = 2000; sigma = 5.670374419e-8; L_theory = sigma * T^4 / pi; L_numeric = integral(@(lambda) planck_lambda(lambda, T), 1e-9, 1e-3); rel_err = abs(L_numeric - L_theory) / L_theory; fprintf('理论值: %.6e W/(m^2 sr)\n', L_theory); fprintf('数值值: %.6e W/(m^2 sr)\n', L_numeric); fprintf('相对误差: %.2e\n', rel_err);

对T=2000K，峰值波长约1.45μm，积分下限取1e-9m比峰值短了三个数量级，上限取1e-3m比峰值长了近三个数量级，两端的剩余贡献已经小到可以忽略。实测下来相对误差通常在1e-6~1e-8量级。如果你算出来的误差明显偏大，那多半不是数值方法的问题，而是函数实现或者单位换算出了问题。这个自检方法也是普朗克公式计算从“能跑”到“可信”的关键一步。

5. 单位、溢出和向量化：我踩过的三个坑

5.1 第一坑：SI版本和微米版本混用

我在一个红外测温项目里吃过一次大亏。当时已经封装好了两个版本的普朗克函数，但因为某个脚本赶时间，一会儿用微米版，一会儿用SI版，而且没在变量名上区分单位。某次标定时发现300K黑体的积分辐射亮度比理论值高了十几个数量级，排查了一整天才发现：积分区间的上下界写成了微米值，但函数调用的是SI版，这个组合等于把一个微米的跨度当成米来积，量级自然全乱了。后来我给自己立了三条规矩：函数名带_si或_um后缀，输入变量名带单位后缀，每个函数头部写清楚输出单位。代码是写给下一次的自己看的，单位信息无论如何强调都不过分。

5.2 第二坑：exp溢出和expm1的使用

普朗克公式里有exp(hc/(λkT))，当λ很短，比如取到1e-9m，同时T又比较低，指数参数很容易超过709，MATLAB会返回Inf。虽然被积函数在这种条件下趋近于0，实际运算中却可能因为Inf参与加减得到NaN，然后顺着向量传染开来。最常见的触发场景是：为了做全谱积分，把波长向量取到1e-9m以下，然后函数输出一堆NaN。

解决方式有两个。一是用expm1(x)替代exp(x)-1，它对很小的x也可以避免灾难性抵消。二是对指数参数做一个截断，把它限制在700以内：

function B = planck_lambda_robust(lambda_m, T) h = 6.62607015e-34; c = 2.99792458e8; k = 1.380649e-23; x = h*c ./ (k .* lambda_m .* T); x = min(x, 700); B = (2*h*c^2 ./ lambda_m.^5) ./ expm1(x); end

700这个值不是随便拍的。double浮点数的最大值约1.8×10³⁰⁸，而exp(709)已经逼近这个边界，留一点余量就能避免在边界附近出现Inf。从物理角度看，当指数参数大于700时，指数项已经大到让整个辐射亮度趋于0，截掉这一段不会对结果产生任何可感知的影响。

5.3 第三坑：漏掉点运算符

这看起来是最基础的问题，但高频使用中真的会反复出现。从C或Java转过来的人尤其容易踩：写出2*h*c^2 / lambda_m^5，遇到数组输入就报“矩阵维度必须一致”。MATLAB里*和^默认是矩阵运算，对数组做逐元素运算必须用.*和.^。普朗克公式恰恰是个典型逐元素计算公式，每个波长点独立计算，完全不涉及矩阵乘法。这个坑不算深，但是每踩一次就会浪费十几分钟去盯错误信息。

以上三个坑，单位问题最难排查，因为它在运行时不报错，结果相对值也可能“看起来合理”，只有和理论值或参考数据对比时才会暴露；溢出问题在宽谱计算中很常见；点运算报错最直接，但次数多了会养成写代码时先检查点运算的习惯。我的经验是：写完函数先跑一个已知点的手算验证，再跑一次全谱积分自检，两步都过了才把函数放心地挪进正式项目。这套流程看起来多花几分钟，实际上能帮你省下后面几天的排查时间。

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