[DJI P4M多光谱数据处理] 辐射校正获取双向反射率
关键词：DJI 精灵4 多光谱 双向反射率 辐射校正
日期：1/7/2021 - 2/7/2021
2019年，DJI大疆创新正式发布了精灵4多光谱版无人机(P4 Multispectral，P4M)，这一多光谱无人机集成系统大大降低了多光谱遥感数据获取的准入门槛和成本水平，大幅提升了农林水草等遥感数据的获取效率，是遥感从业人员的一大利器。
我十分高兴的是，终于不用再分别购买无人机和多光谱相机进行集成了，终于有一个十分稳定的解决方案了。大疆出品，必属精品！
然而，据我了解，DJI多光谱影像的
图像处理指南文档中一些用词并不严谨和专业，甚至有一些模棱两可，而且并没有明确告知双向反射率的计算方式，仅仅介绍了植被指数的计算方式。因此我认为有必要重新系统介绍一下如何对多光谱影像进行辐射校正进而获取双向反射率。
1. DJI P4M多光谱数据介绍
1.1. 影像数据
DJI P4M集成了1 个 可见光相机(RGB合成影像)及 5 个多光谱相机(蓝光，绿光，红光，红边和近红外单波段影像)，所有相机均拥有 200 万像素解析力及配备全局快门。
Blue: 450 ± 16 nm； Green: 560 ± 16 nm； Red: 650 ± 16 nm； Red-edge: 730 ± 16 nm； Near infrared (NIR): 840 ± 26 nm； 需要注意的是：
单次拍摄，不同单波段相机的 曝光时间与增益可能不同；不同次拍摄，不同单波段相机的曝光时间与增益也往往不同。 由于多光谱相机中各个波段相机的位置存在差异，因此不同波段影像的同一位置像元并非对应相同的地物，并不能直接叠加在一起。
需要对不同波段的影像进行
配准处理，即以某一波段的影像为参考位置，其他波段修正位置偏移。由于各波段相机镜头以及内部光路一致，它们的
暗角效应是一致的。也就是说，不同波段相机采用完全相同的暗角效应
校正系数。
1.2. 元数据
令人惊讶的是，DJI提供了非常多的元数据（EXIF metadata，以及 XMP metadata），可以用于影像的几何校正、辐射校正、飞行状态解析等等。
航点位置相关 GNSS RTK：相对高程（RelativeAltitude）、绝对高程(AbsoluteAltitude)、经度(GpsLongtitude)、纬度(GpsLatitude)、RTK状态(RtkFlag)、RTK经度标准差(RtkStdLon)、RTK纬度标准差(RtkStdLat)、RTK高程标准差(RtkStdHgt)。
飞行姿态相关：X方向速度(FlightXSpeed)、Y方向速度(FlightYSpeed)、Z方向速度(FlightZSpeed)、偏航角(FlightYawDegree)、旋转角(FlightRollDegree)、俯仰角(FlightPitchDegree)。
相机姿态相关：其实对应的是云台转角，即云台俯仰角(GimbalPitchDegree)、云台偏航角(GimbalYawDegree)、云台旋转角(GimbalRollDegree)。
辐射校正相关：位数(BitsPerSample)、中心波长(CentralWavelength)、传感器索引(SensorIndex)、暗电流值(BlackCurrent)、曝光时间(ExposureTime)、感光速度(ISOSpeed)、增益(SensorGain)、增益相对校正参数(SensorGainAdjustment)、暗角校正参数(VignettingPolynomial 或 VignettingData)、暗角中心像元位置(VignettingCenter)
辐照度传感器相关：太阳辐照度传感器测量值(SunSensor)、辐照度传感器曝光时间(SunSensorExposureTime)、辐照度传感器增益(IrradianceGain)。
几何相关：影像宽度(ImageWidth)、影像高度(ImageHeight)、像主点(PrincipalPoint)、焦距(PerspectiveFocalLength)、焦距像元数(CalibratedFocalLength)、畸变校正参数(DewarpData 或 PerspectiveDistortion)、相对光心偏移X(RelativeOpticalCenterX)、相对光心偏移Y(RelativeOpticalCenterY)。
注意：在影像的metadata中，以上参数往往存在重复，比如暗角校正参数在在VignettingPolynomial 或 VignettingData中均可找到。
DJI 还提供了很多其他有用的metadata，就不一一列举了，这里有一个示例文件（提取码：1234），可供查看。
2. 逐影像双向反射率计算
我们常说要计算反射率，但是这个称谓并不准确，应该叫做
双向反射率或者双向反射因子（具体为什么，请参考我以前的博客）。有时我们仅仅需要计算每幅影像的逐像元双向反射率，而且往往这也是后处理的输入数据。
2.1. 暗电流校正
我们通常对光谱仪或高光谱相机做暗电流校正，多光谱相机往往被忽视这一问题。既然DJI明确提供了暗电流校正参数，那何不校正一下。公式很简单：
2.2. 曝光时间与增益校正
如前所说，不同次拍摄以及不同波段的曝光时间、增益等参数可能不同 （幸好，不同波段相机的光圈大小是一致的，快门速度(ISO speed)也是一致的），因此为确保它们之间可以相比，需要对DN值进行标准化。
不过，目前我还没搞懂为什么不同次拍摄时光圈一致，快门速度一致时，曝光时间存在差异。有知道的，请打在评论区。
2.3. 暗角校正
暗角效应体现为亮度由中心向边缘进行辐射状的变暗，其通常主要是由于镜头边缘的成像光线与镜头光轴有一个较大的夹角，造成边角失光。当然，还有其他一些原因，但不在本博客的讨论范围之内。
典型的暗角效应模型为
多项式模型，对辐射状暗角进行补偿：其中，是像素点到补偿中心（通常为影像中心）的像素距离：
六次多项式模型的校正系数可以由元数据中的
VignettingPolynomial或VignettingData获得。补偿中心由元数据中的VignettingCenter获得。
2.4. 成像畸变校正
相机在成像过程中，由于真实的相机成像过程并非严格遵从小孔成像模型，会发生
径向畸变和切向畸变(具体可以参考我以前的博客)。通常采用Brown校正模型进行校正。
其中，是图像物理坐标点到图像物理中心（通常为影像中心）的物理距离：，注意其与暗角校正中的坐标系的差异。
上述校正所在的坐标系为
图像物理坐标系，因此最后需要转换到图像像素坐标系中：
其中 和，为影像中心。
成像畸变的校正参数可以由元数据中的
DewarpData获得，注意其顺序为。
如果不采用 运动恢复结构 sfm算法对影像进行三维重建，则需上述公式进行单独的畸变校正；如果采用SfM算法，则不必额外进行畸变校正了，因为SfM过程中已经包含了畸变校正的过程(具体可以参考我以前的博客)。
2.5. 双向反射率计算
经过以上一系列校正，最终我们得到
标准化处理后的DN值即，综合校正公式为：
其中上标代表归一化到[0, 1]区间的像元值。为暗角校正模块。
注意：经上述公式处理后，并非得到双向反射率或者其他物理值，而仍然是DN值，只不过是标准化处理后的DN值。
通常由DN值转换为
双向反射率有两种方式：
单一参考板校正方法：即仅采用一块标准反射率参考板进行校正。 假定参考板与地物具有相同的辐照度条件，则地物双向反射率 和参考板双向反射率 ，联立该方程组，可以得到
注意：上式中假定辐亮度与DN值的关系满足而非，即接受的辐亮度为0时，DN值为0。其实通过暗电流校正，前者是成立的。
多参考板校正方法：即采用多块参考板构建经验校正模型进行校正。 由于某些多光谱传感器自身的非线性问题或者参考板老化等问题，采用多参考板进行校正往往更为合适。可以采用线性校正模型：
即采用多块参考板构建多个方程组，计算和，进而用于地物影像的辐射校正（即双向反射率的计算）
3. 逐影像植被指数计算
如果采用参考板进行辐射校正获取双向反射率影像，则逐影像植被指数的计算非常自然，仅需要将不同波段影像之间的位置偏移进行校正。
DJI规定了以NIR波段影像为参考，以此需要对其他波段的影像进行位置偏移校正。 考虑到不同波段相机的物理位置非常相近，所摄取的地物范围极为相近，因此仅仅通过平移即可。每个波段影像的元数据中的
RelativeOpticalCenterX和RelativeOpticalCenterY给定了需要平移的像元数。图像处理指南中建议飞行时有航点悬停的情况下采用此种方式，不过据我们的测试，飞行速度比较慢（< 5 m/s）时，也可以直接平移。如果确实发现平移不行，那就需要其他的配准方式了。
另外一种计算植被指数的方式是采用
太阳辐照度传感器 sun sensor，但是注意此种方式目前仅适用于计算比值型植被指数。此外，单纯采用DJI的sun sensor是无法计算出双向反射率的，因为欲计算双向反射率，需要采用如下公式，而DJI并未提供足够的参数用于计算。
由于DJI的元数据中并未提供和，因此上式无法求解。
不仅如此，某个波段的相机与其对应的sun sensor的
光谱响应函数也应当保持一直，即在相同的入射辐射通量条件下，它们的响应值是一致的。有关如何用DJI现有的sun sensor计算双向反射率，是一个非常值得研究的问题。因为这对于非直射光条件或者光照强烈变化条件下的双向反射率计算十分有用。我可能后面会专门写一篇博客介绍此事。
4. 双向反射率与植被指数的正射影像
我一般采用Agisoft Metashape软件进行正射影像的生成，原本这并不需要多说什么，但是对于获取DJI P4M的双向反射率和植被指数的正射影像，就需要多说几句了。
Metashape 能够自动识别原始DJI P4M影像，并且读取影像元数据包含EXIF metadata 和 XMP metadata，尤其是可以用XMP metadata中的一些参数提高三维重建的质量。而当原始影像被处理后，在用 exiftool.exe保存TIFF影像时，则会丧失所有的XMP metadata信息。不过倒是可以单独输入一个航点文件（该文件可以自行写脚本生成）需要注意的是，应当先自行写代码输出双向反射率影像，然后在Metashape中处理，进而得到正射影像。如果直接采用植被指数影像重建正射影像比较困难或者失败，则可以先生成双向反射率正射影像，然后再制作植被指数正射影像。
5. 后记
首先我有足够的信心认为，DJI P4M多光谱影像是高质量的可信的数据。但是其图像处理指南写的确实很不专业，这可能是因为P4M团队中并无遥感或者光学专家（我猜的）。其实通过与P4M团队的一个电话会议，我意识到DJI其实是想将P4M做的易用好用，而且一直在不断的努力中。
啥也不说了，等我有钱了，我要买台P4M （目前在借用北林骆校长团队的以及北师大阎老师团队的DJI系列无人机）。
