简介SBDART 是科学界广泛使用的大气辐射传输计算程序该压缩包提供了面向 MATLAB 环境的调用脚本与运行示例面向进行光辐射模拟、遥感和气候研究的学生及科研人员。压缩包共包含 7 个文件其中 6 个为 .m 格式的 MATLAB 源文件主要覆盖模型主控函数、大气与地表参数设置、经纬度辅助计算以及不同情景下的运行样例另附 1 张 PNG 输出示意图便于对比模拟效果。整个压缩包仅 213KB体量轻便适合在个人电脑上快速部署和反复调试。目前已有 386 人浏览学习说明这套脚本对理解 SBDART 的调用流程具有较好的参考价值。通过研究源码中的参数设置和数据传递方式使用者能够学会如何配置大气廓线、地表反照率、太阳几何参数等输入条件掌握辐射通量与辐亮度结果的分析方法并在此基础上开展地表能量平衡、遥感信号模拟以及模型不确定性评估等应用。1. 为什么气象遥感工作者最终都会绕回 SBDART做辐射传输计算的人无论你最初用的是 MODTRAN、6S 还是 libRadtran只要遇到“想快速验证一个大气校正思路、手边又没有商业许可”的场景大概率会在某个深夜搜到 SBDART 这个名字。SBDART 是一套由 Scripps 海洋研究所开发的平面平行大气辐射传输模型计算短波和长波波段的光谱辐射量、辐照度与透过率内部采用离散纵标法求解辐射传输方程。它的最大特点是免费、输入输出全文本、单机单核就能跑通一条完整的大气路径。配合 MATLAB 做封装你能把整条“构建大气场景 → 计算辐射亮度 → 对比实测数据”的链路控制在自己的代码框架内不必依赖商业软件的黑匣子。这篇文章面向两类人一类是要用辐射传输模型做遥感反演或大气校正的研究生另一类是需要在项目里批量计算不同大气条件下传感器入瞳辐亮度的工程师。我会从模型原理讲到 MATLAB 调用细节最后给出几组实际运行中容易翻车的参数坑。2. SBDART 在算什么模型内核与 MATLAB 封装路线2.1 平面平行大气里的辐射传输方程SBDART 的核心是把大气分成若干层每层视为平面平行介质不考虑地球曲率和水平非均匀性。它计算的是在给定太阳天顶角、观测天顶角、地表反射率和大气廓线条件下辐射强度随波长和方向的变化。方程本身是经典的辐射传输方程dI(τ,μ)/dτ I(τ,μ) − ω/2 ∫ P(μ,μ′)I(τ,μ′)dμ′ − (ω/4π)F₀P(μ,−μ₀)exp(−τ/μ₀)其中 τ 是光学厚度μ 是观测天顶角的余弦ω 是单次散射反照率P 是散射相函数F₀ 是太阳辐射通量。SBDART 对散射相函数做勒让德多项式展开展开项数可以在输入文件里控制。理解这个方程对调试参数非常重要云层或气溶胶层设置不当本质上是把 ω 和 P 放错了位置结果出来的辐亮度会偏差几十个百分点。我的习惯是先跑一个“无云、无气溶胶、瑞利散射主导”的小场景确认基线和实测或参照模型一致再加入复杂条件。很多人在 MATLAB 里一上来就塞入气溶胶模式结果光谱曲线抖动得厉害根本分不清是模型问题还是输入格式问题。先把最简场景跑通后面才谈得上调参。2.2 MATLAB 调用的三种路线及其取舍SBDART 官方发行的是 Fortran 源码你拿到手通常需要自己编译。MATLAB 调用它有三种常见方案方案做法优点缺点命令行调用MATLAB 用 system() 调用编译好的 SBDART 可执行文件最稳和 Fortran 解耦每次调用都有 I/O 开销MEX 封装把 SBDART 的核心子程序编译成 mex 文件速度快省去磁盘读写编译环境复杂Fortran 与 C 接口易出错纯 MATLAB 重写把离散纵标法用 MATLAB 重写可单步调试不现实性能差距太大我一般推荐第一种。原因很简单SBDART 的本体是 Fortran它的输入文件是按固定列宽设计的文本格式你几乎不可能绕开文件读写。既然绕不开不如直接用 system() 调用让 MATLAB 负责生成输入文件和解析输出结果把 Fortran 当作一个计算内核。这样后续换版本、换机器只要重新编译一次 Fortran 就行了。2.3 编译 SBDART 时要注意的 Fortran 版本差异你从网上下到的 SBDART 源码不同压缩包里的 Fortran 代码风格差异很大。老版本基于 Fortran 77直接用 gfortran 编译即可。新版有些包含了较多动态数组用法需要 Fortran 90/95 编译器。常见做法是gfortran -O2 -o sbdart sbdart.f如果你在 Windows 上建议安装 MinGW-w64 的 gfortran然后用 MATLAB 自带的 MinGW 工具链来编——这样避免了系统里多套编译器互相干扰。也可以选择 Windows 下用 Cygwin 编译出 exe再拷贝到 MATLAB 工作目录。我自己踩过的坑是用旧版 g77 编译出的 exe 在 64 位 Windows 上运行到一半就崩换 gfortran 后一次通过。提示编译成功后用./sbdart input跑一次自带的测试用例比对输出文件与官方参考结果一致再开始 MATLAB 封装。3. 从零搭一个 MATLAB 调用 SBDART 的最小工作流3.1 准备目录与输入输出文件框架这个流程的目标是在 MATLAB 里生成一份 SBDART 输入文件调用可执行文件再把输出读回来。先建立标准目录结构我的习惯是project/ bin/sbdart // 编译好的可执行文件 inputs/ // 每次运行生成的输入文件 outputs/ // SBDART 输出的原始文本 results/ // MATLAB 解析后的 .mat 文件 run_sbdart_batch.m // 主控脚本这样的好处是批量跑的时候不会把中间文件混在代码里找翻车现场时只用看 inputs 和 outputs 两个目录。3.2 用 MATLAB 生成 SBDART 输入文件SBDART 的输入文件是典型的按列解析的固定格式不同版本支持的变量名略有不同但核心变量基本稳定。这里给出一个覆盖日常使用的最小模板function write_sbdart_input(fname, cfg) % 生成 SBDART 输入文件 % cfg 是 struct包含以下字段 % idim: 1辐射亮度, 2辐照度 % nstr: 流数一般取 4 或 8 % wlinf, wlsup, wlinc: 波长范围与增量单位微米 % sza: 太阳天顶角单位度 % phi: 相对方位角单位度 % nrefl: 地表反射率类型0朗博特 % albedo: 地表反照率0~1 fid fopen(fname, w); fprintf(fid, $INPUT\n); fprintf(fid, idim %d\n, cfg.idim); fprintf(fid, nstr %d\n, cfg.nstr); fprintf(fid, wlinf %f\n, cfg.wlinf); fprintf(fid, wlsup %f\n, cfg.wlsup); fprintf(fid, wlinc %f\n, cfg.wlinc); fprintf(fid, sza %f\n, cfg.sza); fprintf(fid, phi %f\n, cfg.phi); fprintf(fid, nrefl %d\n, cfg.nrefl); fprintf(fid, albedo %f\n, cfg.albedo); fprintf(fid, $END\n); fclose(fid); end这段代码的逻辑是把 MATLAB 端的 struct 数组转换成 SBDART 的文本输入。fprintf的格式控制必须严格对齐因为 SBDART 的 namelist 解析器对空格和等号有特定要求。$INPUT和$END是 namelist 的定界符不能漏。注意不同版本对变量名的支持不一样。老版本只认sza新版可能改成zenith。我先用grep查一下源码里的变量声明再写模板避免把时间花在格式报错上。设置idim1时输出的是某方向的辐射亮度idim2时输出的是上行和下行辐照度。两者的物理意义不同后处理方式也不同建议明确你的应用目标后再选。3.3 主控脚本调用、超时控制与输出解析生成输入文件后调用可执行文件并读回结果是封装的核心。一个健壮的主控脚本长这样function result run_sbdart(cfg, exe_path, input_dir, output_dir) % 单次运行 SBDART返回解析后的结构体 % exe_path: 可执行文件的完整路径 % input_dir/output_dir: 输入输出目录 if ~exist(exe_path, file) error(SBDART 可执行文件不存在: %s, exe_path); end in_file fullfile(input_dir, input.txt); out_file fullfile(output_dir, sbdart.out); write_sbdart_input(in_file, cfg); cmd sprintf(cd %s %s %s %s, ... output_dir, exe_path, in_file, out_file); [status, cmdout] system(cmd); if status ~ 0 error(SBDART 运行失败: %s, cmdout); end result parse_sbdart_output(out_file); end这段代码的关键点是system()调用时用cd切换工作目录确保 SBDART 的临时文件生成在 output_dir 下不会污染当前工作目录。status判断是必要的不要省否则遇到输入文件格式导致的运行时错误你根本不知道模型是否正常结束。另一个容易被忽视的坑是 Windows 下 MATLAB 的system()调用 exe 时路径中不能有中文和空格。如果你把工程放在D:\我的项目\sbdart\运行时会直接报“找不到文件”。统一用英文路径是一个很快就能养成的习惯。3.4 解析输出文件并转成 MATLAB 矩阵SBDART 的输出文件是纯文本表格不同模式下列名和行数都不同。解析脚本要灵活不能硬编码行号function data parse_sbdart_output(out_file) % 解析 SBDART 输出返回波长与辐射量两列 fid fopen(out_file, r); raw textscan(fid, %f %f, HeaderLines, 7); fclose(fid); wl raw{1}; val raw{2}; data.wavelength wl; data.value val; end上面这个解析器假设输出文件前 7 行是表头。但实际表头行数会随idim和输出波段数变化所以更稳的做法是先读前 20 行找到第一个全数字行再开始textscanfid fopen(out_file, r); C textscan(fid, %s, Delimiter, \n); lines C{1}; start_idx 1; for i 1:min(20, length(lines)) if ~isempty(regexp(lines{i}, ^\s*\d\.\d\s\d\.\d, once)) start_idx i; break; end end fclose(fid);这个方法的核心不依赖固定的表头行数而是通过正则表达式找到第一行“两个浮点数开头”的数据行然后用textscan从那一行开始解析。版本升级导致输出格式微调时你的代码大概率仍能正常工作。4. SBDART 关键参数精调大气模式、气溶胶与地表反射率4.1 大气模式的选择逻辑SBDART 内置了多种大气模式热带、中纬度夏季、中纬度冬季、亚极地夏季、亚极地冬季、美国标准大气等它们本质上是不同的温度、湿度、气压廓线。选错模式的影响是水汽吸收带的辐射传输结果偏差巨大特别是在 0.94 μm 和 1.13 μm 附近的水汽吸收谷。判断逻辑很简单你的研究区大致纬度带就是你的选择依据。高纬度冬季场景用亚极地冬季模式热带海洋场景用热带模式。如果你正在做的是多季节时间序列分析建议不选单一模式而是用 MATLAB 循环跑多个模式把不确定性范围量化出来atmos_modes {tropical, midsummer, midwinter}; for i 1:length(atmos_modes) cfg.atm atmos_modes{i}; result(i) run_sbdart(cfg, exe_path, input_dir, output_dir); end这种做法在论文或项目报告里非常有说服力你展示的结果不是“某个大气模式下的一条线”而是“多种大气模式下的一个不确定带”审稿人和甲方都会满意。4.2 气溶胶模式与光学厚度的组合陷阱气溶胶是 SBDART 里最容易出问题的地方。它提供几种默认气溶胶类型乡村型、城市型、海洋型、沙漠型等并且允许你指定 550 nm 处的光学厚度taurer。翻车点在于气溶胶模式不只是一个名字它背后关联的是复折射指数谱和粒径分布不同模式的散射相函数差异极大。气溶胶模式适用场景典型光学厚度乡村型内陆清洁地区0.1 ~ 0.3城市型工业污染区0.3 ~ 0.8海洋型清洁海洋上空气团0.05 ~ 0.2沙漠型沙尘输送期间0.2 ~ 0.6经验是先用边界层能见度反推光学厚度再选模式。能见度 10 km 对应 550 nm 光学厚度大约 0.5能见度 23 km 对应 0.2。如果手里只有 MODIS 气溶胶产品直接读取 550 nm AOD 填入taurer即可不要自己瞎猜量级。注意SBDART 中气溶胶光学厚度是全谱段按波长指数外推的angstrom参数控制这个指数。默认值 1.0 适合大多数场景但沙尘事件时要设到 0.2 左右否则短波段的辐射计算会偏掉。4.3 地表反射率朗博特假设的边界SBDART 默认地表是朗博特体用一个标量反照率描述。这在实际遥感场景中问题很大因为地表反射率随波长的变化非常剧烈。好在 SBDART 允许你通过alb数组逐波段指定反射率这是提高模型真实度的最重要一步。% 植被地表的示例反射率谱波长单位 μm wl_grid [0.4 0.6 0.8 1.0 1.2 1.6 2.0 2.5]; albedo [0.05 0.20 0.45 0.40 0.35 0.25 0.15 0.08];把这段插值后的反射率谱写入输入文件时要确认格式和 SBDART 要求的一致。翻车案例很常见写了新的反射率谱后结果没变化仔细查才发现输入文件里nrefl没改成 1表示使用光谱分辨率反射率模型还在用默认的朗博特值。我的经验是如果做的是海洋遥感直接设albedo0.05就够了如果做的是陆地上空必须按波段给反射率否则蓝波段和近红外波段的模拟辐亮度误差会大到不可用。这个参数是决定模型能不能落到实际反演的关键不要偷懒。5. 避坑指南SBDART MATLAB 运行中的常见故障与修复5.1 路径中存在中文或空格导致调用失败现象system()返回错误码 1cmdout显示“系统找不到指定的路径”但 exe 确实存在。原因Windows 下 MATLAB 的system()对中文路径支持不稳定尤其是经过反斜杠拼接的命令行中文字符被转成 ANSI 后和实际磁盘路径不一致。解决所有工作目录统一用英文、无空格命名。在 MATLAB 开头加一行assert(isempty(regexp(pwd, [^\x00-\x7F])), 工作目录不能含中文)提前拦下来。5.2 输出文件为空但 exe 正常结束现象status0但生成的输出文件只有表头没有数据行。原因输入文件里波长范围设置不合理。比如wlinf wlsup或者wlinc太大导致只有一个采样点SBDART 干脆跳过了计算。解决先打印输入文件的完整内容检查一遍。最常见的是微米和纳米单位搞混把 0.4~2.5 写成 400~2500模型内部解析时自动截断结果输出空表。5.3 输出的辐亮度数值量级异常大或异常小现象算出的入瞳辐亮度比 MODTRAN 结果大 10 倍以上或者全部为零。原因idim2辐照度模式输出的是 W/m²/μm而idim1辐亮度模式输出的是 W/m²/μm/sr。两者数值上差 π 倍或更多。很多人只用其中一种却在后处理时混用公式。解决在做数据对比前先确认输出文件和模型设置一致。在解析函数里做一个断言if cfg.idim 1 mean(data.value) 100 warning(辐亮度值偏大检查是否误用了辐照度输出); end5.4 MATLAB 版本更新后system()行为变化现象同一段代码之前跑得很好升级 MATLAB 后直接报错或者输出结果被缓冲。原因MATLAB 的system()在新版本底层 shell 改成了 PowerShell对输入输出重定向的处理策略不同。解决在调用前显式加上cat或者改用cmd /c前缀cmd sprintf(cmd /c cd %s %s %s %s, ... output_dir, exe_path, in_file, out_file);cmd /c向前兼容老版本和新版本是值得写入默认模板的细节。5.5 输入文件列宽对齐问题现象SBDART 读入参数时不报错但某些参数没有生效比如光学厚度还是默认值。原因namelist 解析按等号后第一个 token 取值模型初始化时可能仍在运行只是把参数跳过。这在 Fortran 的 list-directed I/O 里经常发生——字符串或数值之间缺了空格解析器把两个数并成一个数后续字段全部错位。解决生成输入文件时用固定格式输出比如sza一律用%8.2f不要用紧凑格式。用 MATLAB 生成后先打开文件人工目检一遍或者写一个自动检测脚本% 检查输入文件是否包含可疑的行 txt fileread(in_file); if contains(txt, NaN) || contains(txt, Inf) error(输入文件含非法数值); end6. 进阶用法多场景批量计算与结果可视化当你把单次调用跑通之后真正的生产需求是批量计算比如卫星影像每个像元的太阳天顶角都不同你要对每条观测几何分别运行 SBDART或者 MODIS 影像一幅有 100 万个像元虽然不用逐个跑通常会按角度分档但至少几十种角度组合是跑不掉的。批量计算的关键不是循环本身而是管理好中间文件。我用parfor并行跑多场景每个 worker 有独立的临时目录避免并发写入同一个输入文件scenarios build_scenarios(); % 生成所有角度和大气组合 parfor i 1:numel(scenarios) run_dir fullfile(tempdir, sprintf(sbdart_%04d, i)); mkdir(run_dir); result(i) run_sbdart(scenarios(i), exe_path, run_dir, run_dir); endparfor的注意点MATLAB 并行池里每个 worker 的system()工作目录是独立的但如果你把输出文件写到了同一个文件夹就会互相覆盖。为每次运行创建独立目录是必须的。跑完后统一清理临时目录for i 1:numel(scenarios) rmdir(fullfile(tempdir, sprintf(sbdart_%04d, i)), s); end结果可视化方面我习惯把不同场景的辐亮度谱叠在一张图上按角度或光学厚度渐变着色。这样可以立刻发现异常曲线——比如某条谱线在 0.76 μm 氧气吸收带位置出现奇怪的尖峰那多半是大气模式里的氧气浓度表出了问题赶上早检查早修复。最后说一个个人习惯正式跑批量前我会先构造两个极端场景太阳天顶角 0° 和 80°光学厚度 0.1 和 1.0用 MATLAB 的assert检查结果的单调性是否符合物理直觉。这一步只需要 5 秒但能省掉后面一整轮批量运行后才发现结果的尴尬。这个方案整体投入成本不高——一台普通电脑、一个 gfortran 编译环境、一份空闲的晚上就能把它落地。希望帮到你。本文还有配套的精品资源点击获取