23 - 上下行辐照度图¶
IrradianceMapSensor 把场景在水平方向划分为规则网格，记录每个格子中向下和向上的辐射通量。可用于分析复杂地形上的太阳辐射、林下光合有效辐射的分布，以及地表与冠层之间的辐射交换。
记录方式¶
光子从光源出发在场景中传播。每个格子（称为像元）对应一个竖直的柱体，各有一个下行和一个上行累加量：
- 光子每次与场景发生作用，把当时的通量计入该处像元的下行量；
- 光子从一个像元的作用点传到另一个像元时，把通量计入前一个像元的上行量；同一像元内的连续作用不计上行；
- 光子最终离开场景时，把通量计入最后一次作用所在像元的上行量。
对统计植被，统计叶片上的散射点落在哪个像元，就计入哪个像元。冠层可以跨越任意多个像元。
基本用法¶
import numpy as np
import less
scene = less.Scene()
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],
grid=(40, 40), # 像元数；也可以写 resolution=1.0（m）
photons_per_m2=400, # 每平方米场景的光子数
max_bounces=None, # 记录的作用次数上限，None 表示不限
)
product = scene.simulate(sensor)
product.downwelling # （行, 列, 波段）每个像元柱的总通量，W nm⁻¹
product.upwelling
product.downwelling_irradiance() # 除以像元水平面积后的辐照度，W m⁻² nm⁻¹
product.wavelengths
product.cell_size # 像元边长 (dx, dy)，m
| 参数 | 含义 |
|---|---|
bands |
波段（nm） |
grid |
像元数（nx, ny），默认 (100, 100) |
resolution |
像元边长（m），给出时 grid 不起作用 |
photons_per_m2 |
每平方米的光子数，典型值 400–4000 |
num_photons |
总光子数，给出时代替 photons_per_m2 |
max_bounces |
只记录前若干次作用；1 表示只记录直接到达的太阳和天空辐射 |
换算为辐照度¶
downwelling 和 upwelling 是每个像元柱的总通量。换算为辐照度时需要选择面积：
surface= |
除以 | 含义 |
|---|---|---|
"horizontal"（默认） |
像元水平面积 | 水平面上的辐照度，常用于太阳辐射图 |
"slope" |
水平面积 / cos(坡度) | 贴在坡面上的辐射表读数 |
"surface" |
像元内地表的真实面积 | 按地表实际面积归一化 |
horizontal = product.downwelling_irradiance()
up = product.upwelling_irradiance("horizontal")
par = product.downwelling_integrated((400, 700), surface="horizontal") # 宽波段积分，W m⁻²
起伏地形上，结果中还有每个像元的坡度 terrain_slope（弧度）和地表面积 terrain_surface_area（m²）；平地上它们为 None。
检验：平坦朗伯地表¶
没有天空光时，平坦地表的下行辐照度等于太阳在水平面上的辐照度，上行等于它乘以反射率：
import math
horizontal_sun = np.array([1.6203, 0.7190]) # 期望的水平面辐照度，W m⁻² nm⁻¹
scene.illumination = less.Illumination(
source=less.Sun(zenith=30, azimuth=180, wavelengths=[550, 850],
irradiance=horizontal_sun / math.cos(math.radians(30))),
atmosphere=less.NoAtmosphere())
check = scene.simulate(less.IrradianceMapSensor(bands=[550, 850], resolution=1.0, photons_per_m2=400))
print("下行:", check.downwelling_irradiance().mean(axis=(0, 1))) # 约等于 horizontal_sun
print("上行:", check.upwelling_irradiance().mean(axis=(0, 1))) # 约等于 0.15 × horizontal_sun
Sun.irradiance 是垂直于太阳光方向的辐照度，所以这里除以 cos 30°。
起伏地形¶
from scipy.ndimage import gaussian_filter
rng = np.random.default_rng(42)
dem = np.zeros((128, 128))
for sigma, amplitude in [(18, 6), (9, 3), (4, 1.5), (2, 0.5)]:
dem += amplitude * gaussian_filter(rng.standard_normal((128, 128)), sigma=sigma)
dem = (dem - dem.min()) * 12.0 / (dem.max() - dem.min()) # 起伏 12 m
hills = less.Scene()
hills.size = 100.0
hills.terrain = less.Terrain(property=less.Lambertian(reflectance=0.30), dem=dem.astype(np.float32))
hills.illumination = less.Illumination(
source=less.Sun(zenith=45, azimuth=180), atmosphere=less.NoAtmosphere())
result = hills.simulate(less.IrradianceMapSensor(bands=[550], resolution=0.78, photons_per_m2=2000))
on_horizontal = result.downwelling_irradiance("horizontal")[:, :, 0]
on_slope = result.downwelling_irradiance("slope")[:, :, 0]
on_surface = result.downwelling_irradiance("surface")[:, :, 0]
ratio = result.upwelling_irradiance("horizontal")[:, :, 0] / np.maximum(on_horizontal, 1e-10)
print("上行/下行的平均值:", np.nanmean(ratio))
阳坡和阴影处的下行辐照度差别很大。坡面在像元内被重新截获的辐射会改变上行与下行的关系，因此在起伏地形和冠层中，上行/下行不能直接当作地表反照率，下行减上行也不等于局地吸收。
传感器网格可以与高程网格的分辨率不同。像元内的地表面积按实际三角形累加，同一范围内的总面积不随网格改变。坡度由“水平投影面积 / 三维面积”的反余弦得到；没有地表的像元面积为 0、坡度为 NaN。
与辐射场的区别¶
| 辐照度图 | 辐射场（第 31 章） | |
|---|---|---|
| 维度 | 二维网格 | 三维体素或每个场景元素 |
| 记录内容 | 进入和离开每个像元柱的通量 | 每个体素或元素的吸收 |
| 典型用途 | 地表辐射、地形辐射 | 冠层吸收、APAR、垂直廓线 |
配套脚本：scripts/23_irradiance_map.py