做热辐射仿真和红外系统设计的朋友大概率都绕不开一条曲线普朗克曲线。它的用途很直接——描述黑体在不同温度下辐射能量随波长的分布从太阳光谱模拟到红外测温仪标定都要靠它打底。而曲线上的峰值点更是牵一发动全身峰值波长决定了辐射能量的“主战场”在哪温度一升高它就往短波方向跑这便是维恩位移定律的直观体现。Matlab里拿到这条曲线很容易plot一行就能出图但真正画得严谨、峰值找得准、单位不翻车还是有不少细节要抠。这篇内容就围绕“用Matlab绘制普朗克曲线并寻找峰值点”这件事把公式选择、单位统一、坐标策略、数值峰值检测和常见坑一次讲透适合刚接触热辐射仿真的学生、正在做红外传感器方案验证的工程师以及任何想把黑体辐射曲线用明白的人。1. 先把物理基础理清楚普朗克公式和峰值点意味着什么1.1 普朗克公式的两种写法别在单位上翻车普朗克公式最常见的形态是波长域版本B(λ, T) (2hc²) / λ⁵ × 1 / [exp(hc/(λkT)) - 1]其中 B(λ,T) 表示黑体在温度 T 下、波长 λ 处单位波长间隔内的辐射亮度h 是普朗克常量c 是真空光速k 是玻尔兹曼常量。这个公式描述的是温度越高的物体辐射能量越大而且峰值向短波方向移动。但实际使用时要注意公式还有一个频率域版本B(ν, T) (2hν³) / c² × 1 / [exp(hν/(kT)) - 1]这两个版本物理内涵相同但纵轴物理量不同一个是“每单位波长间隔的辐射亮度”一个是“每单位频率间隔的辐射亮度”。如果混用了转换关系画出来的曲线峰值位置会不同因为波长与频率不是线性比例单位间隔的宽度也不一样。不少新手把两个公式交叉代入结果峰位莫名其妙偏了一截回头还以为是Matlab的问题。更隐蔽的是波长λ的单位选米还是微米、纳米直接影响公式中的指数项量级。教科书里λ一般用米那么指数项里 hc/(λkT) 是一个很大的数容易出现数值溢出。1.2 峰值点背后的物理维恩位移定律怎么来的普朗克曲线峰值点的解析位置不必每次都用数值优化模拟维恩位移定律给出一个简洁的近似λ_max b / T其中 b ≈ 2.897772×10⁻³ m·K这个 b 不是经验拼凑而是从普朗克公式对λ求导并令导数为零推导出来的。推导过程中会出现一个超越方程它的数值根大约是4.965代入常量后就是2.898×10⁻³。所以严格来说维恩位移定律是从普朗克公式这个“老大哥”里长出来的结果两者不是并列关系而是推导关系。这条定律的价值在于只要知道温度就能立刻估算辐射能量集中区的中心波长。太阳表面约5800 K峰值在500纳米附近恰好落在可见光波段人体约310 K峰值在9.3微米附近落在红外热成像仪最关注的8到14微米窗口。这也是为什么红外体温检测设备要选这个波段——不是拍脑袋定的物理规律已经替你圈好范围了。1.3 找峰值点到底有什么用很多人以为找峰值点是为了在图上画个小箭头实际上它的工程价值远超“好看”。第一测温反向标定。测到黑体辐射的峰值波长后反推温度是单色测温的基础。虽然工程上很少只用峰值波长测温因为真实物体不是黑体但它是很多算法最初的出发点。第二波段选择的依据。做红外探测系统时探测器响应波段选在目标辐射峰值附近信噪比最高。想探测低温目标就该选长波红外想跟踪高温尾焰可以选中波红外甚至短波红外。峰值点告诉你该往哪个方向选。第三光源光谱估计。用普朗克曲线去拟合实际光源的测量数据拟合出的温度参数可以近似表示光源的色温哪怕是不完全黑体的白炽灯、卤素灯也能用这个模型做一个可参考的近似。理解了这三点再回过头看“绘制普朗克曲线并寻找峰值点”就不只是一个画图任务而是一个把物理规律转化成工程判断的入口。2. Matlab绘图前的关键准备单位、常量和坐标设置2.1 常量与单位先统一再写公式在Matlab里写物理公式第一件事不是写代码而是定单位。我建议统一用国际单位制波长λ用米温度T用开尔文那么辐射亮度B(λ,T)的单位就是 W·m⁻²·sr⁻¹·m⁻¹。这个单位读起来很拗口但它保证了公式系数不需要额外调整。实用常量建议直接写数值不要每一步都从别的程序里复制物理量数值单位普朗克常量 h6.62607015×10⁻³⁴J·s真空光速 c2.99792458×10⁸m/s玻尔兹曼常量 k1.380649×10⁻²³J/K维恩位移常数 b2.897772×10⁻³m·K如果习惯用微米作为波长单位代入公式前要把λ从微米换算成米也就是乘以1e-6。也有一种更“偷懒”的写法直接把 hc/k 这个组合量先算好因为指数项里总是一起出现可以减少代码中的常量数量降低出错概率。2.2 波长范围动态选择别让区间毁掉曲线初学者常见的错误是把波长范围写死比如恒定为0.1到3微米。这个范围画太阳辐射5800 K没问题峰值在可见光附近但画人体的热辐射310 K就有问题了峰值在9.3微米你画的区间根本没覆盖到峰值曲线看起来单调递减维恩位移定律也无从验证。更好的策略是根据温度动态生成波长范围。一个稳妥的经验是以维恩峰值λ_max为中心向左覆盖0.1倍峰值波长向右覆盖30到50倍峰值波长。为什么是50倍因为普朗克曲线的长波尾巴衰减很慢短波端下降极快长波端相对平缓向右多留一些范围才能展现完整形态。举个例子T 310 K 时 λ_max ≈ 9.3微米那么波长范围可以从0.5微米一直到300微米甚至更远。T 5800 K 时 λ_max ≈ 0.5微米范围可以从0.05微米到20微米。如果同时画多个温度建议用 linspace 或 logspace 分别生成每条曲线的波长向量不要共用同一个向量。2.3 线性坐标与对数坐标峰值位置看起来会“变”普朗克曲线在短波段能量迅速上升、长波段缓慢下降跨越的量级可能超过4个数量级。用线性坐标绘图时峰值尖锐长波尾巴被压缩成一条贴近横轴的线用半对数坐标横轴取log时曲线形态变得更像“钟形”。这里有个非常容易误导视觉的现象把横轴换成对数后峰值波长在图上看起来移动了。其实峰值点还是那个物理峰值点只是坐标变换改变了横轴上不同位置的比例尺度。眼睛看到的“中心”更多是对称中心的感觉不是真正的峰值。如果你要给人汇报我建议线性坐标为主因为物理教科书默认用线性横轴表达峰值如果你想在一张图里同时展示高温和低温曲线的全貌可以考虑画两张图或者用 subplot 分开展示。硬要把多温度曲线塞进同一张图且用同一个线性坐标范围低温曲线多半完全淹没。2.4 向量化写法让Matlab跑得更痛快Matlab号称矩阵实验室处理这种逐项计算的公式最忌讳用 for 循环一个一个点算。向量化的写法既简洁又不给解释器添堵。核心思路是先创建一个等间隔的波长向量然后用公式一次性算出对应的辐射亮度向量。这里要注意点乘与点除也就是在操作符前面加英文句点。写成 2hc^2 ./ (lambda.^5) 而不是 2hc^2 / (lambda.^5)后者变成矩阵运算结果完全不对。指数项同理exp(hc ./ (lambda * k * T))。lambda 是向量hc/(k*T) 是标量用 ./ 做逐元素除法即可。向量化之后即使波长向量有10万个点Matlab也能很快算完完全不需要往循环里钻。3. 手把手实现从普朗克曲线到峰值点检测3.1 完整基础代码画一张像样的普朗克曲线下面这段代码我习惯称为“黑体辐射三件套”生成波长轴、计算辐射亮度、绘图标注。别看结构简单每一步都有讲究。% 物理常量统一用 SI 单位 h 6.62607015e-34; % 普朗克常量J·s c 2.99792458e8; % 光速m/s kB 1.380649e-23; % 玻尔兹曼常量J/K b_wien 2.897771955e-3; % 维恩位移常数m·K T 5800; % 目标温度开尔文 % 动态生成波长范围从0.1倍峰值波长到20倍峰值波长 lambda_max_wien b_wien / T; lambda linspace(0.1 * lambda_max_wien, 20 * lambda_max_wien, 10000); % 普朗克公式波长域版本 % B (2*h*c^2) / lambda^5 * 1 / (exp(h*c/(lambda*kB*T)) - 1) x h * c ./ (lambda * kB * T); % 指数项系数 B zeros(size(lambda)); % 预分配防止溢出导致NaN valid x 700; % 指数过大会溢出 B(valid) (2 * h * c^2 ./ (lambda(valid).^5)) ./ (exp(x(valid)) - 1); % 绘图 figure(Color, w); plot(lambda * 1e6, B, LineWidth, 2); % 横轴转微米 xlabel(波长 \lambda (μm)); ylabel(辐射亮度 B(\lambda,T) (W·m^{-2}·sr^{-1}·m^{-1})); title(sprintf(普朗克曲线 T %d K, T)); grid on;关于溢出处理先说明一下x hc/(λkT)当波长很小或者温度很高时x 可能超过几百exp(x) 会变成 Inf随后 Inf-1 仍是 Inf再除以一个大数可能得到 NaN。代码里用 valid 索引只对不溢出的点计算剩余点保持0。这样做不是数学上绝对严谨但画图、找峰完全够用。计算完成后可以顺手验证一下辐射亮度曲线的基本特征短波端趋近0峰值附近快速上升长波端缓慢下降。3.2 峰值点检测的三种方法我推荐第二种检测峰值点至少有三条路可以走各有取舍。方法一离散网格找最大值。直接在波长向量上用 max 函数找到 B 的最大值和对应索引然后通过 lambda(idx) 读出峰值波长。这个方法最简单但精度受限于波长步长。步长太大峰值可能落在两个网格点之间即使最大点找对了真正峰值波长可能偏差几十纳米。方法二离散网格粗定位加抛物线插值精修。先粗略定位峰值附近的三点再用抛物线插值公式求真正的极值位置。步长不大的情况下这个方法的精度能提高好几个数量级而且实现代码只要几行不需要调用额外工具箱。方法三连续优化。定义函数句柄用 fminbnd 对负辐射亮度求极小值。注意要传负号因为 fminbnd 默认极小化而我们要极大值。这个方法理论精度最高但需要函数连续且初值区间合理性能上也比前两种重。当我们需要批量处理几千个温度点时逐次调用优化器没有优势。我日常用得最多的是方法二理由很简单快、稳、足够精确。峰值附近普朗克曲线用抛物线近似已经非常接近真实曲线三点插值得到的结果和连续优化几乎无差别。% 粗定位 [peak_val, idx] max(B); i0 idx; % 抛物线插值精修 if i0 1 i0 length(B) y1 B(i0-1); y2 B(i0); y3 B(i01); denom y3 - 2*y2 y1; if abs(denom) 1e-30 delta -0.5 * (y3 - y1) / denom; else delta 0; end lambda_peak lambda(i0) delta * (lambda(2)-lambda(1)); else lambda_peak lambda(i0); end fprintf(数值峰值波长: %.6e m\n, lambda_peak); fprintf(维恩定律峰值: %.6e m\n, lambda_max_wien);3.3 峰值点标到图上图表才完整峰值点找到以后光在命令行打印一个数字没有说服力建议直接标到图上。常用的做法是画一条垂直线然后在旁边加文字标注。hold on; plot(lambda_peak * 1e6, peak_val, ro, MarkerSize, 8, LineWidth, 1.5); line([lambda_peak * 1e6, lambda_peak * 1e6], [0, peak_val], ... LineStyle, --, Color, [0.85 0.33 0.1], LineWidth, 1.2); text(lambda_peak * 1e6 * 1.05, peak_val * 0.9, ... sprintf(\\lambda_{max} %.4f μm, lambda_peak * 1e6), ... FontSize, 11); hold off;标注位置需要微调。如果峰值点过于靠左文字放右边可能超出绘图范围可以把坐标乘以0.8放到左边或者直接使用 text 配合水平和垂直对齐属性。这一步不是难点但直接影响图表可读性。我习惯把峰值点和维恩定律理论值画在同一张图上用两种不同颜色区分方便一眼看出两者是否贴合。3.4 和维恩位移定律对照验证误差从哪来画完图、找到峰强烈建议做一步验证把数值检测结果和维恩位移定律的理论值对比一下。正常情况下两者在小步长下应该几乎一致相对误差不超过0.1%。如果误差明显偏大多半是三个原因波长向量范围没有覆盖峰值、步长太大、或者公式单位写错。举个例子T 5800 K维恩理论值是 2.897771955×10⁻³ / 5800 ≈ 4.996×10⁻⁷ m也就是约499.6纳米。如果波长向量是从0.1微米到20微米、网格点数10000那么间距约为1.99纳米抛物线插值后应该能得到499.6纳米左右的结果误差不到0.1纳米。如果只用100个网格点间距200纳米粗定位可能直接跳到500纳米附近但插值精度也会明显下降。所以网格点数太少时峰值点看起来也是对的其实已经偏离了真值。这个对照步骤相当于给整个计算过程装了一个“校准仪”一旦结果差得远说明前面某个环节出了问题越早发现越省时间。4. 常见问题与排查技巧实录4.1 exp溢出与NaN普朗克曲线最容易踩的坑普朗克公式里的指数项最需要提防。当x hc/(λkT)非常大时exp(x)会超出双精度浮点数范围。Matlab中exp(1000)是Inf因为双精度能表示的最大数约1.8×10³⁰⁸而exp(710)就已经接近这个上限了。溢出带来的连锁反应很奇怪有时不是Inf而是NaN。原因是若干个Inf参与后续乘除或者Inf/Inf直接产生NaN。NaN会在向量里“污染”数据导致max、插值、画图全部异常。更麻烦的是某些情况下plot不会报错但图上会出现一条向下扎到负无穷的尖刺看着像数据出错了其实道理很简单——指数项爆了。我的排查经验是三步走。第一步先对公式里的x做检测看看波长向量是否逼近0。第二步用x700作为有效区域掩码。第三步初始化B时用 zeros保证无效区的点自然为0而不是留下NaN。这里阈值700可以再放宽或收紧本质只要保证exp(x)不溢出即可。温度极高或者波长极短时这个方法依然有效。4.2 峰值点找不到或者明显偏了先怀疑坐标范围如果你发现max函数返回的最大值点恰好落在波长向量的端点那么峰值一定超出了你预设的范围。这甚至不算“找不到峰值”而是区间选择失误。举个例子你研究的是3开尔文的宇宙微波背景辐射峰值波长大约1毫米如果波长范围还停留在微米量级曲线在可见光波段基本为0max返回的是向量末端的最大值毫无物理意义。解决方法是动态范围。推荐在绘制之前先用维恩定律估一个λ_max再把范围设为 [λ_max/10, λ_max×30] 或更宽一些。注意长波端要留足余量因为普朗克曲线在长波侧衰减慢切得太早会让曲线看起来像被刀砍了一刀。4.3 曲线短波端异常抖动或出现突刺普朗克曲线在短波端应当单调上升并趋近0但有时候画出来出现锯齿或尖峰尤其是靠近波长下限的位置。最常见的原因是网格点数不足导致短波端急剧上升的部分被采样成折线。另一个可能原因就是你让波长向量里出现了0或负值x变成了Inf或负Inf最终算出的B不是有效数。锯齿问题好解决把波长下限抬高一点避开极端陡区或者增加网格点数到1万以上。我通常用10000到50000个点画出来非常顺滑。如果你用的是对数间隔的点长波端会稀疏需要特别注意峰值附近的采样密度是否足够。4.4 多温度曲线在同一张图里低温和高温的幅值差太悬殊怎么办普朗克曲线的总辐射能量随温度的四次方增长所以6000 K和300 K的曲线峰值幅度相差差不多16万倍。放在同一张线性图里低温曲线就是一条贴着横轴的“平地”什么都看不见。解决办法有几个。最直接的是分图展示用subplot分别画纵轴各自自适应。也可以用归一化处理每条曲线都除以自身的峰值这样形状保留能直观对比不同温度下峰值波长偏移的趋势。归一化后的图不会展示绝对能量差异但对物理趋势的表达非常清楚。如果非要放同一张图且保留绝对量级可以用对数纵轴也就是semilogy。不过要提醒自己纵轴取对数后峰值附近曲线看起来比较平缓不适合用来“目测”峰值位置需要依赖代码计算的数值结果。5. 进阶扩展把画曲线找峰值做成顺手的小工具5.1 封装一个普朗克函数省得每次重写如果你只是临时算一次脚本写法没问题。但如果隔三差五要画不同温度的曲线建议把公式封装成函数。我常用的版本长这样function B planck_lambda(lambda_m, T) % lambda_m: 波长单位米可为向量 % T: 温度单位开尔文 % 输出 B: 黑体光谱辐射亮度 W/(m^2·sr·m) h 6.62607015e-34; c 2.99792458e8; kB 1.380649e-23; x h * c ./ (lambda_m * kB * T); B zeros(size(lambda_m)); valid (x 700) (lambda_m 0); B(valid) (2*h*c^2 ./ (lambda_m(valid).^5)) ./ (exp(x(valid)) - 1); end函数封装完之后画一条曲线只需要三行代码。而且因为输出始终是干净的有效数组之后无论做max还是积分都不会被NaN拖后腿。我在处理多个物理场景时基本靠这一个函数打天下改动率很低。5.2 批量绘制不同温度自动输出峰值表一次分析多个温度时可以写一个循环把每个温度的峰值波长算出来汇总成表格。这是我做红外波段方案对比时最常用的方式。Tlist [300, 500, 1000, 2000, 3000, 5800]; results zeros(length(Tlist), 2); for ii 1:length(Tlist) T Tlist(ii); lambda_max_wien b_wien / T; lambda linspace(0.1*lambda_max_wien, 30*lambda_max_wien, 20000); B planck_lambda(lambda, T); [~, idx] max(B); % 这里可以再补抛物线插值不赘述 results(ii,1) T; results(ii,2) lambda(idx); end disp(温度(K) 峰值波长(m)); disp(results);输出表可以直接导入电子表格工具或者转成图表横轴温度、纵轴峰值波长画出来就是一条反比例曲线和维恩定律的曲线重合。这对理解辐射规律很有帮助。5.3 进一步用积分验证斯蒂芬-玻尔兹曼定律有了普朗克函数后工具箱还能继续扩张。对波长域辐射亮度在0到无穷区间积分理论上应该得到σT⁴/π的结果其中σ是斯蒂芬-玻尔兹曼常数约5.670374419×10⁻⁸ W·m⁻²·K⁻⁴。这个关系是利用黑体辐射计算总辐射出射度的基础。Matlab里用integral做数值积分sigma 5.670374419e-8; T 1000; lambda logspace(-8, -3, 10000); B_all planck_lambda(lambda, T); E_total trapz(lambda, B_all); % 或者 integral((l) planck_lambda(l, T), 1e-9, 1e-2) predicted sigma * T^4 / pi;算出来两个值应当非常接近。这个验证看起来多此一举但如果你把普朗克公式封装成了函数顺手检查一次就能确认公式系数没有写错。我曾经用这个验证抓出过一个麻烦当时把波长单位误写成纳米导致积分结果差了10⁹倍一查就明白了。最后再说一个我自己处理这类问题时的体会画普朗克曲线和找峰值点真正的门槛往往不在Matlab语法而是物理概念和数值稳定性的拉扯。学会用维恩位移定律预先判断结果再用数值方法去验证整个过程会顺畅很多。希望这篇东西能让你少走几步弯路把时间花在更有价值的辐射仿真与数据解读上。