简介基于MATLAB的二维辐射传输方程RTE数值求解源码包面向科学计算与物理模拟方向的研究者、工程师及高年级学生。项目涵盖二维RTE数学建模、散射与吸收介质参数定义、离散坐标法等数值解法实现并附有可运行的主程序与几何示例便于快速理解RTE从方程到离散求解的完整流程。压缩包共30个文件包括22个h头文件、3个m脚本、2个c源文件、2个mexw32编译接口及1个mat示例几何文件整体仅75KB结构紧凑适合直接研读或在MATLAB中扩展应用。已有915人学习或浏览小众但实用尤其对大气辐射、光学成像或传热模拟等场景的建模工作有直接参考价值。通过阅读源码可掌握RTE求解的编程实现技巧并为自定义散射、吸收工况提供可复用框架。1. 为什么是二维RTE从标量辐射传输到MATLAB求解的最小闭环辐射传输方程RTE描述光子在介质中传播时被吸收、散射和发射的过程是遥感、生物组织成像、燃烧诊断和反应堆屏蔽的通用数学模型。一维RTE只沿深度方向演化很多实际问题并不满足这个假设当介质在横向存在强不均匀性时必须把空间拓展到二维。二维RTE把空间坐标限制在(x,y)平面方向角保留在半球或全球立体角内比三维模型少一个自由度却依然能捕捉侧面散射和边界漏射这类关键效应。这里有一个初学者容易反直觉的点二维RTE的数值求解用MATLAB写出来并不需要上千行框架代码几百行就足够跑通但大多数实现都败在角度离散与边界条件交叉处理上而不是方程本身。这篇内容把离散坐标法DOM求解二维RTE的完整路径拆开从方向求积、网格构造、源项迭代到参数调优和验证适合正在写辐射仿真程序、想验证自研代码正确性的工程师和研究者。2. 离散坐标法的核心方程与网格构造2.1 二维RTE控制方程与方向离散对于不含内发射源的稳态辐射传输在二维直角坐标下光谱辐射强度I(x,y,Ω)沿方向Ω在平面内的方向余弦为μ和η满足如下方程μ ∂I/∂x η ∂I/∂y -β_t I σ_s/4π ∫_4π I(x,y,Ω) Φ(Ω,Ω) dΩ其中β_t是消光系数σ_s是散射系数Φ为散射相函数且归一化为在4π立体角上的积分为4π。方程左侧是沿方向的方向导数右侧第一项为吸收和散射造成的衰减第二项为来自其它方向散射的贡献。因为角度变量连续直接求解这个积分微分方程几乎不可能标准做法是用离散坐标法把角度方向离散成有限个方向(μ_m, η_m)每个方向配一个权重w_m。这样角度积分变成加权求和μ_m ∂I_m/∂x η_m ∂I_m/∂y β_t I_m σ_s/4π Σ_m w_m Φ_m,m I_mS_N求积组是应用最广的方向集合。二维平面几何下S_2包含4个方向S_4包含12个方向S_8包含40个方向每个方向都是一对方向余弦和对应的权重。这套命名方式中的N是某方向积分点的数目而不是方向总数。选择阶数时我一般在工程计算里先用S_4跑通再与S_8结果对比如果相对误差在容许范围内就不会继续往高走。S_8以上权重会出现负值容易让迭代产生非物理振荡这不是越高越好。下面给出一个S_2求积组在MATLAB中的表示例子% S_2 求积组对应二维平面几何共4个方向 mu [ 0.57735027, 0.57735027, -0.57735027, -0.57735027]; eta [ 0.57735027, -0.57735027, 0.57735027, -0.57735027]; w pi * ones(4,1); % 每个方向权重为π总和为4π注意这里的权重和为4π对应立体角积分。方向余弦满足μ²η²≤1因为第三维方向余弦不为零时不在平面几何的简化范围内二维平面几何实际上对应无限长柱体沿z方向均匀的情况方向角只在x-y平面内变化。上面的代码里μ和η的绝对值都是1/√3平方和为2/3剩余部分在z方向这种简化是二维RTE特有的。2.2 有限体积法空间网格构造空间离散我通常用有限体积法。计算域划分成Nx×Ny个控制体控制体中心坐标为(x_i, y_j)。对方程在控制体上积分并用迎风差分近似界面上的强度。以μ_m 0且η_m 0的方向为例信息从左边界和下边界进入因此左界面的上游强度取I_{i-1,j}下界面取I_{i,j-1}。离散后得到代数方程μ_m(I_{i,j}^m - I_{i-1,j}^m)/Δx η_m(I_{i,j}^m - I_{i,j-1}^m)/Δy β_t I_{i,j}^m S_{i,j}^m整理成显式计算式I_{i,j}^m ( |μ_m| I_up_x/Δx |η_m| I_up_y/Δy S_{i,j}^m ) / ( |μ_m|/Δx |η_m|/Δy β_t )这里的S_{i,j}^m是散射源项。这个式子就是扫描法的核心。在MATLAB里网格用两个列向量x和y表示物理参数用Nx×Ny的矩阵存储每个方向的方向余弦是标量不随空间改变因此不需要为每个方向分配一套方向矩阵。下面这段代码生成了介质参数场Nx 50; Ny 50; Lx 1.0; Ly 1.0; x linspace(0, Lx, Nx); y linspace(0, Ly, Ny); sigma_t 2.0 * ones(Nx, Ny); % 消光系数矩阵 sigma_a 0.5 * ones(Nx, Ny); % 吸收系数矩阵 sigma_s sigma_t - sigma_a; % 散射系数矩阵 % 绘制散射系数分布确认介质设定 imagesc(x, y, sigma_s); axis xy; colorbar; xlabel(x (m)); ylabel(y (m)); title(散射系数分布);sigma_t、sigma_a、sigma_s都是Nx×Ny的双精度矩阵后续计算直接逐元素参与运算。注意imagesc默认把第一维画在纵轴转置加axis xy可以让图像方向与空间坐标一致。网格不需要在内存里构建二维坐标网格只要知道Δx和Δy扫描时用数组下标即可。2.3 求积组阶数选取与权重归一化求积组的选择直接影响解的精度和稳定性。S_4有12个方向对所有方向权重均为正多数各向同性散射问题够用。S_6和S_8虽然能提高角度分辨但会出现负权重负权重会让散射源项迭代过程出现负强度严重时直接发散。如果必须用高阶梯数建议改用没有负权重的乘积型求积组或者在每次迭代后把负强度强制归零但后者会破坏能量守恒只能作为临时手段。下面给出一个权重检查函数用来确认归一化是否正确function ok check_quad(mu, eta, w) total sum(w); if abs(total - 4*pi) 1e-12 warning(权重之和为%.15f不等于4π, total); ok false; else ok true; end end对S_2、S_4这类标准求积组权重之和应当严格等于4π。如果自定义方向必须先做这个检查否则源项迭代中的系数σ_s/4π会出现整体偏差导致结果系统性偏低。3. 用MATLAB实现2D RTE从离散坐标到矩阵迭代求解3.1 源项迭代法为什么它能收敛离散坐标法把RTE变成一组耦合方程方向与方向之间通过散射源项互相联系。最常用的解算器是源项迭代法也叫“外迭代”。它从猜测的初始强度场出发固定散射源项逐个方向做扫描求解再用新得到的强度更新源项如此反复直到残差足够小。这个方法思路直接适合MATLAB矢量化。每次迭代中每个方向的扫描实际上在解一个没有方向耦合的稳态对流方程因此天然稳定。当散射占主导时源项迭代收敛速度会变得很慢。原因是方程趋近于扩散方程散射源项更新相当于进行一次雅可比迭代谱半径接近1。这个问题在节4.1会单独讨论加速手段基础版本对工程计算足够。3.2 主循环框架代码下面的函数实现了一个二维各向同性散射的RTE求解器骨架。它接收光学参数、方向集和边界条件返回强度场。function [I, iter] solve_rte_2d(sigma_t, sigma_s, mu, eta, w, bc, dx, dy, maxiter, tol) % sigma_t: Nx*Ny 消光系数矩阵 % sigma_s: Nx*Ny 散射系数矩阵 % mu, eta: 每个方向的方向余弦长度 M % w: 方向权重长度 M % bc: 边界条件结构体含 left/right/bottom/top 四个边界 % 每个边界是 M*N 矩阵M方向数N该边界上的空间节点数 % dx, dy: x和y方向的网格间距 % 返回 I: Nx*Ny*M 强度场 [Nx, Ny] size(sigma_t); M length(mu); I zeros(Nx, Ny, M); I_old zeros(Nx, Ny, M); for iter 1:maxiter % 1. 由当前强度场计算散射源项各向同性散射 weighted sum(I .* reshape(w, 1, 1, M), 3); % 对所有方向加权求和 S sigma_s ./ (4*pi) .* weighted; % 源项与方向无关 S repmat(S, 1, 1, M); % 复制到每个方向 % 2. 对所有方向扫描 for m 1:M I_dir I(:,:,m); S_dir S(:,:,m); % 确定扫描起点和方向 ix 1:Nx; iy 1:Ny; if mu(m) 0, ix Nx:-1:1; end if eta(m) 0, iy Ny:-1:1; end % 初始化上游值边界条件 up_x zeros(Nx, Ny); up_y zeros(Nx, Ny); if mu(m) 0 up_x(1,:) bc.left(m,:); % 左边界入射 else up_x(Nx,:) bc.right(m,:); % 右边界入射 end if eta(m) 0 up_y(:,1) bc.bottom(m,:); % 下边界入射 else up_y(:,Ny) bc.top(m,:); % 上边界入射 end % 用迎风差分顺序更新内部节点 for i ix for j iy if i 1 i Nx up_x(i,j) I_dir(i - sign(mu(m)), j); elseif i 1 mu(m) 0 up_x(i,j) I_dir(i1, j); % 反向扫描内部上游 elseif i Nx mu(m) 0 up_x(i,j) I_dir(i-1, j); end if j 1 j Ny up_y(i,j) I_dir(i, j - sign(eta(m))); elseif j 1 eta(m) 0 up_y(i,j) I_dir(i, j1); elseif j Ny eta(m) 0 up_y(i,j) I_dir(i, j-1); end % 更新当前点强度 abs_mu abs(mu(m)); abs_eta abs(eta(m)); coeff abs_mu/dx abs_eta/dy sigma_t(i,j); I_dir(i,j) (abs_mu*up_x(i,j)/dx abs_eta*up_y(i,j)/dy S_dir(i,j)) / coeff; end end I(:,:,m) I_dir; end % 3. 收敛判断 diff max(abs(I - I_old), [], all); I_old I; if diff tol break; end end end代码逻辑分为三个部分先用上一轮的强度场计算每个空间点的加权总强度乘以σ_s/(4π)得到散射源项然后对每个方向根据方向余弦的符号选择扫描起点用迎风差分沿着空间网格顺序更新最后比较新旧强度场的最大绝对差值判断是否收敛。参数说明maxiter是允许的最大迭代次数tol是收敛阈值一般取1e-6即可。这个框架中的边界条件结构体bc每个边界矩阵的行数等于方向数列数等于该边界上的网格节点数。注意对于从边界出发的方向边界上游值在扫描前就确定了对于到达边界的方向边界上计算出的强度自然流出不需要额外处理。这个版本的性能不高因为内层有双循环。实际使用时我会先把边界和上游值的初始化做好再把双循环替换成矩阵切片运算这在节5.2展开。3.3 边界条件初始化与常见错误边界条件在RTE求解中比想象中更容易写错。对于左边界只有μ_m 0的方向需要给定入射强度μ_m 0的方向强度是从计算域内部扫出来的不能提前赋值。同理右边界只有μ_m 0方向需要给定入射强度上边界对应η_m 0方向下边界对应η_m 0方向。构造边界时最常见的错误就是把所有方向的边界强度都赋成常数这会改变能量平衡。下面是一个准直光束入射的边界设置示例M length(mu); bc.left zeros(M, Ny); bc.right zeros(M, Ny); bc.bottom zeros(M, Nx); bc.top zeros(M, Nx); % 左侧入射准直光束方向为 (mu1, eta0) 附近 beam_mu 1.0; beam_eta 0.0; beam_intensity 1.0; for m 1:M if abs(mu(m) - beam_mu) 1e-4 abs(eta(m) - beam_eta) 1e-4 bc.left(m, :) beam_intensity; fprintf(找到准直方向m%d\n, m); end end代码里beam_mu1, beam_eta0对应与x轴重合的入射方向在这个方向上左侧边界各节点强度均为1。如果使用的求积组没有正好包含该方向需要改用带奇异方向修正的方法比如把准直方向单独处理不参与求积。这里不展开修正实际工程中更稳健的做法是采用“有限束”近似将准直方向用一个极窄的高斯锥代替。4. 参数选择与边界条件散射反照率、光学厚度、求积阶数的坑4.1 散射反照率与慢收敛散射反照率定义为ω σ_s / β_t范围在0到1之间。ω接近0时介质以吸收为主源项迭代很快收敛通常不超过10次ω接近1时散射主导光子可以经历大量散射事件源项迭代的收敛就变得极其缓慢。举个例子ω0.99时可能需要迭代上千次才能达到1e-6的残差这在MATLAB里会让人直接失去耐心。提示当ω 0.9时不要无脑增加迭代次数应当先启用扩散合成加速或采用更快的迭代格式。一个相对容易实现的加速手段是给散射源项加一个松弛因子即把新计算的源项与旧源项做加权平均alpha 0.7; % 松弛系数0.5~0.9之间调试 S_new_calc sigma_s ./ (4*pi) .* weighted; S alpha * S_new_calc (1 - alpha) * S_old;松弛因子α小时更稳定但更慢α大时可能振荡。我一般先从0.7开始试观察残差曲线如果出现振荡就调小。这个方法虽然简单但对ω在0.9~0.95区间有明显的加速效果。更高ω时建议改用DSA不过代码复杂度会上升不少。4.2 光学厚度与空间网格约束光学厚度τ β_t × L决定了介质对辐射的总作用。τ远小于1时介质光学薄散射和吸收作用弱辐射基本按几何光学直行低阶求积和较粗网格就能获得好结果。τ大于10时属于光学厚介质平均自由程很小必须细化网格。经验准则是Δx和Δy都小于1/β_t的三分之一。对于光学厚介质迎风差分带来的数值耗散很严重。例如β_t 10 m^-1L 1 mΔx若取0.05 m则每个平均自由程只有2个网格数值结果会出现明显的“数值泄漏”光斑边缘变糊。此时应把网格加密到Nx200以上或者使用二阶格式。4.3 求积阶数对结果的影响下面这份经验表可以帮助初步选型求积阶数方向数权重符号适用场景S_24全正快速测试、光学薄介质S_412全正各向同性散射、工程常规S_624个别负权重适度各向异性、光厚介质S_840负权重常见强各向异性、高精度对比负权重并不是不能用的洪水猛兽但如果迭代过程中出现负强度且持续振荡多半就是负权重导致的。我通常用S_4求解用S_8校验若两者最大相对误差小于1%就以S_4结果作为最终解。这个“低阶求解高阶验证”的双跑流程比一开始就上S_8要高效得多。4.4 边界条件与能量守恒校验边界条件类型决定了系统与外界交换能量的方式。二维RTE的常见边界可归纳为三类真空边界所有入射方向强度均为0模拟开放空间。准直入射指定某个方向的强度为非零值模拟激光或太阳光。反射边界引入反射率r反射方向的强度由入射方向的强度转换而来。设置反射边界时要注意方向配对。镜面反射遵循“入射角等于反射角”在离散方向坐标中反映为μ和η变号。例如在左侧边界发生镜面反射出射方向为(μ, η)的边界强度应等于入射方向(-μ, η)的强度乘以反射率。如果漏掉方向映射边界会变成纯吸收边界能量不守恒。进行能量守恒校验时我一般计算所有边界的净通量% 计算各边界出射通量示意 left_flux_in sum(sum(bc.left .* mu_pos)); % 实际需按方向加权 % ...辐射通量等于边界上的强度乘以该方向的方向余弦再乘以权重对角离散公式求和。若总入射通量与内部吸收功率之差超过5%优先排查权重归一化再检查边界方向是否配齐。4.5 一个可直接修改的参数模板结合上述讨论下面是一个完整算例的参数设置可以直接替换前面的求解器来运行Lx 0.5; Ly 0.5; Nx 50; Ny 50; dx Lx/(Nx-1); dy Ly/(Ny-1); beta_t 5.0; sigma_a 1.0; sigma_s 4.0; % ω 0.8 % 生成光学参数矩阵 sigma_t_mat beta_t * ones(Nx, Ny); sigma_s_mat sigma_s * ones(Nx, Ny); % 方向与权重使用S_4 [S_mu, S_eta, S_w] get_quad_S4(); % 该函数返回12个方向 % 边界条件 M length(S_mu); bc.left zeros(M, Ny); bc.right zeros(M, Ny); bc.bottom zeros(M, Nx); bc.top zeros(M, Nx); % 左侧入射准直光束方向近似(1,0) for m 1:M if abs(S_mu(m)-1) 1e-3 abs(S_eta(m)) 1e-3 bc.left(m,:) 1.0; end end % 调用求解器 [I, iter] solve_rte_2d(sigma_t_mat, sigma_s_mat, S_mu, S_eta, S_w, bc, dx, dy, 100, 1e-6); fprintf(迭代次数%d\n, iter);这段代码中的get_quad_S4()需要自己实现或从资料中获取。参数模板的好处是修改介质分布或边界条件时不需要改动求解器只用改矩阵和结构体。5. 验证与加速技巧用解析解校验、向量化与并行5.1 用纯吸收介质解析解验证程序正确性写完求解器后第一件事不是跑复杂工况而是把散射系数设为0让方程退化为纯吸收问题。此时RTE变成一组独立的单向传输方程沿光束方向强度按朗伯-比尔定律指数衰减。取左侧中心一条水平线理论解为I(x) I0 exp(-β_t x)。与数值解对比可以直接定位扫描方向或边界条件的错误。% 纯吸收验证sigma_s 0 sigma_s zeros(Nx, Ny); [I, ~] solve_rte_2d(sigma_t_mat, sigma_s, S_mu, S_eta, S_w, bc, dx, dy, 1, 1e-6); % 取中心水平线上的强度找到最接近水平方向的m [~, mid] min(abs(S_mu - 1) abs(S_eta)); I_num squeeze(I(:, Ny/2, mid)); x_line x(:); I_exact 1.0 * exp(-beta_t * x_line); err norm(I_num - I_exact) / norm(I_exact); fprintf(纯吸收验证相对误差%e\n, err);如果误差在1e-2以下说明扫描顺序和边界方向基本正确若误差达到0.5以上先检查mid对应的方向是否真正沿x轴很多S_4方向并不包含μ1此时需要把准直方向单独处理或选用更细方向集。5.2 把空间双循环改成矩阵切片在迭代循环中空间扫描的双循环是瓶颈。以方向μ0、η0为例内部节点(i2:Nx, j2:Ny)的更新可以完全用数组运算完成% 假设mu0, eta0 I_dir(2:Nx, 2:Ny) (abs(mu)/dx * I_dir(1:Nx-1, 2:Ny) ... abs(eta)/dy * I_dir(2:Nx, 1:Ny-1) ... S_dir(2:Nx, 2:Ny)) ./ ... (abs(mu)/dx abs(eta)/dy sigma_t(2:Nx, 2:Ny));边界上的点单独处理。这样空间双循环被完全消除MATLAB底层的矩阵运算速度比逐元素循环快一个数量级。注意I_dir(2:Nx, 1:Ny-1)得到的是下边界的内部点前提是η0若η0则切片方向要反过来。对于不同符号组合只需要把切片方向和顺序对应调整。我通常写四个分支μ±, η±每个分支内都做矩阵切片这样可以避免在循环内判断符号效率最高。5.3 方向循环用parfor并发源项迭代法的每个方向扫描使用旧迭代步的源项因此不同方向之间没有数据依赖。可以把方向循环改为parfor但必须把每次迭代的中间结果先存入新的数组防止并行写入冲突。parfor m 1:M I_dir I_old(:,:,m); % 读旧值 % ... 执行扫描更新得到I_dir I_new(:,:,m) I_dir; % 写入新数组 end I I_new;使用parfor需要提前把I_new初始化为与I同样大小的数组并且保证每个m只访问I_old(:,:,m)。在本地多核机器上这个改动通常能把单次迭代时间缩小到原来的1/3到1/2。如果机器有GPU也可以把I声明为gpuArray让矩阵切片自动在GPU上执行但网格太小时GPU传输开销会掩盖收益建议在Nx×Ny≥2500时再考虑。5.4 用profile找出真正的热点对求解器做性能分析时我习惯在启动求解前加上profile on跑完一次后执行profile viewer。热点几乎总是集中在散射源项计算和方向扫描中。对于前者sum(I .* reshape(w,1,1,M), 3)已经是比较高效的方式如果还嫌慢可以改用pagemtimes或把三维数组展开成二维矩阵做矩阵乘法。对于后者优先确认是否已经使用矩阵切片而不是把时间花在循环提速上。一个实用的长期做法是把误差计算和绘图放进同一个脚本每当调整求积阶数或网格尺寸后就画出强度分布并记录收敛曲线。通过对比S_4和S_8的结果能快速找出离散误差的主要来源再决定是加密网格还是提高方向阶数这样调试路径比盲目尝试要短得多。本文还有配套的精品资源点击获取