【图像处理基石】遥感图像高度信息提取：Python实战全流程+常用库汇总 - 教程

在地形测绘、灾害监测、生态评估等场景中，从遥感图像提取高度信息（高程数据）是核心需求。本文将从技术原理出发，用Python实现完整的高度提取流程，并汇总遥感处理中常用的工具库，帮你快速上手遥感数据处理。

一、核心实战：遥感图像高度提取（基于立体像对）

遥感图像提取高度的核心是立体像对（Stereo Pair）+ 视差计算——通过同一区域不同视角的两张图像，计算地物的像素位置差异（视差），再结合摄影测量原理推导高程。

1.1 技术原理速览

• 立体像对：同一区域的左、右两张遥感图（重叠度需>60%）；
• 视差（Disparity）：同一地物在左右图中的像素偏移，距离越近，视差越大；
• 高程公式：已知相机参数（焦距f、基线距B，从图像元数据获取），高程计算为：
  $h=\frac{B\cdot f}{d}+h_0$
  其中，$d$为视差，$h_0$为参考高度（如海平面）。

1.2 环境准备

需安装基础依赖库，用于图像处理、遥感数据读写：

pip install opencv-python numpy gdal matplotlib rasterio  # 核心库

1.3 完整代码实现（分步骤）

步骤1：读取遥感立体像对与地理信息

遥感图像常用TIFF格式，需保留地理坐标（用于后续结果映射），这里用gdal读取图像数据和地理变换参数：

import cv2
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt
def read_rs_image(path):
"""读取遥感图像，返回灰度数组+地理变换参数"""
ds = gdal.Open(path)  # 打开TIFF文件
if ds is None:
raise ValueError(f"无法读取图像：{path}")
# 读取图像数组（多通道转灰度）
arr = ds.ReadAsArray()  #  shape: (通道数, 高度, 宽度)
if len(arr.shape) == 3:  # 多通道（如RGB）转灰度
arr = cv2.cvtColor(arr.transpose(1,2,0), cv2.COLOR_BGR2GRAY)
else:
arr = arr.transpose(1,2)  # 单通道调整为 (高度, 宽度)
# 获取地理变换参数（用于后续保存高程图的坐标映射）
geo_transform = ds.GetGeoTransform()  # (左上角x, x分辨率, 0, 左上角y, 0, -y分辨率)
return arr, geo_transform
# 读取左、右立体像对（需替换为你的文件路径）
left_img, left_geo = read_rs_image("left.tif")
right_img, right_geo = read_rs_image("right.tif")
# 查看图像尺寸（需确保左右图尺寸一致）
print(f"图像尺寸：左图{left_img.shape}，右图{right_img.shape}")

步骤2：立体匹配计算视差图

用OpenCV的SGBM（Semi-Global Block Matching） 算法计算视差——比基础块匹配更鲁棒，适合遥感图像的复杂地形：

def compute_disparity_map(left, right):
"""SGBM算法计算视差图"""
# SGBM参数（需根据图像分辨率调整，关键参数已标注）
window_size = 3  # 匹配窗口大小（奇数，分辨率高则调小）
min_disp = 0     # 最小视差值
num_disp = 128   # 视差范围（必须是16的倍数，范围越大越慢）
# 初始化SGBM匹配器
left_matcher = cv2.StereoSGBM_create(
minDisparity=min_disp,
numDisparities=num_disp,
blockSize=window_size,
P1=8 * 3 * window_size**2,  # 平滑项：相邻像素视差差异惩罚（小值）
P2=32 * 3 * window_size**2, # 平滑项：远距离像素视差差异惩罚（大值）
disp12MaxDiff=1,            # 左右视差一致性检查阈值
uniquenessRatio=15,         # 视差唯一性阈值（过滤模糊匹配）
speckleWindowSize=100,      # 噪声过滤窗口大小
speckleRange=32             # 噪声过滤视差范围
)
# 计算视差（结果需除以16转为实际视差值）
disparity = left_matcher.compute(left, right).astype(np.float32) / 16.0
return disparity
# 计算视差图
disparity_map = compute_disparity_map(left_img, right_img)
# 可视化视差图（用jet色表，值越大代表视差越大）
plt.figure(figsize=(8, 6))
plt.imshow(disparity_map, cmap="jet")
plt.colorbar(label="视差值（像素）")
plt.title("SGBM立体匹配视差图")
plt.axis("off")
plt.show()

步骤3：视差图转高程图

根据摄影测量公式计算高程，需过滤无效视差（如负值、0，避免除以零错误）：

def disparity_to_elevation(disparity, focal_len, baseline, h0=0):
"""视差图转换为高程图"""
# 1. 过滤无效视差（视差需为正，排除遮挡、噪声区域）
valid_mask = disparity > 1e-6  # 避免除以零
elevation = np.zeros_like(disparity, dtype=np.float32)
# 2. 高程公式：h = (基线距*焦距)/视差 + 参考高度
elevation[valid_mask] = (baseline * focal_len) / disparity[valid_mask] + h0
elevation[~valid_mask] = np.nan  # 无效区域标记为NaN（后续可插值填充）
return elevation
# 相机参数（从遥感图像元数据获取，此处为示例值，需替换为实际值！）
focal_length = 1500  # 焦距（像素单位，从相机标定或图像元数据获取）
baseline = 50        # 基线距（米，两相机中心的距离）
h0 = 0               # 参考高度（如海平面，米）
# 计算高程图
elevation_map = disparity_to_elevation(disparity_map, focal_length, baseline, h0)

步骤4：保存高程图为GeoTIFF（保留地理坐标）

将高程结果保存为标准GeoTIFF，方便后续用ArcGIS、QGIS等工具分析：

def save_elevation_tif(elevation, geo_transform, output_path):
"""保存高程图为GeoTIFF（保留地理坐标）"""
rows, cols = elevation.shape
# 1. 创建GeoTIFF驱动
driver = gdal.GetDriverByName("GTiff")
# 2. 创建输出文件（宽度、高度、波段数、数据类型）
out_ds = driver.Create(
output_path, cols, rows, 1,
gdal.GDT_Float32  # 浮点型，支持NaN
)
# 3. 设置地理变换参数（继承左图坐标）
out_ds.SetGeoTransform(geo_transform)
# 4. 写入高程数据
out_ds.GetRasterBand(1).WriteArray(elevation)
# 5. 标记NoData值（NaN对应GDAL的NoData）
out_ds.GetRasterBand(1).SetNoDataValue(np.nan)
# 6. 释放资源
out_ds.FlushCache()
print(f"高程图已保存至：{output_path}")
# 保存高程图
save_elevation_tif(elevation_map, left_geo, "elevation.tif")
# 可视化高程图（用terrain色表，模拟地形起伏）
plt.figure(figsize=(8, 6))
plt.imshow(elevation_map, cmap="terrain")
plt.colorbar(label="高程（米）")
plt.title("遥感图像提取的高程图")
plt.axis("off")
plt.show()

1.4 关键注意事项（避坑指南）

1. 数据质量是前提：
   • 立体像对应先做辐射校正（消除光照差异）和几何校正（确保同名点对齐）；
   • 避免云雾、阴影区域（会导致视差计算错误）。
2. 参数优化技巧：
   • 图像分辨率高 → 减小window_size（如3），分辨率低 → 增大（如5）；
   • 若视差图有大量空洞 → 增大num_disp或减小uniquenessRatio。
3. 精度提升方法：
   • 用地面控制点（GCP）校准相机参数（焦距、基线）；
   • 对无效区域（NaN）用反距离加权（IDW）插值填充；
   • 结合LiDAR数据或已有DEM验证精度。

二、扩展：遥感图像处理常用库汇总

除了上文的opencv-python、gdal，还有很多专用库可应对不同场景，按功能分类整理如下：

2.1 数据读写与基础处理

库名	核心功能	优势与适用场景	示例代码片段
rasterio	遥感影像读写、坐标处理	API更Python化，支持上下文管理器，替代GDAL	with rasterio.open("img.tif") as src: arr = src.read()
geopandas	矢量数据（Shapefile/GeoJSON）处理	与Pandas兼容，可裁剪影像至矢量边界	gdf = gpd.read_file("roi.shp"); rasterio.mask.mask(src, gdf.geometry)
xarray	多维遥感数据（时序/多光谱）处理	支持标签索引，适合MODIS/Landsat时序分析	ds = xarray.open_rasterio("time_series.tif"); ds.sel(band=1)

2.2 光谱与高光谱处理

库名	核心功能	优势与适用场景	示例代码片段
spectral	高光谱数据降维、端元分析、矿物识别	专为高光谱设计，支持ENVI格式	img = spectral.open_image("hyper.hdr"); pca = spectral.principal_components(img)
Py6S	大气校正（消除大气散射/吸收影响）	基于6S辐射传输模型，提升定量精度	from Py6S import SixS; s = SixS(); s.atmos_profile = SixS.AtmosProfile.PredefinedType("MidlatitudeSummer")
pyhdf	HDF格式数据（MODIS/OMI卫星）读写	支持NASA卫星Level-1/2数据	from pyhdf.SD import SD; hdf = SD("modis.hdf"); data = hdf.select("reflectance")

2.3 地理投影与坐标转换

库名	核心功能	优势与适用场景	示例代码片段
pyproj	坐标系转换（WGS84→UTM、高斯投影等）	基于PROJ库，支持所有主流坐标系	transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650"); x,y = transformer.transform(30, 120)
cartopy	地理数据可视化（带投影的地图绘制）	替代basemap，支持UTM/墨卡托等投影	ax = plt.axes(projection=ccrs.UTM(50)); ax.imshow(arr, transform=ccrs.UTM(50))

2.4 深度学习与语义分割

库名	核心功能	优势与适用场景	备注
torchgeo	遥感专用深度学习库（数据集+模型）	PyTorch生态，支持EuroSAT/LoveDA数据集	适合快速搭建土地覆盖分类、变化检测模型
segmentation-models-pytorch	语义分割模型（U-Net/DeepLab）	支持多光谱输入，可加载预训练权重	需调整输入通道数（如4通道 Sentinel-2）
tensorflow-graphics	3D视觉工具（立体匹配、相机几何优化）	TensorFlow生态，适合提升视差计算精度	需一定3D数学基础

2.5 专业遥感工具包

库名	核心功能	优势与适用场景	备注
Orfeo Toolbox (OTB)	全流程遥感处理（辐射校正、分类、变化检测）	开源替代ENVI/ERDAS，支持Python绑定	需要单独安装OTB，再导入otbApplication
scikit-image	图像处理（滤波、分割、特征提取）	API清晰，与NumPy无缝兼容	适合遥感图像预处理（去噪、边缘检测）

三、结语

本文从实战出发，覆盖了遥感图像高度提取的全流程，并汇总了不同场景的专用库。实际应用中，需根据数据类型（如高光谱、时序遥感）和需求（如精度、速度）选择合适的工具组合。