DJI无人机影像地理坐标系校正

介绍

代码简介

代码的主要功能是将无人机拍摄的图像中的每个像素位置转换为地理坐标（经纬度），并计算图像的实际物理尺寸。具体包括以下几个步骤：

1. 坐标转换：

   • 从图像的像素坐标开始，转换为归一化的图像平面坐标。
   • 将这些坐标转换为相机坐标系中的坐标，考虑到焦距的影响。

2. 坐标系转换：

   • 利用无人机的姿态（俯仰角、滚转角、偏航角）将相机坐标系的坐标转换到无人机的世界坐标系。

3. 地理坐标计算：

   • 将世界坐标系中的坐标转换为地理坐标（经纬度）。这里采用简单的地球模型来近似计算地理位置。

4. 物理尺寸计算：

   • 基于图像的分辨率、传感器尺寸和焦距，计算图像的实际物理宽度和高度（以米为单位）。

代码函数

• pixel_to_latlon_matrix(image_resolution, camera_params, drone_params)：
  主要函数，接受图像分辨率、相机参数和无人机参数，返回图像中每个像素的经纬度矩阵以及图像的实际物理宽度和高度。

• rotation_matrix(pitch, roll, yaw)：
  计算给定俯仰角、滚转角和偏航角的旋转矩阵，用于将相机坐标系转换到无人机的世界坐标系。

• world_to_latlon_batch(dX, dY, drone_lat, drone_lon)：
  将世界坐标系的偏移量转换为经纬度的偏移量，基于地球的半径和无人机的当前位置计算。

示例使用

• 定义了一个示例的相机参数和无人机参数，包括焦距、传感器尺寸、图像分辨率、无人机的经纬度、海拔高度以及姿态角度。
• 通过调用 pixel_to_latlon_matrix 函数计算图像的经纬度矩阵及其物理尺寸，并打印结果。

这段代码适用于将无人机拍摄的高分辨率图像中的每个像素位置准确地映射到地理坐标系中，并提供图像的实际物理尺寸，有助于进行地理信息分析和处理。

代码

import numpy as np

def pixel_to_latlon_matrix(image_resolution, camera_params, drone_params):
    # 创建坐标网格
    x_coords, y_coords = np.meshgrid(np.arange(image_resolution[0]), np.arange(image_resolution[1]))
    
    # 1. 像素坐标转换为归一化图像平面坐标
    x_pixel = (x_coords - image_resolution[0] / 2) * camera_params['sensor_size'][0] / image_resolution[0]
    y_pixel = (y_coords - image_resolution[1] / 2) * camera_params['sensor_size'][1] / image_resolution[1]
    
    # 2. 归一化图像平面坐标转换为相机坐标系
    Z_c = camera_params['focal_length'] / 1000  # 转换为米
    X_c = x_pixel * Z_c / camera_params['focal_length']
    Y_c = y_pixel * Z_c / camera_params['focal_length']
    
    # 3. 相机坐标系转换到无人机的世界坐标系
    pitch, roll, yaw = np.radians(drone_params['angles'])
    R = rotation_matrix(pitch, roll, yaw)
    world_coords = np.dot(R, np.array([X_c.flatten(), Y_c.flatten(), Z_c * np.ones_like(X_c.flatten())]))
    
    # 4. 计算地理坐标系（经纬度）
    drone_lat = drone_params['latitude']
    drone_lon = drone_params['longitude']
    
    latlon_matrix = np.zeros((image_resolution[1], image_resolution[0], 2))
    
    # 利用批量计算的世界坐标
    dX = world_coords[0, :].reshape(image_resolution[1], image_resolution[0])
    dY = world_coords[1, :].reshape(image_resolution[1], image_resolution[0])
    
    lat, lon = world_to_latlon_batch(dX, dY, drone_lat, drone_lon)
    
    latlon_matrix[:, :, 0] = lat
    latlon_matrix[:, :, 1] = lon
    
    # 计算图像的实际物理尺寸
    image_width_meters = camera_params['sensor_size'][0] * (image_resolution[0] / camera_params['focal_length']) * 1000
    image_height_meters = camera_params['sensor_size'][1] * (image_resolution[1] / camera_params['focal_length']) * 1000
    
    return latlon_matrix, image_width_meters, image_height_meters

def rotation_matrix(pitch, roll, yaw):
    R_x = np.array([
        [1, 0, 0],
        [0, np.cos(roll), -np.sin(roll)],
        [0, np.sin(roll), np.cos(roll)]
    ])
    
    R_y = np.array([
        [np.cos(pitch), 0, np.sin(pitch)],
        [0, 1, 0],
        [-np.sin(pitch), 0, np.cos(pitch)]
    ])
    
    R_z = np.array([
        [np.cos(yaw), -np.sin(yaw), 0],
        [np.sin(yaw), np.cos(yaw), 0],
        [0, 0, 1]
    ])
    
    return np.dot(R_z, np.dot(R_y, R_x))

def world_to_latlon_batch(dX, dY, drone_lat, drone_lon):
    a = 6378137.0  # Semi-major axis in meters
    f = 1 / 298.257223563  # Flattening
    b = a * (1 - f)  # Semi-minor axis
    
    lat_rad = np.radians(drone_lat)
    
    N = a / np.sqrt(1 - (2*f - f**2) * np.sin(lat_rad)**2)
    
    m_per_deg_lat = (np.pi / 180) * (a * (1 - f)) / (1 - (2*f - f**2) * np.sin(lat_rad)**2)**1.5
    m_per_deg_lon = (np.pi / 180) * (N * np.cos(lat_rad))
    
    delta_lat = dY / m_per_deg_lat
    delta_lon = dX / m_per_deg_lon
    
    object_lat = drone_lat + delta_lat
    object_lon = drone_lon + delta_lon
    
    return object_lat, object_lon

# 示例使用
camera_params = {
    'focal_length': 24,  # in mm
    'sensor_size': (36.0, 24.0),  # in mm, (width, height)
    'image_resolution': (8064, 6048)  # in pixels, (width, height)
}

drone_params = {
    'latitude': 40.052235,  # degrees
    'longitude': 116.243683,  # degrees
    'altitude': 179,  # meters
    'angles': (10, 5, 30)  # degrees, (pitch, roll, yaw)
}

latlon_matrix, image_width_meters, image_height_meters = pixel_to_latlon_matrix(camera_params['image_resolution'], camera_params, drone_params)

print(f'LatLon Matrix Shape: {latlon_matrix.shape}')
print(f'Image Width: {image_width_meters:.2f} meters')
print(f'Image Height: {image_height_meters:.2f} meters')

打印结果

LatLon Matrix Shape: (6048, 8064, 2)
Image Width: 41.47 meters
Image Height: 30.24 meters