跳转至

21 - 双向反射

地表反射率随太阳和观测方向变化。植被冠层在观测方向与太阳方向一致时最亮，这时几乎看不到阴影，称为热点；稠密冠层在大观测角时通常也比星下点亮。双向反射率因子（BRF）是某一方向上地表的反射辐亮度与同样光照下理想朗伯白板反射辐亮度之比。

多角度反射率传感器

BRFSensor 一次计算一组观测方向上整个场景的平均反射率因子，不生成影像：

import less

sensor = less.BRFSensor(
    bands=[650, 850],
    directions=[(0, 0), (30, 180), (60, 0)],   # （观测天顶角, 观测方位角），度
    num_photons=500_000,
)

光子从光源出发，所有方向共用同一批光子路径，因此适合大量方向。结果的主要内容：

属性 形状 含义
.brf （方向数, 波段数） 各方向、各波段的反射率因子
.directions 列表 观测方向
.wavelengths （波段数,） 波段中心波长

玉米田的主平面

主平面是包含太阳方向和天顶的竖直平面。太阳方位角为 180° 时，方位角 180° 的方向与太阳同侧（后向散射），方位角 0° 的方向与太阳异侧（前向散射）。热点在后向散射一侧、观测天顶角等于太阳天顶角处。

import numpy as np

scene = less.Scene()
scene.size = 10.0
scene.repetitive = True       # 看作无限延伸的均匀农田
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(
    source=less.Sun(zenith=30, azimuth=180), atmosphere=less.NoAtmosphere())

maize = less.Object("maize", mesh=less.asset_path("maize.obj"))
maize.set_property(less.Prospect(cab=45, car=10, cw=0.012, cm=0.006, N=1.55))
positions = less.place.grid(scene, spacing=0.75, jitter=0.04, seed=42)
rng = np.random.default_rng(42)
scene.add(maize, positions=positions,
          scales=rng.uniform(0.9, 1.1, len(positions)),
          rotations=rng.uniform(0, 360, len(positions)))

zeniths = range(0, 65, 5)
directions = [(z, 180) for z in zeniths] + [(z, 0) for z in zeniths if z > 0]
result = scene.simulate(less.BRFSensor(bands=[650, 850], directions=directions,
                                       num_photons=500_000))

print(" 观测角     红光     近红外")
for (zenith, azimuth), (red, nir) in zip(result.directions, result.brf):
    signed = zenith if azimuth == 180 else -zenith      # 正值为后向，负值为前向
    print(f"{signed:6.0f}  {red:8.4f}  {nir:8.4f}")

result.save("maize_brf.csv")

result.show(plane="principal", style="line") 用 matplotlib 画出主平面曲线。结果也可以保存为 .json 或 .npy。

半球分布

在天顶角和方位角的网格上计算，可以画出整个半球的反射率分布：

import itertools
import matplotlib.pyplot as plt

zenith_grid = list(range(0, 75, 10))
azimuth_grid = list(range(0, 360, 30))
hemisphere = scene.simulate(less.BRFSensor(
    bands=[850], directions=list(itertools.product(zenith_grid, azimuth_grid)),
    num_photons=500_000))

values = np.zeros((len(azimuth_grid), len(zenith_grid)))
for (zenith, azimuth), brf in zip(hemisphere.directions, hemisphere.brf[:, 0]):
    values[azimuth_grid.index(azimuth), zenith_grid.index(zenith)] = brf

theta = np.deg2rad(azimuth_grid + [360])
values = np.vstack([values, values[:1]])          # 首尾相接
fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
mesh = ax.contourf(theta, zenith_grid, values.T, levels=20, cmap="viridis")
fig.colorbar(mesh, label="BRF (850 nm)")
fig.savefig("maize_brf_polar.png", dpi=150, bbox_inches="tight")

影像还是多角度反射率

需要 使用
若干方向上的平均反射率，例如双向反射曲线 BRFSensor
某个方向上的影像，关心空间纹理和阴影 OpticalImager 加 Orthographic(view_zenith=…)

BRFSensor 的方向越多，计算量也越大，但比为每个方向单独成像省得多。

有限场景与周期场景

本章的玉米田开启了周期边界，结果代表无限延伸的均匀农田。关闭周期边界时，场景外是空的，结果会包含地块边缘的影响。与观测数据或其他模型比较时，两边应使用相同的边界条件。

参数化双向反射模型

不需要三维结构时，可以直接用 RPV 模型描述地表的双向反射（见第 9 章）：

flat = less.Scene()
flat.size = 10.0
flat.terrain = less.Terrain(property=less.RPV(rho0=0.08, k=0.7, theta=-0.15, rhoc=0.05))
flat.illumination = less.Illumination(
    source=less.Sun(zenith=30, azimuth=180), atmosphere=less.NoAtmosphere())
rpv = flat.simulate(less.BRFSensor(bands=[650], directions=directions, num_photons=200_000))

配套脚本：scripts/21_brdf.py

下一步