从太阳常数到地表热辐射：手把手教你用Python计算遥感中的辐射能量

从太阳常数到地表热辐射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, label5778K黑体辐射) plt.title(太阳辐射光谱分布, fontsize14) plt.xlabel(波长 (μm), fontsize12) plt.ylabel(光谱辐射率 (W·sr⁻¹·m⁻²·μm⁻¹), fontsize12) plt.xscale(log) plt.yscale(log) plt.grid(True, whichboth, 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-, alpha0.7, label实测太阳光谱) plt.fill_betweenx([0, max(measured_irradiance)], VISIBLE_RANGE[0]*1e6, VISIBLE_RANGE[1]*1e6, coloryellow, alpha0.2, label可见光波段) plt.title(太阳辐射理论模型与实测对比, fontsize14) plt.xlabel(波长 (nm), fontsize12) plt.ylabel(光谱辐照度 (W·m⁻²·nm⁻¹), fontsize12) 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, emissivity0.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, labelf{temp}K) plt.title(不同温度地表热红外辐射(8-14μm), fontsize14) plt.xlabel(波长 (μm), fontsize12) plt.ylabel(辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹), fontsize12) 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, altitude0): 简化的大气透过率和路径辐射模型 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(地表辐射与表观辐射对比, fontsize14) plt.xlabel(波长 (μm), fontsize12) plt.ylabel(辐射亮度 (W·sr⁻¹·m⁻²·μm⁻¹), fontsize12) 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(城市热岛效应的温度与辐射特征模拟, fontsize14) 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(f10μ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, labelBand 10 (10.9μm)) rects2 ax.bar(x width/2, band_results[Band 11], width, labelBand 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()