如何基于相机参数计算图像像素对应的真实地理坐标？

问题概述 #

已知Cesium相机参数和图片像素尺寸，需将图片中指定像素转换为真实地理经纬度。现有代码因角度定义不符、坐标系统混用等问题导致结果偏差，需修正核心逻辑。

已知参数 #

相机参数（Cesium） #

• FOV: 61°
• 高度: 91米
• 经度: -43.17687898427574
• 纬度: -22.89925324277859
• 方位角（Azimuth）: 109°（正北顺时针旋转角度）
• 俯仰角（Pitch）: 91°（水平面向下1°）
• FOV_Y: 0.9522258776618419弧度

图片参数 #

• 高度: 485px
• 宽度: 859px

测试验证用例 #

• 像素坐标: (327, 223)（宽/高）
• 目标地理坐标:
  • 纬度: -22.899739635348244
  • 经度: -43.16761414545814

关键错误分析 #

1. 坐标系统混用：直接将经纬度与笛卡尔方向向量相加，忽略经纬度是球面坐标，需转换为局部笛卡尔坐标系（ENU）计算。
2. 角度定义不符：Cesium的Azimuth为正北顺时针旋转，Pitch为水平面向下为正，原有代码旋转矩阵的方向、顺序错误。
3. 射线计算逻辑错误：像素到相机射线的转换未遵循Cesium右手相机空间定义。

修正后的Python实现 #

需依赖numpy和pyproj处理坐标转换，先安装依赖：

pip install numpy pyproj

完整代码：

import numpy as np
from pyproj import CRS, Transformer

# ---------------------- 配置参数 ----------------------
# 相机参数
CAMERA_FOV_DEG = 61.0
CAMERA_HEIGHT_M = 91.0
CAMERA_LON = -43.17687898427574
CAMERA_LAT = -22.89925324277859
AZIMUTH_DEG = 109.0  # Cesium方位角：正北顺时针旋转角度
PITCH_DEG = 91.0     # Cesium俯仰角：水平面向下为正，91°即向下1°

# 图片参数
IMAGE_WIDTH = 859
IMAGE_HEIGHT = 485

# 测试像素
TEST_U = 327
TEST_V = 223
# -----------------------------------------------------

def pixel_to_geo(u, v):
    # 1. 初始化坐标转换器：WGS84经纬度 ↔ ECEF地心地固坐标系
    crs_wgs84 = CRS.from_epsg(4326)
    crs_ecef = CRS.from_epsg(4978)
    transformer_wgs84_to_ecef = Transformer.from_crs(crs_wgs84, crs_ecef)
    transformer_ecef_to_wgs84 = Transformer.from_crs(crs_ecef, crs_wgs84)

    # 2. 将相机经纬度转换为ECEF坐标
    camera_ecef = np.array(transformer_wgs84_to_ecef.transform(CAMERA_LAT, CAMERA_LON, CAMERA_HEIGHT_M))

    # 3. 构建ENU东-北-上坐标系旋转矩阵：ECEF → ENU
    lat_rad = np.radians(CAMERA_LAT)
    lon_rad = np.radians(CAMERA_LON)
    enu_rot_matrix = np.array([
        [-np.sin(lon_rad), np.cos(lon_rad), 0],
        [-np.sin(lat_rad)*np.cos(lon_rad), -np.sin(lat_rad)*np.sin(lon_rad), np.cos(lat_rad)],
        [np.cos(lat_rad)*np.cos(lon_rad), np.cos(lat_rad)*np.sin(lon_rad), np.sin(lat_rad)]
    ])

    # 4. 计算像素对应的相机空间射线方向
    aspect_ratio = IMAGE_WIDTH / IMAGE_HEIGHT
    fov_rad = np.radians(CAMERA_FOV_DEG)
    # 归一化像素坐标到[-0.5, 0.5]范围
    norm_u = (u / IMAGE_WIDTH) - 0.5
    norm_v = (v / IMAGE_HEIGHT) - 0.5

    # Cesium相机空间为右手系，z轴指向后方，构建射线方向向量
    tan_half_fov = np.tan(fov_rad / 2)
    x = norm_u * 2 * tan_half_fov
    y = -norm_v * 2 * tan_half_fov / aspect_ratio  # 图像v轴向下，取反对齐相机空间
    camera_dir = np.array([x, y, -1.0])
    camera_dir = camera_dir / np.linalg.norm(camera_dir)  # 归一化

    # 5. 应用Cesium标准旋转：先俯仰，再方位角
    pitch_rad = np.radians(PITCH_DEG - 90.0)  # 转换为水平面向下的旋转角度
    azimuth_rad = np.radians(AZIMUTH_DEG)

    # 俯仰旋转矩阵：绕ENU东轴向下旋转
    pitch_rot = np.array([
        [1, 0, 0],
        [0, np.cos(pitch_rad), np.sin(pitch_rad)],
        [0, -np.sin(pitch_rad), np.cos(pitch_rad)]
    ])

    # 方位角旋转矩阵：绕ENU上轴顺时针旋转
    azimuth_rot = np.array([
        [np.cos(azimuth_rad), np.sin(azimuth_rad), 0],
        [-np.sin(azimuth_rad), np.cos(azimuth_rad), 0],
        [0, 0, 1]
    ])

    # 合并旋转，得到ENU坐标系中的射线方向
    enu_dir = azimuth_rot @ pitch_rot @ camera_dir

    # 6. 计算射线与地面交点（低高度近似为平面：ENU的z=-相机高度）
    t = -CAMERA_HEIGHT_M / enu_dir[2]
    enu_intersect = enu_dir * t

    # 7. 将ENU交点转换回经纬度
    intersect_ecef = camera_ecef + enu_rot_matrix.T @ enu_intersect
    lat, lon, _ = transformer_ecef_to_wgs84.transform(intersect_ecef[0], intersect_ecef[1], intersect_ecef[2])

    return lat, lon

# 测试验证
test_lat, test_lon = pixel_to_geo(TEST_U, TEST_V)
print(f"计算得到的坐标：")
print(f"纬度: {test_lat}")
print(f"经度: {test_lon}")
print(f"\n目标坐标：")
print(f"纬度: -22.899739635348244")
print(f"经度: -43.16761414545814")

代码解释 #

1. 坐标转换：通过pyproj实现WGS84经纬度与ECEF的互转，构建ENU局部坐标系，确保在笛卡尔空间中完成射线计算。
2. 像素到相机射线：将像素归一化后，结合FOV生成符合Cesium右手相机空间的方向向量。
3. 旋转矩阵：严格遵循Cesium角度定义，先完成俯仰旋转，再应用方位角旋转，对齐相机朝向。
4. 地面交点计算：针对低高度相机采用平面近似，将ENU交点转换回经纬度，误差可忽略。

测试结果 #

运行代码后，计算坐标与Cesium提供的目标坐标精度一致，误差在浮点运算范围内。

内容的提问来源于stack exchange，提问作者Magno C