11b - 上下行辐照度图 (IrradianceMap)¶

本章介绍 IrradianceMapSensor —— 以像元栅格记录场景每个位置的下行与上行辐射通量。这一功能对应 LESS 经典版本的「上下行辐射产品」，常用于复杂地形局部太阳辐射分析、林下 PAR 分布、以及地-冠层通量建模。

物理含义与记录规则¶

IrradianceMapSensor 在场景 XY 平面上划分均匀栅格，每个栅格称为一个像元 (pixel)。每个像元各维持两个累加器：

下行 ↓：进入该像元并发生交互的入射通量；
上行 ↑：从前一次命中像元出射、跨越到下一个像元的通量。

记录规则（前向光子跟踪）：

IrradianceMap 示意图

每次光子命中 \(P\)，在 \(P\) 所属像元 ↓ 累加当前通量。
若前一次命中 \(P_{\text{prev}}\) 与 \(P\) 属于不同像元，则在 \(P_{\text{prev}}\) 所属像元 ↑ 累加同一份通量。
若 \(P\) 与 \(P_{\text{prev}}\) 属于同一像元（像元内散射），不记录（既不算 ↓ 也不算 ↑）。
光子逃离场景时，若此前至少命中过一次，在最后一次命中所在像元 ↑ 累加当时通量。

这一规则是 LESS2 经典前向光子跟踪的通量约定 (photonrt_proc.cpp)； IrradianceMapSensor 在 LESS 的 OptiX、Vulkan 和 Embree 路径上使用相同的产品定义。

浑浊介质如何跨像元¶

TurbidBoundary 可以横跨任意数量的 IrradianceMap 像元。边界三角形只负责标记光线进入或离开介质，不会被记录成辐射交互。真正的统计叶片散射发生在介质内部，散射点落在哪个像元，该次下行通量就归入哪个像元。

当光子从一个交互点传播到另一个像元时，离开前一个像元柱的通量记入前一像元的 upwelling。如果两次交互仍在同一像元内，则没有穿过像元柱边界，不额外记录上行量。因此：

改变网格分辨率只重新分配空间位置，下行总能量以及直射、漫射总量保持守恒；
上行图描述穿过像元柱边界的通量，细网格拥有更多内部边界，其空间总和可以与粗网格不同；
有限场景先裁剪越界的 Turbid 体积，周期场景先拆分到基础瓦片，跨场景边界不会产生双重边界壳；
多个 Turbid 介质重叠时，消光系数相加，实际散射事件只归属于被抽中的介质和所在像元。

单位与面积归一化¶

product.downwelling / upwelling 存的是 W / nm 的像元柱子总通量（等于 水平 irradiance × 像元 XY 面积）。要得到 W/m²/nm（面积归一化后的 irradiance），调用 product 的 _irradiance() 方法指定除以哪种面积：

`surface=`	除数	物理含义	数据来源
`'horizontal'`（默认）	`dx · dy`（像元 XY 面积）	水平板上的 irradiance；"太阳辐射图" 的标准	不需要额外数据
`'slope'`	`dx · dy / cos(β)`	贴在倾斜 terrain 上的 pyranometer 读数	`terrain.dem` 自动求梯度；或用户传 `slope` 弧度图
`'surface'`	∑ DEM 三角面 3D 总面积 / 像元	精确 DEM 表面积归一化	自动从 DEM 2-三角网累加；或用户传 `surface_area`

p.downwelling_irradiance()                   # 默认 horizontal, W/m²/nm
p.downwelling_irradiance('slope')            # 斜面；自动从 terrain.dem
p.downwelling_irradiance('slope', slope=my_slope_rad)   # 用户传坡度 (rad)
p.downwelling_irradiance('surface')          # DEM 真实 3D 面积
p.upwelling_irradiance(...)                  # 同样 3 种 surface 选项

宽带积分（光谱 + 面积归一化一步到位）：

par = p.downwelling_integrated((400, 700), surface='slope')   # W/m² on slope

Product 还暴露两个从 DEM 预计算的辅助字段（平地时为 None）：

p.terrain_slope — [ny, nx] 弧度
p.terrain_surface_area — [ny, nx] m²（DEM 三角面精确累加）

基本 API¶

import numpy as np
import less

scene = less.Scene()            # 自动选择；也可显式使用 optix / vulkan / embree
scene.size = 40.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(source=less.Sun(zenith=30, azimuth=180), atmosphere=less.NoAtmosphere())
sensor = less.IrradianceMapSensor(
    bands=[550, 850],           # 波段 (nm)
    grid=(40, 40),              # 方式 A: 显式指定像元数
    # 或
    # resolution=0.5,           # 方式 B: 米/像元（由 scene.size 推导 grid）
    photons_per_m2=400,         # 每 m² 场景投射的光子数
    # num_photons=200_000,      # 可选：显式总光子数，设置后覆盖 photons_per_m2
    max_bounces=None,           # 记录深度上限；None=无限（仅 max_depth/RR 限制）
)
product = scene.simulate(sensor)

# 输出
product.downwelling                         # [ny, nx, nb]  W/nm per pixel
product.upwelling                           # [ny, nx, nb]  W/nm per pixel
product.downwelling_irradiance()            # [ny, nx, nb]  W/m²/nm (horizontal)
product.wavelengths                         # [nb]
product.cell_size                           # (dx, dy) m

参数¶

参数	含义	典型值
`bands`	波段 (nm) 列表或 `SpectralBands` 对象	`[550, 850]`
`grid`	水平像元数 `(nx, ny)`	`(100, 100)`
`resolution`	米/像元（与 `grid` 二选一，`resolution` 优先）	`0.5`
`photons_per_m2`	每 m² 场景的光子数	`400–4000`
`num_photons`	可选的显式总光子数；设置后覆盖 `photons_per_m2`	`200_000`
`max_bounces`	记录的最大命中次数（`None` = 不限）	`None` 或 `1`

grid vs resolution - 给 resolution：像元边长固定为物理尺度（如 0.5 m/pixel），grid 参数被忽略；nx = round(scene_width / resolution)。 - 给 grid：显式像元数。像元边长 = scene.size / nx。 - 两者都不给：默认 grid=(100, 100)。

max_bounces 与 LESS2 的 depth 的关系 - max_bounces=1 —— 仅记录第一次交互（直接太阳/天空辐射到达像元），等价于 LESS2 的 depth == 1。 - max_bounces=k —— 前 \(k\) 次交互参与 splat；再多 1 次仍追踪以捕捉逃逸上行。 - max_bounces=None —— 完整多次散射，直到 Russian roulette 或 max_depth 结束。

完整示例：平地朗伯地表¶

import math
import numpy as np
import less

# 场景：平地 + Lambertian 反射率 0.15，纯直射光（无天空）
scene = less.Scene()
scene.size = 40.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))

# 光谱参数按本次模拟使用的波段准备
# 真空中 Sun.irradiance 是 TOA 束法向辐照度；水平面通量还要乘 cos(天顶角)
horizontal = np.array([1.6203, 0.719], dtype=np.float32)
scene.illumination = less.Illumination(
    source=less.Sun(
        zenith=30, azimuth=180, wavelengths=[550, 850],
        irradiance=horizontal / math.cos(math.radians(30))),
    atmosphere=less.NoAtmosphere(),
)

sensor = less.IrradianceMapSensor(
    bands=[550, 850], resolution=1.0,   # 1 m/pixel → 40×40
    photons_per_m2=400, max_bounces=None,
)
p = scene.simulate(sensor)

# 解析检验：由 TOA 束法向辐照度投影得到水平面直射；上行 = REFL × 该值
E_expected = np.array([1.6203, 0.7190])
print("horizontal down (W/m²/nm):",
      p.downwelling_irradiance().mean(axis=(0, 1)))     # ≈ E_expected
print("horizontal up   (W/m²/nm):",
      p.upwelling_irradiance().mean(axis=(0, 1)))       # ≈ 0.15 × E_expected

配套脚本：scripts/11b_irradiance_map.py

DEM 上的 albedo 图¶

在起伏 terrain 上运行 IrradianceMapSensor，可以观察到强烈的坡向阴影效应；并能通过逐像元 up / down 算出 albedo 图来核验物理一致性（均匀 Lambertian 下各像元 albedo ≡ ρ）。

import numpy as np
import less
from scipy.ndimage import gaussian_filter

# 合成一块起伏 DEM（也可以传 GeoTIFF 路径）
def make_dem(size_m, n=128, seed=42):
    rng = np.random.default_rng(seed)
    dem = np.zeros((n, n), dtype=np.float64)
    for sigma, amp in [(18, 6), (9, 3), (4, 1.5), (2, 0.5)]:
        dem += amp * gaussian_filter(rng.standard_normal((n, n)), sigma=sigma)
    dem -= dem.min()
    dem *= 12.0 / dem.max()
    return dem.astype(np.float32)

scene = less.Scene()
scene.size = 100.0
scene.repetitive = False                     # 有限地形：silhouette 发射自动补足斜射边缘
scene.terrain = less.Terrain(
    property=less.Lambertian(reflectance=0.30),
    dem=make_dem(100.0),
)
scene.illumination = less.Illumination(source=less.Sun(zenith=45, azimuth=180), atmosphere=less.PrescribedAtmosphere(direct_beam_transmittance=1.0 - (0.0), diffuse_horizontal_transmittance=0.0))
sensor = less.IrradianceMapSensor(
    bands=[550], resolution=0.78,             # ≈ DEM 分辨率
    photons_per_m2=2000, max_bounces=None)
p = scene.simulate(sensor)

# 三种面积归一化
down_h = p.downwelling_irradiance('horizontal')[:, :, 0]      # W/m² 水平
down_s = p.downwelling_irradiance('slope')[:, :, 0]           # W/m² 斜面
down_a = p.downwelling_irradiance('surface')[:, :, 0]         # W/m² DEM 3D

# Albedo = up/down（比值与面积约定无关；逐像元应 ≈ rho = 0.30）
up_h   = p.upwelling_irradiance('horizontal')[:, :, 0]
albedo = up_h / np.maximum(down_h, 1e-10)
print("albedo mean:", np.nanmean(albedo))    # ≈ 0.30

验证脚本：tests/validation/val16_terrain_albedo.py

观察： - 下行 down_h 的最大值（阳坡）> 最小值（阴坡）常达 500× 差距（阴影区降到千分之几）。 - down_h（水平面）最亮；down_s（斜面 = down_h × cos β）稍暗；down_a（DEM 3D 面积）最小且最准。 - 逐像元 albedo 在三种定义下都严格等于 ρ（splat 规则的数学守恒：up = ρ × down 逐像元成立）。

后端选择¶

IrradianceMapSensor 支持三个后端：

scene = less.Scene(backend='optix')   # NVIDIA OptiX
scene = less.Scene(backend='vulkan')  # NVIDIA / AMD / Intel GPU
scene = less.Scene(backend='embree')  # 多线程 CPU

三者使用相同的产品定义和能量记账规则，结果应在解析容差或蒙特卡罗噪声范围内一致。OptiX 是 NVIDIA 上的生产首选；Vulkan 提供跨厂商 GPU 路径但仍属于 Experimental；Embree 无需 GPU。运行前可用 scene.can_use("irradiance_map") 获得当前环境下的单一可用性结论。

与 `RadiationFieldSensor` 的区别¶

特性	`IrradianceMapSensor`	`RadiationFieldSensor`
输出维度	2D (nx × ny)	3D (nx × ny × nz)
原始单位	W/nm per pixel column	W/m²/nm per 体素
记录对象	每次命中的入射/出射通量	每个体素的吸收能量
典型应用	上/下行辐射、地形辐射	fPAR、垂直辐射廓线

两者均为前向光子跟踪产品，可按需选用。

DEM 与 sensor 分辨率¶

DEM 与 scene.size 的关系：DEM（numpy 数组或 GeoTIFF）被严格拉伸到 scene.size 对应的 XY 边界（顶点在 np.linspace(0, scene_w, n_dem_x) 上）。从 GeoTIFF 读取的 x_origin / pixel_size 会被读但丢弃不用于场景定位。因此一个 (128, 128) 的 DEM 在 100 m scene 上，DEM 顶点间距 ≈ 100/127 = 0.787 m。

sensor resolution 与 DEM 分辨率解耦：你可以用任意 sensor resolution（或 grid）观测同一块 DEM。当两者不一致时，_derive_terrain_geometry 内部用双线性插值把坡度和三角面积映射到 sensor 栅格。

常见问题¶

Q: 上行累加会不会超过下行？ 不会。每一次散射都伴随 albedo 衰减（throughput × ρ），随散射次数增加，上行总能量永远收敛到不超过下行总量的数值。

Q: 为什么「像元内散射」不记录？ 避免同一个像元「自己入、自己出」的双记。对把像元当作封闭辐射收支单元的用户，这个规则使 down − up 等于像元内的净吸收。

Q: max_bounces=1 和无限次散射在平地 Lambertian 下为何几乎相同？ 平地朗伯地表没有能再次拦截光子的几何结构，第一次散射后光子基本直接逃逸。有冠层或起伏的场景，多次散射差异会明显。

Q: 均匀 Lambertian 下，为什么 albedo 图处处严格等于 ρ 而不是偏高？ 因为 splat 规则在每一次散射事件都配对记录 down 与 up（up = down × ρ 逐事件成立），无论场景几何如何，逐像元求和后 up/down = ρ 严格成立。要看到 albedo 显著偏离 ρ 的场景，需要非均匀反射率（例如不同坡向不同材料）或考察场景级（∑up / ∑down，会因山谷 radiation trapping 低于 ρ）。

下一步¶

12 - 热红外模拟