从太阳常数到地表热辐射：Python实战遥感辐射能量计算
遥感技术中的辐射能量计算是理解地球系统能量平衡的核心技能。当我们谈论太阳常数、辐照度或黑体辐射时，这些概念不再是抽象的物理公式，而是可以通过Python代码直接量化的真实数据。本文将带您从基础物理公式出发，逐步构建完整的辐射能量计算流程，最终生成可交互的辐射曲线可视化。
1. 辐射度量基础与Python实现
在遥感能量计算中，我们首先需要明确几个关键物理量及其数学表达。太阳常数定义为地球大气层外垂直于太阳光线方向单位面积接收的辐射功率，国际公认值为1361 W/m²。这个数值看似简单，却蕴含着整个地球能量系统的基准。
让我们用Python定义这些基础物理量：
import numpy as np # 基本常数定义 SOLAR_CONSTANT = 1361 # W/m² SPEED_OF_LIGHT = 2.998e8 # m/s PLANCK_CONSTANT = 6.626e-34 # J·s BOLTZMANN_CONSTANT = 1.381e-23 # J/K SUN_EFFECTIVE_TEMP = 5778 # K # 波长范围定义 (μm转换为m) VISIBLE_RANGE = (0.38e-6, 0.76e-6) INFRARED_RANGE = (0.76e-6, 1000e-6)辐射度量体系中的核心公式是普朗克黑体辐射定律，它描述了理想黑体在特定温度下辐射能量随波长的分布：
$$ B_\lambda(T) = \frac{2hc^2}{\lambda^5} \frac{1}{e^{hc/\lambda kT} - 1} $$
对应的Python实现如下：
def planck_law(wavelength, temperature): """计算黑体辐射光谱辐射率""" wavelength = np.asarray(wavelength) exponent = PLANCK_CONSTANT * SPEED_OF_LIGHT / (wavelength * BOLTZMANN_CONSTANT * temperature) numerator = 2 * PLANCK_CONSTANT * SPEED_OF_LIGHT**2 denominator = wavelength**5 * (np.exp(exponent) - 1) return numerator / denominator注意：波长单位必须统一为米(m)，温度单位为开尔文(K)，输出结果为W·sr⁻¹·m⁻³
2. 太阳辐射模拟与可视化
太阳辐射可以近似为5778K的黑体辐射。利用前面定义的普朗克函数，我们可以生成太阳辐射的理论曲线：
import matplotlib.pyplot as plt # 生成波长数组 (0.1nm到3μm) wavelengths = np.logspace(-10, -5, 500) # 对数间隔更符合物理实际 # 计算太阳辐射曲线 solar_spectrum = planck_law(wavelengths, SUN_EFFECTIVE_TEMP) # 可视化 plt.figure(figsize=(10, 6)) plt.plot(wavelengths*1e6, solar_spectrum*1e-6, label='5778K黑体辐射') plt.title('太阳辐射光谱分布', fontsize=14) plt.xlabel('波长 (μm)', fontsize=12) plt.ylabel('光谱辐射率 (W·sr⁻¹·m⁻²·μm⁻¹)', fontsize=12) plt.xscale('log') plt.yscale('log') plt.grid(True, which='both', linestyle='--') plt.legend() plt.show()实际太阳辐射与理想黑体的差异主要体现在夫琅和费吸收线上。我们可以通过实验数据来修正理论模型：
# 加载实测太阳光谱数据 (示例) def load_solar_spectrum(): # 这里假设我们已经有一个包含实测数据的CSV文件 data = np.loadtxt('solar_spectrum.csv', delimiter=',') return data[:,0], data[:,1] # 波长(nm), 辐照度(W/m²/nm) measured_wl, measured_irradiance = load_solar_spectrum() # 将理论值与实测值对比 plt.figure(figsize=(12, 6)) plt.plot(wavelengths*1e6, solar_spectrum*1e-6, label='理论黑体辐射') plt.plot(measured_wl, measured_irradiance, 'r-', alpha=0.7, label='实测太阳光谱') plt.fill_betweenx([0, max(measured_irradiance)], VISIBLE_RANGE[0]*1e6, VISIBLE_RANGE[1]*1e6, color='yellow', alpha=0.2, label='可见光波段') plt.title('太阳辐射理论模型与实测对比', fontsize=14) plt.xlabel('波长 (nm)', fontsize=12) plt.ylabel('光谱辐照度 (W·m⁻²·nm⁻¹)', fontsize=12) plt.xlim(200, 2500) plt.legend() plt.show()3. 地球辐射能量平衡计算
地球接收的太阳辐射与自身发射的热辐射构成了能量平衡系统。我们可以计算地球接收的总太阳能：
# 地球参数 EARTH_RADIUS = 6371e3 # m EARTH_AREA = 4 * np.pi * EARTH_RADIUS**2 # 地球接收的总太阳辐射功率 total_solar_power = SOLAR_CONSTANT * np.pi * EARTH_RADIUS**2 print(f"地球接收的总太阳辐射功率: {total_solar_power/1e15:.2f} PW (拍瓦)")地球自身也向外辐射能量，根据斯特藩-玻尔兹曼定律：
$$ P = \sigma T^4 A $$
Python实现如下：
STEFAN_BOLTZMANN = 5.670374419e-8 # W·m⁻²·K⁻⁴ def earth_emission(temperature): """计算地球热辐射总功率""" return STEFAN_BOLTZMANN * temperature**4 * EARTH_AREA # 假设地球平均辐射温度 EARTH_AVG_TEMP = 255 # K (大气层顶有效温度) earth_radiated_power = earth_emission(EARTH_AVG_TEMP) print(f"地球辐射总功率: {earth_radiated_power/1e15:.2f} PW")我们可以创建一个函数来计算不同波段的辐射贡献：
def band_contribution(wavelength_range, temperature): """计算特定波长范围内的辐射能量占比""" wl = np.linspace(wavelength_range[0], wavelength_range[1], 1000) spectrum = planck_law(wl, temperature) band_power = np.trapz(spectrum, wl) total_power = STEFAN_BOLTZMANN * temperature**4 / np.pi return band_power / total_power # 计算各波段能量占比 bands = { "紫外": (0.01e-6, 0.38e-6), "可见光": (0.38e-6, 0.76e-6), "近红外": (0.76e-6, 3e-6), "中红外": (3e-6, 6e-6), "远红外": (6e-6, 15e-6), "微波": (1e-3, 1) } band_ratios = {name: band_contribution(range_, SUN_EFFECTIVE_TEMP) for name, range_ in bands.items()} # 显示结果 for name, ratio in band_ratios.items(): print(f"{name}波段能量占比: {ratio*100:.2f}%")4. 地表热辐射的实用计算模型
地表热辐射计算需要考虑大气透过率和地表比辐射率。我们可以构建一个简化的地表辐射模型：
def surface_radiation(wavelength, surface_temp, emissivity=0.95): """ 计算地表热辐射 参数: wavelength: 波长(m) surface_temp: 地表温度(K) emissivity: 地表比辐射率(0-1) """ blackbody = planck_law(wavelength, surface_temp) return emissivity * blackbody # 示例: 计算不同地表温度的热红外辐射(8-14μm) thermal_bands = np.linspace(8e-6, 14e-6, 100) temps = [280, 300, 320] # 不同地表温度(K) plt.figure(figsize=(10, 6)) for temp in temps: radiance = surface_radiation(thermal_bands, temp) plt.plot(thermal_bands*1e6, radiance*1e-6, label=f'{temp}K') plt.title('不同温度地表热红外辐射(8-14μm)', fontsize=14) plt.xlabel('波长 (μm)', fontsize=12) plt.ylabel('辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹)', fontsize=12) plt.grid(True) plt.legend(title='地表温度') plt.show()对于遥感应用，我们经常需要计算传感器接收到的表观辐射亮度，这需要考虑大气影响：
def apparent_radiance(surface_radiance, atmospheric_transmittance, path_radiance): """ 计算传感器接收的表观辐射亮度 参数: surface_radiance: 地表辐射亮度 atmospheric_transmittance: 大气透过率(0-1) path_radiance: 大气路径辐射 """ return surface_radiance * atmospheric_transmittance + path_radiance # 大气参数模型 (简化示例) def atmospheric_model(wavelength, altitude=0): """简化的大气透过率和路径辐射模型""" transmittance = 0.8 - 0.1 * (wavelength*1e6 - 10) # 经验公式 path_radiance = 0.2 * planck_law(wavelength, 280) # 假设大气辐射相当于280K黑体 return np.clip(transmittance, 0, 1), path_radiance # 计算表观辐射 wavelengths = np.linspace(8e-6, 14e-6, 100) surface_temp = 300 # K surface_emis = 0.95 surface_rad = surface_radiation(wavelengths, surface_temp, surface_emis) transmittance, path_rad = atmospheric_model(wavelengths) apparent_rad = apparent_radiance(surface_rad, transmittance, path_rad) # 可视化对比 plt.figure(figsize=(12, 6)) plt.plot(wavelengths*1e6, surface_rad*1e-6, label='地表辐射') plt.plot(wavelengths*1e6, apparent_rad*1e-6, 'r--', label='表观辐射(传感器接收)') plt.plot(wavelengths*1e6, path_rad*1e-6, 'g:', label='大气路径辐射') plt.title('地表辐射与表观辐射对比', fontsize=14) plt.xlabel('波长 (μm)', fontsize=12) plt.ylabel('辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹)', fontsize=12) plt.legend() plt.grid(True) plt.show()5. 实际案例：城市热岛效应分析
利用上述模型，我们可以模拟城市热岛效应的辐射特征。假设城区温度比郊区高3K：
# 生成模拟场景 x = np.linspace(0, 10, 100) # 10km横断面 y_rural = 290 + 2 * np.sin(x) # 郊区温度波动 y_urban = y_rural + 3 # 城区温度高3K # 计算10μm波段的辐射亮度 wavelength = 10e-6 rad_rural = surface_radiation(wavelength, y_rural) rad_urban = surface_radiation(wavelength, y_urban) # 可视化 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) ax1.plot(x, y_rural, 'g-', label='郊区温度') ax1.plot(x, y_urban, 'r-', label='城区温度') ax1.set_ylabel('温度 (K)') ax1.legend() ax1.grid(True) ax2.plot(x, rad_rural*1e6, 'g--', label='郊区辐射') ax2.plot(x, rad_urban*1e6, 'r--', label='城区辐射') ax2.set_xlabel('距离 (km)') ax2.set_ylabel(f'{wavelength*1e6:.1f}μm辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹)') ax2.legend() ax2.grid(True) plt.suptitle('城市热岛效应的温度与辐射特征模拟', fontsize=14) plt.tight_layout() plt.show()对于定量分析，我们可以计算热岛强度指数：
# 计算热岛强度 urban_avg_temp = np.mean(y_urban) rural_avg_temp = np.mean(y_rural) heat_island_intensity = urban_avg_temp - rural_avg_temp # 计算辐射差异 urban_avg_rad = np.mean(rad_urban) rural_avg_rad = np.mean(rad_rural) rad_difference = urban_avg_rad - rural_avg_rad print(f"热岛强度: {heat_island_intensity:.2f} K") print(f"10μm波段辐射亮度差异: {rad_difference*1e6:.2f} W·sr⁻¹·m⁻²·μm⁻¹")6. 高级应用：多光谱辐射能量分析
在实际遥感中，我们经常需要分析多个波段的辐射能量。下面展示如何计算Landsat 8热红外波段的辐射：
# Landsat 8 TIRS波段中心波长 landsat_bands = { 'Band 10': 10.9e-6, 'Band 11': 12.0e-6 } # 计算各波段辐射 surface_temp = 305 # K radiances = {band: surface_radiation(wl, surface_temp) for band, wl in landsat_bands.items()} # 显示结果 print("Landsat 8 TIRS波段地表辐射亮度:") for band, rad in radiances.items(): print(f"{band}: {rad*1e6:.4f} W·sr⁻¹·m⁻²·μm⁻¹") # 计算亮度温度 def brightness_temperature(radiance, wavelength): """从辐射亮度反演亮度温度""" numerator = PLANCK_CONSTANT * SPEED_OF_LIGHT / (wavelength * BOLTZMANN_CONSTANT) denominator = np.log((2 * PLANCK_CONSTANT * SPEED_OF_LIGHT**2) / (radiance * wavelength**5) + 1) return numerator / denominator # 反演亮度温度 b_temp = {band: brightness_temperature(rad, landsat_bands[band]) for band, rad in radiances.items()} print("\n反演亮度温度:") for band, temp in b_temp.items(): print(f"{band}: {temp:.2f} K")我们可以进一步分析不同地类的辐射特征：
# 定义典型地类参数 land_classes = { '水体': {'temp': 290, 'emissivity': 0.98}, '植被': {'temp': 295, 'emissivity': 0.96}, '裸土': {'temp': 305, 'emissivity': 0.94}, '建筑': {'temp': 310, 'emissivity': 0.92} } # 计算各波段辐射 band_results = {} for band, wl in landsat_bands.items(): band_results[band] = [] for cls, params in land_classes.items(): rad = surface_radiation(wl, params['temp'], params['emissivity']) band_results[band].append(rad*1e6) # 可视化 x = np.arange(len(land_classes)) width = 0.35 fig, ax = plt.subplots(figsize=(10, 6)) rects1 = ax.bar(x - width/2, band_results['Band 10'], width, label='Band 10 (10.9μm)') rects2 = ax.bar(x + width/2, band_results['Band 11'], width, label='Band 11 (12.0μm)') ax.set_ylabel('辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹)') ax.set_title('不同地类在Landsat TIRS波段的辐射特征') ax.set_xticks(x) ax.set_xticklabels(land_classes.keys()) ax.legend() plt.tight_layout() plt.show()
