跳转至

19 - 多光谱与高光谱成像

本章用同一块玉米田模拟 RGB、Sentinel-2 和高光谱影像，并计算植被指数和红边位置。

场景

import numpy as np
import less

scene = less.Scene()
scene.size = 10.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(
    source=less.Sun(zenith=30, azimuth=150),
    atmosphere=less.SimpleSpectralAtmosphere(turbidity=2.0))

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, 0.40), jitter=0.04, margin=0.3, 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)))

RGB

rgb = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=128), bands=[650, 550, 450], quality=16))
rgb.save("maize_rgb.png")

卫星多光谱

LESS 内置 Sentinel-2A MSI 和 Landsat 8 OLI 的光谱响应函数：

bands = less.Sentinel2A()
print(bands.band_names)      # B1 … B12
print(bands.band_centers)    # 中心波长（nm）

s2 = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=128), bands=bands,
    spectral_resolution=20, quality=16))
s2.save("maize_sentinel2.tif")

每个波段的结果是在其光谱响应范围内加权积分得到的，宽波段传感器尤其需要这样计算。spectral_resolution 是积分时的光谱采样间隔（nm）：这里用 20 nm 以便快速运行，定量研究时应逐步减小到 10 nm、5 nm，直到结果不再变化。间隔越小，需要计算的波长越多。

用波段名找到红光和近红外波段，计算 NDVI：

data = s2.to_brf().data
red = data[:, :, bands.band_names.index("B4")]
nir = data[:, :, bands.band_names.index("B8")]
ndvi = (nir - red) / (nir + red)
print(f"平均 NDVI: {ndvi.mean():.3f}")

NDVI 的典型范围：裸土 0.05–0.15，稀疏植被 0.15–0.4，茂密植被 0.4–0.9。

less.Landsat8_OLI() 的用法相同。

高光谱

Hyperspectral 生成等间隔的连续波段，包含起止波长：

hyper_bands = less.Hyperspectral(start=400, stop=1000, step=10)
hyper = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=64), bands=hyper_bands, quality=16))
hyper.save("maize_hyperspectral.tif")

wavelengths = hyper_bands.get_band_centers()
mean_spectrum = hyper.to_brf().data.mean(axis=(0, 1))
print(len(wavelengths), "个波段")

波段越多，计算时间和内存越多。

红边位置

红边（约 700–750 nm）是植被光谱从红光强吸收到近红外高反射的陡升段，常用一阶导数最大处的波长表示：

edge_bands = less.Hyperspectral(start=650, stop=800, step=2)
edge = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=64), bands=edge_bands, quality=32))

edge_wavelengths = edge_bands.get_band_centers()
spectrum = edge.to_brf().data.mean(axis=(0, 1))
derivative = np.gradient(spectrum, edge_wavelengths)
print(f"红边位置: {edge_wavelengths[np.argmax(derivative)]:.0f} nm")

光照和大气

本章使用简化光谱大气，光照覆盖 300–2500 nm。模拟航空或卫星高度上的辐亮度时，使用 less.Atmosphere，并设置 observation_altitude_km 或 observation_level="toa"（见第 16 章）。

配套脚本：scripts/19_multispectral_imaging.py

下一步