大疆M3M/P4M 航拍图像辐射定标流程及python实现

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

Github：https://github.com/xhdidi/TIFF/blob/master/Radiometric_calibrations.py

前言

为什么要对航拍数据进行辐射定标

单波段相机获取的是DN值，即像元亮度值，其大小与传感器的辐射分辨率、地物发射频率、大气透过率、散射率有关。无人机载相机作为农业巡航设备，其值大小受太阳角度、双向反射分布函数效应、云的阴影、相机增益及曝光时间等因素影响，不能直接反映作物生长状况，因此需要转化为反射率，即物体表面反射的辐射量/接收的辐射量，才能保证所建立的检测模型在不同设备、不同时间和不同天气条件下的鲁棒性。

---

一、辐射定标的定义

传感器辐射定标分相对辐射定标和绝对辐射定标两种，前者标定和纠正传感器不同探元间的相对响应差异，后者则是利用绝对参照体建立图像DN值（或灰度信息）与实际地物辐亮度的响应关系，将图像DN值转换成辐亮度、反射率等绝对地物信息（具有具体意义，如反射率信息）。

二、DJI M3M/P4M的辐射定标原理（反射率的计算）

图像的像素值代表阳光（入射光）照在地物（target）上的反射值，同时无人机顶端配有光强传感器以获取入射光的信号值。有研究表明，可见光图像的辐射定标方程是非线性的，而多光谱图像的辐射定标方程是线性的[2]。因此，某波段的反射率公式可定义如下，其中${\rho }_X$是调节图像信号与多光谱光强传感器信号之间相互转化的参数，使入射光（多光谱光强传感器信号值）与反射光（图像信号值）的量级（单位）保持统一。
$$X_{ref}=\frac{X_{reflected}}{X_{incident}}=\frac{X_{camera}}{X_{LS}}\times {\rho }_x$$
由于光强传感器在不同波段存在感度差异，即在相同光源的条件下不同波段的相机和不同波段的多光谱光强传感器会得到不同的信号值，因此需要对相机和感光器各自做感度校准。
在此以NIR波段为基准，所有其他波段的相机都参照 NIR 波段的相机的感度做校准，得到校准参数$pCa{m}_x$和$pL{S}_x$，各波段与NIR的调节参数定义如下：
$${\rho }_x={\rho }_{NIR}\times \frac{pCa{m}_x}{pL{S}_x}$$
则各波段的反射率定义可写作：
$$X_{ref}=\frac{X_{reflected}}{X_{incident}}=\frac{X_{camera}}{X_{LS}}\times {\rho }_x=\frac{X_{camera}\cdot pCa{m}_x}{X_{LS}\cdot pL{S}_x}\times {\rho }_{NIR}$$
$pCa{m}_x$在元数据中保存为：Xmp.drone-dji.SensorGainAdjustment
$X_{LS}\cdot pL{S}_x$在元数据中保存为：Xmp.drone-dji.Irradiance
${\rho }_{NIR}$没有给出具体数值，需要参照植被指数计算软件如ENVI等进行手动线性校正，考虑到植被指数的数值本身并无意义，只是作为作物生长状态的特征表达，故此步骤可忽略。
下面对$X_{camera}$进行计算（图像校正）。

1. 相机位置误差校正

多个单波段相机拍摄的物理位置存在微小差异，为了能够将不同相机拍摄的图像相匹配（以计算植被指数），需要进行位置误差校正。

在此以NIR相机为基准，TIF文件的元数据保存有该波段相机相对于NIR波段相机的位置偏移量（像素），分别为 Xmp.drone-dji.RelativeOpticalCenterX ， Xmp.drone-dji.RelativeOpticalCenterY，直接对原始图像进行padding后裁剪即可。

2. 相机曝光时间误差校正

不同相机的曝光时间存在差异，为了减轻曝光时间差异导致的位置偏移影响，首先对图像进行平滑处理（高斯滤波等），再对图像进行配准（sobel滤波器进行特征提取，Enhanced Correlation Coefficient Maximization做对齐）。

3. 黑电平校正

用归一化的像素值减去黑电平值即可，黑电平值在元数据里保存为Xmp.Camera.BlackCurrent，归一化前值为4096。

4. 暗角补偿

暗角是指画面四角有变暗的现象，从照片来看暗角其实就是指拍摄亮度均匀的场景时，画面四角却出现扇形向外延伸的渐暗区域。

暗角补偿公式如下，其中$x,y$为像素坐标，$r$为当前像素到暗角补偿中心的距离，$k$为暗角补偿的多项式系数：
$$I(x,y)=I_0(x,y)\times (k_5\cdot r^6+k_4\cdot r^5+k_3\cdot r^4+k_2\cdot r^3+k_1\cdot r^2+k_0\cdot r+1.0)).$$
$$r=\sqrt{(x-centerX{)}^2+(y-centerY{)}^2}.$$
暗角补偿中心坐标在元数据中保存为：Xmp.drone-dji.CalibratedOpticalCenterX, Xmp.drone-dji.CalibratedOpticalCenterY
暗角补偿系数在元数据中保存为：Xmp.drone-dji.VignettingData

5. 畸变校正

相机的径向畸变和切向畸变模型定义如下：
$$\{\begin{array}{rl}{u}^{\prime }=u(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)& \\ v^{\prime }=v(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)& \end{array}$$
$$\{\begin{array}{rl}{u}^{\prime }=u+[2p_{1}v+p_2(r^2+2u^2)]& \\ v^{\prime }=v+[p_1(r^2+2v^2)+2p_{2}u]& \end{array}$$

像素坐标系与相机坐标系之间的转换由内参决定（相机参数），相机坐标系到世界坐标系的转换由外参决定（相机位姿）。因此畸变校准只需要计算相机内参。
相机内参定义如下：
$$\left(\begin{array}{ccc}{f}_x& 0& centerX+cx\\ 0& f_y& centerY+cy\\ 0& 0& 1\end{array}\right)\text{\hspace{1em}(2)}$$

$fx,fy,cx,cy,k1,k2,p1,p2,p3$在元数据中保存为Xmp.drone-dji.DewarpData，其中$fx,fy,cx,cy$为相机参数，$k1,k2,p1,p2,p3$为畸变校准的多项式系数，分别用于校正径向畸变和切向畸变。
采用cv2.undistort函数实现，校正原理在此不详细阐述。

Summary

对以上畸变进行校正，并考虑曝光增益、曝光时间，得到图像信号值定义如下：
$$X_{camera}=(I_x-I_{blackcurrent})\times Corretion/(X_{gain}\ast \frac{X_{etime}}{1e6})$$
其中$Correction$包括暗角校正和畸变校正，
$X_{gain}$在元数据中保存为：Xmp.drone-dji.SensorGain
$X_{etime}$在元数据中保存为：Xmp.drone-dji.ExposureTime

三、代码实现

class TIF_Image():
    def __init__(self, filename):
        self.filename = filename
        self.pyexiv2 = Image(self.filename)
        self.xmp = self.pyexiv2.read_xmp()
        self.exif = self.pyexiv2.read_exif()

        self.W, self.L = int(self.exif['Exif.Image.ImageWidth']), int(self.exif['Exif.Image.ImageLength'])

        t_raster = gdal.Open(filename)
        self.image = t_raster.ReadAsArray()  # <class 'numpy.ndarray'>

    def physical_position_calibration(self):
        '''相机物理位置误差校准'''
        pos_x = int(np.around(eval(self.xmp['Xmp.drone-dji.RelativeOpticalCenterX']), decimals=0, out=None))
        pos_y = int(np.around(eval(self.xmp['Xmp.drone-dji.RelativeOpticalCenterY']), decimals=0, out=None))
        
        # pos_x<0, pos_y>0
        image_pad = np.pad(self.image, ((0, pos_y), (-pos_x, 0)), 'constant', constant_values=(0, 0))
        self.image = image_pad[pos_y:self.L+pos_y, :self.W]

    def exposure_time_calibration(self):
        '''相机曝光时间误差校准(高斯滤波)'''
        self.image = cv2.GaussianBlur(self.image, (3,3), 0)

    def cam_correction(self):
        # 暗角补偿
        center_x = eval(self.xmp['Xmp.drone-dji.CalibratedOpticalCenterX'])
        center_y = eval(self.xmp['Xmp.drone-dji.CalibratedOpticalCenterY'])
        k0, k1, k2, k3, k4, k5 = self.xmp['Xmp.drone-dji.VignettingData'].split(', ')
        k0, k1, k2, k3, k4, k5 = eval(k0), eval(k1), eval(k2), eval(k3), eval(k4), eval(k5)

        max_distance = math.sqrt(
            max((vertex[0] - center_x) ** 2 + (vertex[1] - center_y) ** 2 for vertex in
                [[0, 0], [0, W], [L, 0], [L, W]]))
        for x in range(0, L):
            for y in range(0, W):
                distance = math.sqrt(pow(x - center_x, 2) + pow(y - center_y, 2))
                r = distance / max_distance
                k = k5 * r**6 + k4 * r**5 + k3 * r**4 + k2 * r**3 + k1 * r**2 + k0 * r + 1.0
                self.image[x, y] *= k

        # 畸变校准
        fx, fy, cx, cy, k1, k2, p1, p2, p3 = self.xmp['Xmp.drone-dji.DewarpData'].split(';')[-1].split(',')

        distcoeffs = np.float32([eval(k1), eval(k2), eval(p1), eval(p2), eval(p3)])
        cam_matrix = np.zeros((3, 3))
        cam_matrix[0, 0] = eval(fx)
        cam_matrix[1, 1] = eval(fy)
        cam_matrix[0, 2] = eval(cx) + center_x
        cam_matrix[1, 2] = eval(cy) + center_y
        cam_matrix[2, 2] = 1

        self.image = cv2.undistort(self.image, cam_matrix, distcoeffs)


    def cam(self):
        bit = eval(self.exif['Exif.Image.BitsPerSample'])
        self.image /= pow(2.0, bit)
        # black_current = eval(self.xmp['Xmp.Camera.BlackCurrent']) / 4096
        # self.image -= black_current

        self.cam_correction()

        pcam = eval(self.xmp['Xmp.drone-dji.SensorGainAdjustment'])
        gain = eval(self.xmp['Xmp.drone-dji.SensorGain'])
        etime = eval(self.xmp['Xmp.drone-dji.ExposureTime'])
        denominator = eval(self.xmp['Xmp.drone-dji.Irradiance'])

        self.image *= pcam/(gain*etime/1e6)/denominator

    def calibration(self):
        '''总校准函数'''
        self.physical_position_calibration()
        self.exposure_time_calibration_cv()
        self.cam()

---

参考文献

[1] P4M图像处理指南
[2] 章凌翔. 基于无人机多光谱影像油菜产量预测及倒伏风险评价[D].西北农林科技大学,2022.DOI:10.27409/d.cnki.gxbnu.2022.000191.