简介:本资源是面向大气科学、遥感与气象建模方向的科研学习者及高年级本科生/研究生的MATLAB版SBDART辐射传输模型实践包,用于快速开展光谱辐射传输模拟,解决太阳辐射在大气-地表系统中吸收、散射与反射过程的定量分析问题。压缩包共7个文件(6个MATLAB脚本+1张示例截图),总大小213KB;其中sbart.m为核心驱动程序,example1b.m至example3b.m为典型场景案例(如不同大气廓线、地表反照率与太阳天顶角组合),latlon.m与live_example.m支持地理坐标输入与交互式演示,screenshot1.png直观展示运行界面与输出效果。目前已有385人学习下载,资源结构简洁、即装即用,附带可直接运行的完整参数配置与结果可视化逻辑,显著降低SBDART模型入门门槛,助力用户掌握辐射传输数值模拟的关键流程与物理内涵。 我去年年底在一个遥感反演项目里被辐射传输计算卡了整整三周,原因很简单:手头只有 MODTRAN 的学生版授权,波段范围和分辨率都受限,想算一下大气顶辐亮度都得排队等 license。后来组里前辈甩给我一个 SBDART 的离线包,说“这个不用买授权,Matlab 里自己调就行”,我一开始将信将疑,结果用下来发现这玩意儿虽然不是商业软件那种开箱即用,但胜在源码开放、物理过程透明、改起来也顺手。今天就把我在 Matlab 环境下调用 SBDART 做辐射传输计算的经验完整梳理一遍,从物理基础到代码实现再到踩坑记录,希望能帮你少走弯路。
SBDART 全称是 Santa Barbara DISORT Atmospheric Radiative Transfer,它本质上是一套基于平面平行大气假设的辐射传输求解程序,核心求解器用的是 DISORT(Discrete Ordinate Radiative Transfer)算法。简单说,你给它大气剖面、气溶胶光学厚度、地表反射率这些输入,它就能算出你关心的辐亮度、辐照度、透过率等物理量。它非常适合做遥感器仿真、大气订正算法验证、辐照度估算这类工作,尤其适合学校里没有商业辐射传输软件授权的研究组。这篇文章适合遥感、大气物理、环境科学这些方向的研究生和工程师阅读,也适合刚接触辐射传输建模、想快速上手一套开源工具的人。
我需要说明一点:SBDART 官方原版是 Fortran 写的,网上流传的SBDART_matlab.rar这类压缩包里,通常有两种东西,一种是直接把 Fortran 源码编译成可执行文件,然后用 Matlab 去调用;另一种是有人做过 mex 封装,把核心计算函数包装成 Matlab 可以直接调用的接口。两种方式各有优劣,下面我会分别讲。
1. 方案选型:为什么在 Matlab 里调 SBDART,而不是直接用 Fortran 或 Python
1.1 三种调用方式的对比
我在决定用 Matlab 调 SBDART 之前,其实纠结过一阵子:是直接写 Fortran 主程序,还是转用 Python 的 PySBDART,或者干脆用 Matlab 调可执行文件。三者的差别还是挺明显的,我整理了一个对比表:
| 方案 | 优点 | 缺点 | 适合场景 |
|---|---|---|---|
| 原生 Fortran 编译调用 | 执行效率最高,无中间层开销 | 输入文件写起来麻烦,数据结构不直观,调试循环参数很痛苦 | 大批量长期跑,且不介意 Fortran 编程 |
| Python + PySBDART | 生态好,numpy/xarray 配合方便,社区活跃 | 需要另外配置 Python 环境,遥感老团队不一定有 | 新项目、团队以 Python 为主 |
| Matlab 调用 SBDART | 与遥感图像处理、矩阵运算无缝衔接;直接出图;学校常有正版授权 | 封装层可能带来一定性能损失;网上流传的包质量参差不齐 | 需要频繁交互调试、可视化、与影像处理流程耦合 |
对我个人来说,最大的痛点不是计算本身,而是输入参数的批量构造和结果的可视化。Matlab 的脚本环境刚好把这两点都解决了:我可以把 SBDART 的输入文件写成字符串模板,用循环批量替换参数,一次跑几十组大气条件;跑完直接imagesc出辐亮度分布图,再叠加个色标就是论文插图了。这套流程用 Fortran 写会非常折磨,用 Matlab 就舒服得多。
1.2 为什么 SBDART 适合做辐射传输仿真
相比 MODTRAN、6S 这类工具,SBDART 的定位一直是“科学研究的便捷工具”,它的物理过程覆盖得很全,包括瑞利散射、气溶胶吸收与散射、水汽和臭氧吸收、云层多次散射、地表双向反射等。它的几个特性特别戳中我:
- 平面平行大气假设下的辐射传输方程求解非常成熟,DISORT 算法的数值稳定性和精度都是经过几十年检验的;
- 内置了多种标准大气模型(热带、中纬度夏季、中纬度冬季、亚北极夏季、亚北极冬季、1976 美国标准大气),不需要自己去 ECMWF 拉探空数据;
- 气溶胶模型至少支持对流层、平流层、海洋型、沙漠型等多种预设,还可以自定义光学参数;
- 波长范围覆盖 0.25 到 100 微米,从紫外到远红外都行,遥感常用波段全包含。
如果你只是想在项目前期做一个“大气效应有多大”的敏感性分析,或者需要给传感器设计提供一组不同大气条件下的辐亮度曲线,SBDART 是性价比极高的选择。
2. 核心机制:SBDART 的输入文件、运行逻辑和物理参数
2.1 SBDART 是怎么工作的
SBDART 的运行方式并不是那种图形界面点一点就出结果,而是通过一个文本文件(通常是INPUT)传递参数,程序读入后调用 DISORT 求解辐射传输方程,最后输出各个波段的辐亮度、辐照度、透过率等物理量。
我画一个简单的流程理解:
构造输入文件(大气模型、观测几何、波长、气溶胶等) ↓ SBDART 主程序读取并解析输入 ↓ DISORT 求解平面平行大气辐射传输方程 ↓ 输出多个物理量文件(辐亮度、辐照度、透过率等) ↓ Matlab 读取输出文件,可视化或后处理这个链路里最核心的是输入文件里那几十个参数的物理含义。你如果不理解这些参数,就很容易出现“模型跑通了但结果完全是错的”这种尴尬情况。下面我把关键的参数逐个说明。
2.2 输入参数的核心物理含义
SBDART 的输入文件本质是一系列 Fortran namelist 风格的变量赋值。每个变量控制辐射传输计算的一个维度。我把常用的参数按功能分组整理如下:
大气与几何参数
| 参数名 | 含义 | 典型值/说明 |
|---|---|---|
nstr | 离散纵标法的流数 | 通常取 4 或 8,越大精度越高、越慢;遥感辐亮度计算建议至少 8 |
nlyr | 大气分层数 | 默认 60 层左右,可以改,但一般不用动 |
idatm | 大气模型编号 | 0=美国标准大气,1=热带,2=中纬度夏季,3=中纬度冬季,4=亚北极夏季,5=亚北极冬季 |
sza | 太阳天顶角(度) | 0 表示太阳在头顶,90 表示水平,一般 0~80 之间 |
phi | 相对方位角(度) | 太阳和观测方向之间的方位角差 |
umu | 观测天顶角的余弦值 | 比如天底观测 umu = 1,斜视 60 度时 umu = cos(60°) = 0.5 |
波长与光谱参数
| 参数名 | 含义 | 典型值/说明 |
|---|---|---|
wlmmin/wlmmax | 波长范围(微米) | 比如可见光 0.4 到 0.7 |
wlres | 波长分辨率(微米) | 0.005 对应约 5 nm,精度越高计算越慢 |
isat | 卫星/传感器波段模式 | 0 表示按波长扫描,1 表示用内置传感器响应函数 |
iout | 输出物理量类型 | 1=辐亮度,2=辐照度,3=透过率,4=路径辐射分量 |
气溶胶与云参数
| 参数名 | 含义 | 典型值/说明 |
|---|---|---|
iaer | 气溶胶模型 ID | 0=无气溶胶,1=对流层,2=平流层,3=海洋型,4=沙漠型等 |
taer55 | 550 nm 处气溶胶光学厚度 AOD | 0.1=清洁大气,0.5=中度污染,1.0=重污染 |
baer55 | 550 nm 处气溶胶单次散射反照率 | 默认 0.9 左右,吸收性强的烟尘会低于 0.8 |
caer55 | 550 nm 处气溶胶不对称因子 | 默认 0.7 左右,影响前向散射强度 |
wlin/wlout | 云层的水滴有效半径 | 典型 10 微米左右,用于云光学计算 |
tcloud | 云的光学厚度 | 晴空为 0,云越厚值越大 |
地表参数
| 参数名 | 含义 | 典型值/说明 |
|---|---|---|
alb | 地表反照率(单值,假设朗伯体) | 0.05=水体,0.15=植被,0.3=裸土,0.8=雪 |
uvflx | 是否计算紫外辐射通量 | 1=计算 0.28~0.4 μm 的紫外通量 |
lamb/lambc | 地表光谱反照率选项 | 可以用多个波长的反照率拟合地表反射谱 |
这些参数看着多,但实际常用组合就那么几套。后面我给的模板可以直接复制改参数,不用每次从零构造。
2.3 输出文件的物理含义
SBDART 的输出文件不止一个,常用的有这几个:
OUTPUT:主输出文件,包含辐亮度/辐照度随波长的变化,格式是两列(波长,物理量值);TAPE7、TAPE8: 中间计算文件,存了各层通量和 radiance 的中间结果;- 如果你设置了
iout=3,输出的是透过率而不是辐亮度,需要看清单位。
在 Matlab 里处理输出文件最常用的方式就是textread或load,把两列数据读完直接画图。但这里有个坑:SBDART 输出文件的行首可能带有空格或者 Fortran 格式的空格分隔,直接load有时会因为列数判断错误而失败,建议用textscan按格式读取。
3. 搭建在 Matlab 中调用 SBDART 的完整环境
3.1 环境准备与文件解压
网上流传的SBDART_matlab.rar解压后,一般会有这些内容:
SBDART_matlab/ ├── sbdart.exe # 编译好的 SBDART 可执行程序 ├── sb2mat.m # Matlab 封装函数(读取输出) ├── run_sbdart.m # 示例运行脚本 ├── INPUT # 示例输入文件 ├── OUTPUT # 示例输出文件 └── TAPE5, TAPE7, TAPE8 # 中间文件拿到压缩包后,我建议你先做三件事:
- 把整个文件夹放到一个纯英文路径下,不要有中文和空格,SBDART 的 Fortran 程序对路径里的特殊字符非常敏感;
- 先跑一遍
run_sbdart.m,确认能正常在当前环境运行,再动手改自己的参数; - 检查
sbdart.exe能否在命令行直接执行,如果不能,可能是缺少动态链接库或权限问题。
如果你拿到的是 Fortran 源码而不是编译好的 exe,那就需要自己编译。Windows 下可以用 gfortran(MinGW-w64 自带),Linux/macOS 下系统装过 gfortran 就能直接编。编译命令大致是:
gfortran -O3 -o sbdart sbdart.f disort.f *.f注意 SBDART 源码依赖 DISORT 的子程序文件,必须一起编译,很多初学者只编了主文件,结果链接报错一堆未定义引用。
3.2 Matlab 调用 exe 的基本方式
Matlab 调用外部可执行程序,用的是system或unix函数。基本模式是这样的:
% 切换到 SBDART 工作目录 cd('D:\projects\sbdart_matlab'); % 调用 SBDART 执行计算 [status, cmdout] = system('sbdart.exe < INPUT');上面这行的逻辑是:把INPUT文件内容作为标准输入喂给sbdart.exe,程序跑完会生成新的OUTPUT文件,status返回 0 表示运行成功,cmdout会捕获输出信息方便排查错误。
有人会问:为什么要用< INPUT重定向,而不是直接把参数作为命令行参数传进去?因为 SBDART 的设计就是读标准输入,它的执行流程是硬编码的,不接受命令行参数。所以“构造输入文件 - 调用 exe - 读取输出文件”是标准姿势。
3.3 用脚本批量生成输入文件的技巧
你要是跟我一样需要跑几十组不同 AOD、不同太阳天顶角的组合,手动改 INPUT 文件会改到崩溃。这里我分享一个我自己用的批量生产方法:用 Matlab 的字符串模板 +fprintf写文件。
你可以把 INPUT 文件想象成一个填空模板,里面只替换需要变的参数。我一般是这样写的:
function write_sbdart_input(filename, sza, aod, wlmin, wlmax) % 生成 SBDART 输入文件 fid = fopen(filename, 'w'); fprintf(fid, '$INPUT\n'); fprintf(fid, ' nstr = 8\n'); fprintf(fid, ' idatm = 2\n'); % 中纬度夏季大气 fprintf(fid, ' sza = %.2f\n', sza); % 太阳天顶角 fprintf(fid, ' phi = 0.0\n'); fprintf(fid, ' umu = 1.0\n'); % 天底观测 fprintf(fid, ' wlmmin = %.4f\n', wlmin); fprintf(fid, ' wlmmax = %.4f\n', wlmax); fprintf(fid, ' wlres = 0.005\n'); fprintf(fid, ' iaer = 1\n'); % 对流层气溶胶 fprintf(fid, ' taer55 = %.2f\n', aod); % 550nm AOD fprintf(fid, ' baer55 = 0.90\n'); fprintf(fid, ' caer55 = 0.70\n'); fprintf(fid, ' alb = 0.15\n'); % 植被地表 fprintf(fid, ' iout = 1\n'); % 输出辐亮度 fprintf(fid, '$end\n'); fprintf(fid, '\n'); fclose(fid); end这样写好函数之后,批量运行就变成了几行代码的事:
sza_list = [0 20 40 60]; aod_list = [0.1 0.3 0.5 1.0]; for i = 1:length(sza_list) for j = 1:length(aod_list) fname = sprintf('INPUT_sza%d_aod%d', sza_list(i), round(aod_list(j)*10)); write_sbdart_input(fname, sza_list(i), aod_list(j), 0.4, 0.7); % 运行 SBDART [status, ~] = system(sprintf('sbdart.exe < %s', fname)); if status ~= 0 warning('运行失败: %s', fname); end % 复制输出文件,避免被下一次运行覆盖 copyfile('OUTPUT', sprintf('OUTPUT_sza%d_aod%d.txt', sza_list(i), round(aod_list(j)*10))); end end这段代码是我在项目里实际用的简化版。核心思想就一条:每次运行前后都要备份输入和输出文件,因为 SBDART 的输出文件名是固定的,不备份下一次运行就会把上一次的结果覆盖掉。这个坑我刚开始踩过很多次,辛辛苦苦跑了一晚上,第二天发现 OUTPUT 文件只剩最后一组的结果。
4. 核心调用逻辑与 Matlab 封装函数解析
4.1 自己写一个sb_run函数,少走弯路
网上那个压缩包的sb2mat.m只能读取输出,不能帮你跑程序。我建议你在此基础上再封装一个更高层的函数,把“写输入 - 运行 - 读输出”三个步骤串起来,这样后面每次调用就只需要传物理参数,不用关心文件细节。
我写的函数大致长这样:
function result = sb_run(sza, aod, wlmin, wlmax, varargin) % SB_RUN 执行SBDART辐射传输计算,返回结构体结果 % 输入: % sza - 太阳天顶角(度) % aod - 550nm处气溶胶光学厚度 % wlmin - 起始波长(微米) % wlmax - 结束波长(微米) % 可选参数: % 'iaer' 气溶胶模型,默认1 % 'alb' 地表反照率,默认0.15 % 'umu' 观测天顶角cos值,默认1 % 'nstr' 离散纵标流数,默认8 % 返回: % result.wavelength - 波长向量 % result.radiance - 辐亮度向量 % result.status - 运行状态 % 解析可选参数 p = inputParser; addParameter(p, 'iaer', 1); addParameter(p, 'alb', 0.15); addParameter(p, 'umu', 1.0); addParameter(p, 'nstr', 8); addParameter(p, 'idatm', 2); parse(p, varargin{:}); opt = p.Results; % 使用临时文件名,避免同名冲突 input_file = sprintf('INPUT_%s.tmp', datestr(now, 'HHMMSS')); output_file = 'OUTPUT'; % 写输入文件 fid = fopen(input_file, 'w'); fprintf(fid, '$INPUT\n'); fprintf(fid, ' nstr = %d\n', opt.nstr); fprintf(fid, ' idatm = %d\n', opt.idatm); fprintf(fid, ' sza = %.3f\n', sza); fprintf(fid, ' phi = 0.0\n'); fprintf(fid, ' umu = %.3f\n', opt.umu); fprintf(fid, ' wlmmin = %.4f\n', wlmin); fprintf(fid, ' wlmmax = %.4f\n', wlmax); fprintf(fid, ' wlres = 0.005\n'); fprintf(fid, ' iaer = %d\n', opt.iaer); fprintf(fid, ' taer55 = %.3f\n', aod); fprintf(fid, ' baer55 = 0.90\n'); fprintf(fid, ' caer55 = 0.70\n'); fprintf(fid, ' alb = %.3f\n', opt.alb); fprintf(fid, ' iout = 1\n'); fprintf(fid, '$end\n'); fprintf(fid, '\n'); fclose(fid); % 运行 [status, cmdout] = system(sprintf('sbdart.exe < %s', input_file)); if status ~= 0 error('SBDART 执行失败: %s', cmdout); end % 读取输出文件 data = load(output_file); result.wavelength = data(:,1); result.radiance = data(:,2); result.status = status; % 清理临时输入文件 delete(input_file); end这个函数虽然简单,但有两个好处:一是所有临时文件名都带时间戳,避免多进程或重复运行时文件覆盖;二是返回的是结构体,后面画图、存数据都比较方便。你在实际项目里可以直接在这个基础上扩展,比如把iout参数也做成可选参数,让它既能算辐亮度又能算透过率。
4.2 编译 mex 接口的高级玩法
如果你追求的是执行效率,不满足于每次用system调 exe 的开销(虽然这个开销通常只有几十毫秒,但批量跑几百次也挺烦人),可以考虑用 mex 封装。原理是把 SBDART 的 Fortran 源码改造成一个可以被 Matlab 直接调用的函数,比如:
subroutine sbdart_mex(sza, aod, wlmin, wlmax, radiance, wavelength)然后用 mex 命令编译成sbdart_mex.mexw64。这样做的好处是省掉了文件 I/O 和进程创建的时间,缺点是得懂一点 Fortran 和 mex 的接口规则,而且调试起来比较费劲。
我个人的建议是:如果你只是做研究和论文仿真,用system调用就足够了,没必要折腾 mex;如果是开发一个需要反复调用的产品原型,再考虑 mex 封装。
5. 实操复现:一个从辐亮度计算到出图的完整例子
5.1 研究场景设定
我拿一个实际做过的场景来说明:假设你有一个可见光波段的传感器,想评估不同气溶胶负荷下大气顶(TOA)辐亮度的变化。这里的关键是想知道,当 AOD 从 0.1 升到 1.0 时,传感器接收到的信号会怎样变,进而评估气溶胶对植被指数计算的干扰。
物理上,这个问题的核心是辐射传输方程中大气散射和吸收对目标信号的调制:气溶胶增多会增强大气程辐射(path radiance),同时削弱地表反射信号,两个效应叠加,导致 TOA 辐亮度发生变化,且不同波段变化幅度不同。
我用中纬度夏季大气模型(idatm=2)、太阳天顶角 40 度、天底观测(umu=1.0),地表设为植被(alb=0.15)。波段设在 400 到 700 nm,步长 5 nm,分别计算 AOD = 0.1、0.3、0.5、1.0 四组。
5.2 逐行运行代码
先把上面的sb_run函数存成sb_run.m,然后新建一个脚本来批量跑:
% 清空环境 clear; clc; % 设置工作目录 cd('D:\projects\sbdart_matlab'); % 设定参数组合 sza = 40; aod_list = [0.1 0.3 0.5 1.0]; wlmin = 0.4; wlmax = 0.7; % 存储结果 results = cell(1, length(aod_list)); for k = 1:length(aod_list) fprintf('正在计算 AOD = %.1f ...\n', aod_list(k)); results{k} = sb_run(sza, aod_list(k), wlmin, wlmax); % 保存数据到文件,方便后面分析 T = table(results{k}.wavelength, results{k}.radiance, ... 'VariableNames', {'Wavelength_nm', 'Radiance'}); writetable(T, sprintf('TOA_AOD%.1f.csv', aod_list(k))); end跑完后,当前目录下会生成四个 CSV 文件,每个文件里都有波长和对应的 TOA 辐亮度值。
5.3 画图并分析结果
把四组数据画到一张图上观察趋势,这是我的标准操作:
figure('Color', 'w', 'Position', [100 100 800 600]); colors = lines(4); for k = 1:length(aod_list) plot(results{k}.wavelength * 1000, results{k}.radiance, ... 'Color', colors(k, :), 'LineWidth', 1.5); hold on; end xlabel('波长 (nm)'); ylabel('TOA 辐亮度 (W m^{-2} sr^{-1} μm^{-1})'); legend('AOD=0.1', 'AOD=0.3', 'AOD=0.5', 'AOD=1.0', 'Location', 'best'); grid on; title('不同AOD下大气顶辐亮度对比');你大概率会看到这样的现象:在蓝光波段(400-500 nm),AOD 增大导致辐亮度显著增加,因为气溶胶散射增强了程辐射;在红光波段(650-700 nm),变化相对较小,因为瑞利散射随波长的四次方衰减,气溶胶的散射也随波长降低。这正是气溶胶会“污染”植被指数这类比值型遥感指数的原因之一。
5.4 计算并验证植被指数偏差
为了量化这个影响,我接着计算归一化植被指数 NDVI。NDVI 的定义是:
NDVI = (ρ_nir - ρ_red) / (ρ_nir + ρ_red)在真实 SBDART 计算中,我们需要同时输出 NIR 和 RED 波段的反射率或者辐亮度。如果我用 650 nm 辐亮度代替红光,用 800 nm 辐亮度代替近红外,就可以模拟“不经过大气校正、直接拿 TOA 辐亮度算 NDVI”的效果。代码如下:
% 提取 650nm 和 800nm 附近的辐亮度值 lambda_target = [0.65, 0.80]; ndvi_toa = zeros(1, length(aod_list)); for k = 1:length(aod_list) [~, idx650] = min(abs(results{k}.wavelength - lambda_target(1))); [~, idx800] = min(abs(results{k}.wavelength - lambda_target(2))); r_red = results{k}.radiance(idx650); r_nir = results{k}.radiance(idx800); ndvi_toa(k) = (r_nir - r_red) / (r_nir + r_red); fprintf('AOD=%.1f NDVI(TOA)=%.4f\n', aod_list(k), ndvi_toa(k)); end你会发现 AOD 从 0.1 变化到 1.0 时,基于 TOA 辐亮度算出的 NDVI 出现系统性下降或上升(具体方向取决于波段响应关系),这就是气溶胶对植被遥感监测造成干扰的直接证据。你可以把这个结果作为论文中的敏感性分析图,也可以用来论证“为什么遥感反演前需要做大气校正”。
6. 常见问题与排查技巧实录
6.1 SBDART 运行报错“Error opening file”或“Cannot open INPUT”
这个错误我在换路径后经常遇到,原因基本是工作目录不对。system函数默认在当前 Matlab 工作目录下找sbdart.exe和INPUT文件,如果你改了目录而没切换,就会出现找不到文件的错误。
排查顺序:
- 先用
cd切到 SBDART 所在目录; - 确认
sbdart.exe和INPUT文件都在当前目录下; - 用
dir命令看看文件是否存在、权限是否正确; - 如果 SBDART 读取的文件名是固定的
INPUT,就不要改成INPUT_sza40这种名字去跑,否则程序找不到。
6.2 输出文件全是 NaN 或 0
如果你发现load出来的数据全是 0 或者 NaN,大概率是输入参数里的问题。常见原因如下:
umu设成了负值或大于 1,观测天顶角余玄不合理,导致几何计算出错,这时 SBDART 可能给了无效输出;sza大于 90 度,太阳在地平线以下,辐亮度理论上接近 0,但数值计算会出现奇怪结果;wlres设置太小,导致波长网格数过多,内存溢出或计算时间爆炸;在 Windows 32 位系统上尤其容易出问题。
检查思路是先把参数调回一个基准值(比如sza=30、umu=1.0、wlres=0.01),确认程序能跑出正常结果,再逐步改参数,每改一步就验证一次,定位问题出在哪个参数上。
6.3 同一次运行多次结果不一致
有时候你发现连续跑两次同样输入,结果却不同。这种情况我遇到过,原因多半是忘记备份输出文件,或者 SBDART 读取的环境变量有残留。但 SBDART 是确定性程序,输入相同理应结果相同,如果你的结果波动很大,建议检查:
- 是否有其他程序在同时读写
OUTPUT文件; - 是否用了多线程 Matlab,并行任务之间相互覆盖了输入/输出文件;
- 是否不小心改了参数但没改脚本里的显示值,导致“看起来一样、其实不一样”。
6.4 辐亮度单位换算的问题
SBDART 输出的辐亮度单位通常是W m^-2 sr^-1 μm^-1(即单位波长间隔的辐亮度),这是遥感常用的单位。但有些传感器输出的是mW cm^-2 sr^-1 μm^-1或者W m^-2 sr^-1 nm^-1,互相之间换算要小心。网上有经验丰富的同学发过换算技巧:1W m^-2 sr^-1 μm^-1等于 0.1mW cm^-2 sr^-1 μm^-1;如果按纳米归一化,需要除以 1000。单位写错会把结果差出几个数量级,这是最隐蔽的低级错误。
6.5 SBDART 计算速度慢怎么办
SBDART 本身比 MODTRAN 快很多,但如果设了非常细的波长分辨率和很高的nstr,计算量会成倍增加。我实测的经验数据如下:
| 参数组合 | 耗时(单次运行) |
|---|---|
| nstr=4, wlres=0.01, 0.4-0.7 μm | <1 秒 |
| nstr=8, wlres=0.005, 0.4-0.7 μm | 约 2-3 秒 |
| nstr=16, wlres=0.001, 0.4-2.5 μm | 几十秒 |
如果要做敏感性试验几百组,建议先用粗分辨率跑一遍,找到趋势后再对关键区间加密。不要一上来就全波段最高精度,否则时间全耗在计算上了。
6.6 常见错误速查表
我把实操中常见的错误整理成一张表,方便你对照排查:
| 错误现象 | 可能原因 | 解决方法 |
|---|---|---|
system返回 1,无输出 | 工作目录不对,或 exe 权限不足 | 检查 cd 目录和文件权限 |
| 输出文件全为 0 | umu或sza几何参数非法 | 检查几何参数范围 |
| 结果出现负值 | 波长范围包含大气吸收带且没有设置合适输出 | 检查iout类型,尝试iout=2 |
报错Segmentation fault | 某个参数数组越界(如nstr过大) | 降低nstr,检查版本兼容性 |
| 计算结果与文獻差较大 | 气溶胶模型或大气模型选择不当 | 检查idatm和iaer是否符合场景 |
| 读取 OUTPUT 时 textscan 失败 | Fortran 格式的空格和制表符混用 | 改用textscan(fid,'%f %f')读取 |
7. 更多应用场景:从辐亮度到多角度遥感模拟
SBDART 不仅能算单一的观测几何,还可以通过循环改变umu和phi来模拟多角度观测。这个功能在做 BRDF 敏感性分析或者多角度遥感器设计时非常有用。比如我把umu依次设为 1.0、0.8、0.6、0.4,对应观测天顶角 0、36.9、53.1、66.4 度,看不同观测角度下气溶胶的影响是否不同。这个计算如果用手写代码会很繁琐,有了sb_run封装函数,就是一个循环而已:
umu_list = [1.0 0.8 0.6 0.4]; for j = 1:length(umu_list) r = sb_run(40, 0.3, 0.4, 0.7, 'umu', umu_list(j)); plot(r.wavelength*1000, r.radiance); hold on; end你会发现,观测天顶角越大,大气路径越长,气溶胶散射贡献越明显。这个趋势在遥感大气校正常常被用来修正观测几何效应,也是很多气溶胶反演算法依赖于多角度信息的原因。
另外,SBDART 还可以用来做热红外波段的仿真。只要把波长范围改成 8-14 μm,大气模型选热带或中纬度,就能模拟地表温度反演中大气透过率的影响。虽然 SBDART 在热红外波段的精度不如 MODTRAN,但对于传感器通道选择和误差敏感性分析,精度完全够用。
8. 几个提高效率的小技巧
最后分享三个我在项目里用得很顺手的技巧,可能不常见但非常管用。
第一个是并行计算。Matlab 的parfor可以直接套在批量运行的循环外面,但要注意你的 SBDART 每次运行时输出的文件名是固定的,多个并行 worker 同时写OUTPUT文件会互相覆盖。解决方案是给每个 worker 分配不同的临时目录:
parfor k = 1:length(aod_list) workdir = sprintf('D:\projects\sbdart_matlab\tmp%d', k); mkdir(workdir); cd(workdir); % 拷贝 SBDART 可执行文件到临时目录 copyfile('D:\projects\sbdart_matlab\sbdart.exe', workdir); % 在这里执行 sb_run 逻辑 % ... end第二个是用 SBDART 做传感器波段响应加权。大多数遥感传感器的波谱响应不是矩形,而是类似高斯分布。SBDART 的isat参数内置了一些常见传感器响应函数,但如果你的传感器不在内置列表里,就需要自己计算加权平均。方法是在传感器响应范围内做细网格 SBDART 计算,然后按响应函数加权平均:
% 假设 sensor_response 是 400-700 nm 的响应向量 % sbdart_output 是相同波长范围的辐亮度向量 effective_radiance = trapz(wavelength, sbdart_output .* sensor_response) / trapz(wavelength, sensor_response);第三个建议是,记住 SBDART 的输入文件里$INPUT和$end之间的每个条目后面尽量不要留多余空格,虽然 Fortran 的 list-directed 读取比较宽容,但某些旧版本会因为多余字符报错。如果你的程序在别人机器上能跑、在你机器上报错,先检查一下是不是输入文件格式有细微差别。
我个人在实际操作中最深的体会是,辐射传输工具的价值不在于它有多先进,而在于你有多理解它背后的物理假设。SBDART 的平面平行大气假设决定了它在大天顶角、水平非均匀场景下的局限性,但只要你清楚这个边界,在遥感影像仿真、传感器设计论证、气溶胶影响评估这类场景里,它完全够用。而且因为源码开放,你可以随时翻看某个物理过程到底是怎么算的——这一点是商业黑盒软件给不了的。最后再分享一个小建议:跑任何一组参数之前,先跑一个已知结果的基准算例(比如无气溶胶的瑞利散射场景),确认自己的环境没问题,再开始批量计算,这能避免你把所有时间都耗在调试一个低级错误上。
本文还有配套的精品资源,点击获取