利用无人机搭载的多光谱相机直接测定地表反射率的新方法

孙华生, 张 远, 史云飞, 赵 敏

1. 临沂大学资源环境学院(山东省水土保持与环境保育研究所)， 山东 临沂 276000

2. 华东师范大学地理科学学院， 崇明生态研究院， 地理信息科学教育部重点实验室， 上海 200241

引 言

无人机遥感具有操作简便、 数据获取成本低、 时空分辨率高、 可在云下实现数据采集等诸多优势， 弥补了卫星或大型航空遥感系统的缺点， 所以在小范围、 高时间和高空间分辨率遥感监测中， 发挥了极为重要的作用。 目前， 无人机遥感技术已被广泛应用于农业、 林业、 环境领域的定量监测[1]。 反射率是地物较为稳定的物理属性， 也是定量遥感模型的重要输入参数之一。 在遥感图像定量化应用过程中， 需要将相机记录的数字量化(digital number， DN)值转换为地表反射率。 因此， 得到准确的地表反射率是实现遥感数据定量化应用的基础。

在环境光照稳定的条件下， 通常只要在无人机起飞前或航拍完成后， 拍摄两个或更多具有不同反射率的参考板， 通过对每个波段记录的DN值进行线性拉伸， 即可得到每个波段的地表反射率[2-4]。 然而， 实际的状况却是异常复杂的[5-7]， 这是因为： (1)受云层遮挡和太阳高度角变化的影响， 数据获取时的光照条件可能会随时发生变化； (2)为了得到高质量的图像， 相机曝光时间、 光圈大小和ISO值， 需要根据拍摄时的光照条件自动调整； (3)在很多情况下， 在地面上铺设参考板是很难或无法实现的(如森林、 水域或沼泽地中)， 而且要想在每幅图像中都拍摄到一组参考板也是不现实的。

为了摆脱对参考板的依赖， 一些研究者采用经验线法得到近似的地表反射率[2]； 还有的建议在无人机上安装一个光强传感器， 来同步记录相机拍摄时每个波段对应的光照强度， 并将该测量结果用于计算地表反射率[3-5, 7-9]， 所以获取可靠的地面接收的太阳辐照度是得到准确地表反射率的前提。 如果光强传感器不是安装在一个精密的增稳平台上， 受飞行器姿态变化的影响， 导致传感器的测量结果与到达地面的太阳辐射强度并不一致。 对此， 一些研究者采用区块校正法来降低传感器的倾斜效应[8-10]， 但其得到的只是一种近似的结果。 还有采用余弦函数对直接测量的结果进行倾斜校正， 但实验测试表明， 余弦校正的结果并不理想[3, 6]。 这是因为： 太阳辐射不仅包括直射， 还包括散射， 在不同天气状况下， 太阳辐射的组成有着很大的差异， 而传感器倾斜对太阳直射和散射的影响是完全不同的。 为了获取更加准确的太阳辐照度测量结果， 可将光强传感器安装在一个自动调平装置上。 例如， Markelin等利用一个光学平衡系统， 可以对光强传感器进行15°以内的倾斜补偿[4]， 但对多旋翼无人机来说， 在飞行状态下， 其倾斜角度有时甚至会超过30°。 如果只是将光强传感器搭载在一个简易的调平装置上， 并不能完全消除飞行器倾斜造成的影响； 而如果单独为光强传感器安装一个精密的调平装置(如一个三轴增稳云台)， 将会大大增加飞行器的负担和制造成本。 在实际应用中， 在没有精密调平装置的条件下， 为了使光强传感器的测量结果更加可靠， 常用的做法包括： (1)要求航拍时太阳高度角较大， 并且航线与太阳方位角垂直， 以降低传感器倾斜对太阳辐照度测量结果的影响， 但这种方法并不能消除传感器倾斜的影响； (2)采用悬停方式进行拍摄， 但这种方式只有在没有风的情况下才能使传感器保持水平， 否则， 并不能完全消除传感器倾斜的影响， 而且其作业效率极低。 因此， 目前利用无人机光强传感器对太阳辐照度的测量， 仍然存在一定的缺陷， 或者只能在一定的限定条件下使用。

如何在倾斜状态下得到准确的太阳辐照度， 进而得到准确的地表反射率， 是无人机遥感领域多年来一直未能很好解决的技术难题。 本研究对该问题进行了深入探讨， 并提出了一种太阳辐照度准确测量的新方法。 该方法可以在任意光照条件下， 直接将多光谱图像记录的DN值转换为准确的地表反射率。

1 实验部分

1.1 无人机平台和传感器

利用大疆精灵4 RTK无人机作为遥感数据获取平台。 该无人机搭载了一台一体式的多光谱相机， 共有6个1/2.9 英寸的CMOS传感器(详见图1)， 包含1个可见光波段的RGB传感器和5个多光谱传感器。 5个多光谱波段参数的描述， 如表1所示。 此外， 在该无人机的顶部还集成了一个光强传感器， 可同步检测5个波段的太阳辐照度， 并将其记录在图像的XMP(Extensible Metadata Platform)元数据中[11]。

1.2 多光谱图像转换为地表反射率的方法

地表反射率的计算方法如式(1)所示

Rλ=Er(λ)/Eg(λ)

(1)

式(1)中，Rλ为波长为λ的地表反射率，Er(λ)为地物对波长为λ的电磁波在单位波长单位面积的反射辐射功率(W·m-2·nm-1)，Eg(λ)为到达地面的入射辐射辐照度(W·m-2·nm-1)。

对Er(λ)的定量描述如式(2)所示

(2)

式(2)中，L为每个波段每个像元对应的地物反射辐照度(W·m-2·Sr-1·nm-1)，θ为太阳辐射的天顶角，φ为太阳辐射的方位角。

在计算反射率时， 需要首先将DN值转换为对应像元的辐射亮度， 然后再将辐射亮度进一步转换为地表反射率。 由于无人机多光谱相机拍摄的原始图像记录的DN值会受曝光时间、 ISO值、 暗角效应等因素的影响， 所以需要事先对原始图像记录的DN值进行修正处理， 以得到统一标准下的结果。 具体转换方法如式(3)所示

DN′=(DN-DNblack)/N/Gsensor/Te×A×V

(3)

式(3)中，DN′为转换后的每个波段的数字量化值，DN为原始图像记录的每个波段的数字量化值，DNblack为黑电平的数字量化值，N为图像可表示的最大灰度值，Gsensor为传感器增益值，Te为曝光时间，A为传感器敏感度修正系数，V为相机镜头的暗角效应修正值。 以上参数都是在实验室中测定的， 其中DNblack，Gsensor，Te和A可以从XMP元数据中查找到， 而暗角效应修正值V可按照式(4)计算得到

V(x,y)=(1+k1r+k2r2+k3r3+k4r4+k5r5+k6r6)

(4)

修正后的DN′(λ)与每个波段的辐射亮度L(λ)之间满足简单的线性转换关系， 其表达式如式(5)所示

L(λ)=Gλ×DN′(λ)+Bλ

(5)

式(5)中，Gλ为增益值，Bλ为偏移值，DN′为修正后的数字量化值。

因此， 根据式(1)、 式(2)和式(5)， 地表反射率Rλ可表示为式(6)的形式

Rλ=Er(λ)/Eg(λ)=πL(λ)/Eg(λ)=

π[Gλ×DN′(λ)+Bλ]/Eg(λ)

(6)

如果能够准确地得到图像拍摄时到达地面的入射辐射辐照度Eg(λ)， 通过在一张图像中拍摄到两张或者更多具有不同反射率的参考板(其反射率是已知的)， 即可得到每个波段的定标参数Gλ和Bλ， 其计算方法如式(7)所示

DN′(λ)×Gλ+Bλ=RλEg(λ)/π

(7)

由此可见， 利用倾斜的光强传感器获取准确的Eg(λ)， 是需要解决的关键问题， 而这也是长期以来一直未能很好解决的技术难题。 在本研究中， 提出了一种实现太阳辐照度准确测量的新方法， 对其具体描述如下：

首先， 需要对光强传感器接收的太阳辐射进行定量化描述。 其接收的太阳辐射包括： 太阳直射、 天空散射和少量的地表反射。 对太阳直射来说， 其辐照度测量结果受传感器倾斜的影响， 如式(8)所示

(8)

(9)

(10)

根据以上分析可知， 在光强传感器平面倾斜时， 直接测量的每个波段的太阳辐照度Em的表达式， 如式(11)所示

(11)

根据式(11)可知， 如果能够通过某种方法获取朝向两个或者更多方向的太阳辐照度， 就可以直接求出直射Edir和散射Esca。 需要注意的是： 如果z≥90°， 即cos(z)≤0， 就会导致直射辐射无法到达光强传感器的正面， 从而使传感器正面处于阴影中。 此时， 需要令cos(z)=0。

在得到直射Edir和散射Esca之后， 地面接收的太阳辐照度Eg可表示为式(12)的形式

Eg=Edircosθ+Esca

(12)

在得出地面接收的太阳辐照度Eg后， 即可根据式(6)得到各个波段每个像素对应的地表反射率Rλ。

实际上， 如果能够得到直射辐射占总太阳辐射的比例p(p=Edir/Esun， 而Esun=Edir+Esca)， 那么就能将光强传感器测量值Em， 直接转换为地面接收的太阳辐照度Eg。 把Edir=pEsun，Esca=(1-p)Esun代入式(11)， 即可得到Em和Esun之间的关系式， 其结果如式(13)所示

Em=[pcosz+(1-p)cos2(s/2)+

(13)

因此， 地面接收的太阳辐照度Eg可以进一步转换为式(14)的形式

Eg=Edircosθ+Esca=(pcosθ+1-p)Esun=

(14)

1.3 实验测试

1.3.1 实验区

实验区位于山东省临沂市， 其地理位置为35°9′53.5″N—35°10′02.2″N， 118°16′00.5″E—118°16′08.5″E。 该实验区的面积约250 m×200 m。 航线规划是通过DJI GS Pro应用软件实现的， 设定的航高为35 m， 航线重叠度为67%， 旁向重叠度为55%， 共生成14条航线， 每条航线拍摄31张图像。 飞行速度约为3.8 m·s-1， 每次航拍约持续16 min。 从2020年7月开始， 一直持续到11月底， 共获取了4个时期的航拍图像。

1.3.2 相机定标

选取黑、 灰、 白三张参考板(尺寸为18 cm×25 cm)作为靶标， 以实现无人机多光谱图像的反射率定标和反射率测量结果的验证。 各个参考板在不同波段的反射率通过PSR-1100高光谱仪测量得到。 该高光谱仪的波谱覆盖范围为320～1 100 nm， 在600 nm处的光谱分辨率为3 nm。 黑、 灰、 白参考板及其高光谱测试结果， 如图2所示。 黑、 灰、 白三张参考板在不同波段的平均反射率， 如表2所示。

表2 利用高光谱仪测量的参考板的平均反射率

图2 黑、 灰、 白参考板及其高光谱测试结果

在得到无人机拍摄的参考板的多光谱图像之后， 根据式(7)构建Rλ和DN′(λ)之间的线性拉伸关系， 即可得到每个波段的Gλ和Bλ的结果。 对一台多光谱成像系统来说， 一旦得到其定标参数， 该定标参数通常是较为稳定的， 所以在实际应用中就不必再依靠参考板， 即可得到准确的地表反射率。

1.3.3 验证方案设计

如果飞行器自身携带两个或者更多的朝向不同方向的光强传感器， 即可根据式(11)实现直射Edir和散射Esca的求解。 然而， 由于本研究采用的无人机只有一个光强传感器， 而且是固定在无人机顶部的， 所以无法直接实现以上解算。 为了验证本方法的可行性， 设计了如下的实验方案：

首先， 选取一块开阔的实验场地， 确定其周围没有高大的树木或建筑物。 然后， 用手举起飞行器并手动调整其姿态， 以让光强传感器平面朝向不同的方向(在操作时可以使传感器平面分别处于水平、 前倾、 后倾、 左倾、 右倾的状态)， 并拍照以记录不同方向上各波段接收的太阳辐照度， 从而可实现对各个波段的直射Edir和散射Esca的求解。 为了检验得到的地表反射率的准确度， 在实验场中还布置了一组黑、 灰、 白参考板。

分别在不同日期不同光照条件下获取航拍数据， 而每次实验时的天气都较为稳定。 在这种情况下， 在一段较短的时间内， 可认为太阳辐射的组分是不变的(即使太阳辐照度有小幅度的变化， 但其组分的变化较小)。 因此， 可利用本方法， 得到在航拍前或航拍后的直射比例p， 并根据式(14)计算出所有图像在拍照时刻地面接收的太阳辐照度Eg， 从而可根据该结果得到所有图像的地表反射率。 最后， 利用参考板来验证计算结果的准确度。

2 结果与讨论

由于每次实验都是在光照条件较为稳定的情况下进行的(但不同日期的光照条件是不同的)， 所以航拍时段的直射比例可以在航拍前或航拍后， 利用本方法来确定。 航拍时段不同波段直射比例的计算结果， 如表3所示。

表3 航拍时段不同波段的直射比例p

在得到直射比例p之后， 通过式(14)即可将光强传感器直接测量的结果Em， 转换为地面接收的太阳辐照度Eg。 不同日期的具体计算结果如图3所示。

图3 光强传感器直接测量的太阳辐照度Em及对其进行校正后得到的地面接收的太阳辐照度Eg的结果

由图3可以看出， 光强传感器的直接测量结果受飞行器自身姿态变化(因航向发生了改变)的影响很大， 尤其是当太阳高度角较小(即入射角较大)时， 其影响极为显著。 利用本研究提出的方法可以很好地消除测量结果的剧烈波动。 为了进一步检验该校正结果的准确度， 利用事先铺设的黑、 灰、 白参考板， 来检验无人机多光谱图像得到的反射率的准确度， 其结果如表4所示。

表4 利用无人机多光谱图像得到的参考板在不同波段的平均反射率

利用平均绝对误差(MAE)， 以及标准差(S)两种指标对测量结果的准确度进行评价。 计算结果显示， 蓝(B)、 绿(G)、 红(R)、 红边(RE)和近红外(NIR)波段在所有日期所有参考板的MAE的结果分别为3.34%， 1.36%， 1.45%， 1.81%和1.63%， 而各自的标准差S分别为2.11%， 0.85%， 1.21%， 1.88%和0.59%。 黑、 灰、 白三张参考板在所有日期的所有波段中的MAE的结果分别为0.82%， 1.99%和2.94%， 而各自的标准差S分别为0.55%， 1.31%和1.84%。 引起地表反射率测量误差的因素主要包括以下几个方面： (1)参考板不是严格的漫反射体， 得到的反射率会因太阳入射方向和相机观测方向的不同而存在一定的误差； (2)参考板的高光谱测量结果， 以及得到的相机定标参数存在一定的误差； (3)光照条件不是绝对稳定的， 因为在航拍时段内太阳高度角和云量发生了一定的变化， 从而导致真实的直射比例p不是绝对稳定的。

尽管得到的地表反射率存在着一定的误差， 但从以上统计结果可以看出， 在不同光照条件下得到的不同波段的地表反射率， 仍然具有较高的准确度。 由此可以推断出， 利用本方法得到的地面接收的太阳辐照度Eg的结果是可靠的。

3 结 论

提出了一种全新的可以利用多光谱无人机图像实现地表反射率直接测量的方法。 与传统的基于光强传感器实现太阳辐照度测量方法不同的是： 本方法不是直接消除光强传感器倾斜的影响， 而是借助传感器的倾斜效应来实现对光照强度的校正。 即使在光照强度变化的条件下， 利用该方法仍然能够得到准确的测量结果， 从而可据此进一步得到准确的地表反射率。 本研究解决了多年来一直未能得到很好解决关键技术问题， 对多光谱无人机的设计， 以及无人机遥感数据的定量化应用， 都具有重要意义。 受硬件条件的限制， 本研究并没有直接利用无人机搭载的多个传感器， 来获取每个拍照时刻的太阳辐照度。 在未来的研究中， 将对飞行器进行硬件上的改进， 使其携带多个朝向不同方向的光强传感器， 以实现每个拍照时刻的直射和散射辐照度的准确测量， 并检验其实际测量效果。

附注1: 几个角度的计算方法

(1)太阳直射的天顶角θ可以通过式(A1)来确定

cos(θ)=sin(Ψ)sin(δ)+cos(Ψ)cos(δ)cos(H)

(A1)

式(A1)中，Ψ为某一点的纬度(弧度)；δ为太阳赤纬，δ=-23.44/180πcos(2π/365(N+10))，N为某日期在一年中的天数；H为时角(以南为参考)。H可以通过式(A2)计算出来

H=π/12(Tsolar-12)

(A2)

式(A2)中，Tsolar为地方太阳时， 即根据太阳的具体位置所确定的时刻， 与地方标准时间TLst的关系为Tsolar=TLst+12/π(Lonloc-Lonst)+E，Lonloc为某一点的经度(弧度)，Lonst为地方标准时间所采用的标准经度(弧度)，E为地球绕太阳公转时的进动和转速变化对地方太阳时的修正，E=[9.87sin(2B)-7.53cos(B)-1.5sin(B)]/60， 而B=2π(N-81)/364。

(2)传感器倾斜的坡度s可以利用式(A3)来确定

(A3)

式(A3)中，vZ为Z轴的方向向量(vZ=(0, 0, 1)T)；vn为斜面的法线向量，vn=RvZ(R为旋转矩阵， 可以利用俯仰角、 侧滚角和偏航角计算出来， 它们都记录在XMP文件中)。

(3)太阳直射在倾斜传感器坡面上的入射角z可以利用式(A4)来确定

cos(z)=cos(θ)cos(s)+sin(θ)sin(s)cos(φ-a)

(A4)

式(A4)中，θ为太阳直射的天顶角；s为传感器倾斜的坡度；φ为太阳方位角，a为传感器倾斜的坡向。 而φ和a可以分别通过式(A5)和式(A6)计算出来。 注：φ和a的参考方向需要与式(A2)中的时角H一致， 即以南为参考。

φ=arctan2(sinφ, cosφ)

(A5)

a=arctan2(sina, cosa)+π

(A6)