-
随着社会经济的发展，滩涂地区的开发利用逐步走向科学化、精细化，这对滩涂地区的监测与规划提出了更高要求。传统摄影测量和地面常规测量手段在陆地地形探测、地物分类方面有较好的应用实践，但滩涂地区潮汐变化较快，底质松软易下陷，难以布设有效的地面控制点，传统方法无法高效地获得高精度的测绘成果。
相比陆地地形地物，滩涂地区地形相对平缓，地物构成相对简单（通常只分布着泥、沙和水以及少量船只和建筑物等）。但准确的滩涂地物分类结果可为滩涂的合理开发与规划、登陆/抗登陆作战中登陆点的选取与评估[1]、瞬时水边线提取[2]、滩涂淤积和侵蚀的监控[3]以及潮沟发育演变的监测[4]等工作提供重要依据。
机载激光雷达（light detection and ranging，LiDAR）技术不仅能够获得高精度、高密度的滩涂三维空间信息，而且可以有效克服传统摄影测量和地面常规测量在滩涂测绘中的不足。但由于滩涂地势平坦，因此仅利用LiDAR数据中的高程信息难以实现滩涂地物的精细分类。而LiDAR系统采集的点云数据除了包含高精度的三维坐标信息外，还记录了回波强度信息。回波强度信息是地物光照反射的光强信息，由于不同性质的地物对光的吸收和折射特性不同，因此强度信息在一定程度上能反映地物的物理特性，可基于此对地物进行更为精细的分类。目前，已有利用LiDAR数据的强度信息对滩涂进行地物分类、表面分析和特征提取等方面的研究，并取得了较好的效果[2, 5]。
LiDAR数据的强度信息虽然可以辅助高程信息进行地物分类，但受目标属性、采集方式、仪器性能和测量环境等多种因素的影响，在使用前需要进行数据修正。目前，国内外提出的LiDAR强度修正方法[6-8]中大部分研究均假设LiDAR的强度符合朗伯反射模型。陆地地形中的城市和山地等地物由于分布较复杂，地面坡度范围变化较大，因此适用于该假设。但滩涂地形中的地势相对平坦、地表粗糙度相对较小，其LiDAR强度并不完全符合简单的朗伯反射模型。尤其在航带中心位置附近，由于激光入射角较小，使得LiDAR系统发射的激光产生强烈的镜面反射，造成回波强度在航带中心区域出现明显的高光现象。目前，针对LiDAR强度信息高光修正方面的研究较少[9-10]。本文通过实验分析发现，滩涂LiDAR强度信息中高光区分布不仅与激光入射角有关，还与飞机姿态角有关，基于此构建信息修正模型。
-
在假设扫描目标具有朗伯特性的条件下，可以将激光雷达测距方程简化为[11]：
$${P_r} = \frac{{{P_t}D_r^2\rho {\rm{cos}}\theta }}{{4{R^2}}}{\eta _{{\rm{atm}}}}{\eta _{{\rm{sys}}}}$$ (1) 式中，${P_r}$为接收到的光功率；${P_t}$为发射的光功率；${D_r}$为接收孔径；$R$为测距，即传感器与目标点之间的距离；ρ为目标反射率；θ为激光入射角，即激光脉冲方向与目标表面法向之间的夹角；${\eta _{{\rm{atm}}}}$为单程大气传输因子；${\eta _{{\rm{sys}}}}$为光学系统传输因子。
按文献[11]方法，进一步假设扫描系统对接收功率线性放大，并对测距R和入射角θ分离变量，可将强度修正至参考距离和参考角度，使修正后的强度只反映目标反射率，即：
$${I_c} = \frac{I}{{{f_2}\left( \theta \right){f_3}\left( R \right)}}{f_2}\left( {{\theta _s}} \right){f_3}\left( {{R_s}} \right)$$ (2) 式中，${I_c}$为修正后的强度；I为LiDAR数据的原始强度值；${\theta _s}$为参考角度；${R_s}$为参考测距；f2(θ)=cosθ；f3(R)=R-2×10-2aR/10 000，a表示大气衰减系数。为便于计算，一般标定${R_s}$为飞机平均飞行高度；${\theta _s}$为0°，即垂直入射。式（2）也称为归一化校正模型[12]，其具体表达式有很多种[6-14]，有些忽略单程大气传输因子的影响[7]，有些忽略入射角的影响[8]，还有一些增加了大气透射比、传输能量等因素的影响[13-14]，但各种方法之间在形式上非常接近，只是假设条件不同。
-
本文以经典的强度修正方法为基础，综合考虑测距、入射角和大气密度的影响，并从原始公式出发，分别通过独立地修正测距影响和入射角影响达到修正强度的目的：
$${I_c} = I\frac{{{f_2}\left( {{\theta _s}} \right){f_3}\left( {{R_s}} \right)}}{{{f_2}\left( \theta \right){f_3}\left( R \right)}} = {I_d}\frac{{{f_2}\left( {{\theta _s}} \right)}}{{{f_2}\left( \theta \right)}} = {I_d}\frac{1}{{{\rm{cos}}\theta }}$$ (3) 式中，${I_d}$为消除测距影响后的强度，只依赖于目标反射率和入射角：
$${I_d} = I\frac{{{f_3}\left( {{R_s}} \right)}}{{{f_3}\left( R \right)}} = I\frac{{{R^2} \times {{10}^{2aR/10{\rm{\;}}000}}}}{{R_s^2 \times {{10}^{2a{R_s}/10{\rm{\;}}000}}}}$$ (4) 式中，a为大气衰减系数。
利用式(4)的经典强度修正方法虽然可以消除环境因素的影响，但该方法需要假设扫描目标具有朗伯特性，即近似认为目标物为理想的漫反射面，且各个散射方向的辐射强度相同。但由于激光照明的单色性和连贯性，LiDAR应用中几乎不会发生理想的漫反射情况[15-16]。特别是对于滩涂这种地势平坦、地表潮湿平整的区域，很容易产生镜面反射现象。因此当漫反射和镜面反射同时存在时，朗伯反射模型就不再适用[10]。目前，常见的光照模型将光分为4部分：物体自发光、环境光、漫反射光和镜面反射光。对于LiDAR系统，物体自发光和环境光可忽略不计，因此，在排除测距和大气密度影响的条件下，光照模型可表示为：
$${I_d} = {I_{{\rm{diff}}}} + {I_{{\rm{spec}}}}$$ (5) 式中，${I_{{\rm{diff}}}}$为理想漫反射分量，可用朗伯模型进行拟合；${I_{{\rm{spec}}}}$为镜面反射分量，可用多种模型进行拟合。文献[10]中已经对该项的各种拟合模型的适应性进行了分析，并认为Torrance-Sparrow模型匹配效果最好。但利用该模型进行全局性修正时，在参数确定和反射率获取上存在一定的难度。
而Phong光照模型作为一种经典理想光照模型，由于需确定的参数较少，方法简单直观，被广泛应用于各领域。本文采用该模型对${I_{{\rm{spec}}}}$分量进行拟合。Phong光照模型的基本形式为：
$$ {I_d} = {I_{{\rm{in}}}}\left[ {{k_d}\left( {\mathit{\boldsymbol{l}} \cdot \mathit{\boldsymbol{n}}} \right) + {k_s}{{\left( {\mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{r}}} \right)}^n}} \right] $$ (6) 式中，${I_{{\rm{in}}}}$为LiDAR发射的光强度；${k_d}$为漫反射系数；${k_s}$为镜面反射系数（${k_s} + {k_d} = 1$）；l为目标指向光源的向量；n为目标表面的法线方向；v为视线方向；r为光反射方向。
式(6)中等式右边第一项即为漫反射分量，第二项为镜面反射分量。对于LiDAR系统，光信号发射位置和接收位置相同，即v与l共线，因此文献[9-10]中认为v·r = cos2θ，由此可得：
$${I_d} = {I_{{\rm{in}}}}\left[ {\left( {1 - {k_s}} \right){\rm{cos}}\theta + {k_s}{\rm{co}}{{\rm{s}}^n}2\theta } \right]$$ (7) $${I_{{\rm{spec}}}} = {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta $$ (8) 式(7)的光照模型是以θ=0°为轴对称分布的，即强度值总体上应以航带中心为轴呈对称分布。但通过观察实际数据发现，在航带中心附近，即使对于同种地物，其强度值有时也不是以航带中心为轴对称分布，如图 1所示。在此情况下，如果仍以入射角θ为基础进行高光修正，会出现强度修正中心位置偏离高光区位置的修正偏移现象，因此需要进一步分析产生修正偏移的原因，并对入射角进行适当校正，使高光现象得到合理的修正。
图 1 入射角分布和强度叠加图
Figure 1. Overlay of Incidence Angle Distribution and Intensity Image
由于${I_d}$已经消除测距和大气密度的影响，因此${I_d}$只与目标反射率和入射角相关。通过比对航空影像数据得知，图 1中航带中心区域为同一种地物，即该区域内各点的目标反射率应相同，因此可以进一步排除目标反射率对高光区分布的影响，从而推断高光区域偏离航带中心只与入射角有关。但由于滩涂地区地势平坦，地面坡度接近于0°，可以忽略地面坡度和地物倾角对强度的影响，因此又可以进一步排除地物的倾角和趋势角对高光区分布的影响。
反观视点一端，即LiDAR系统的激光发射端或反射信号的接收端，会随着飞机姿态的不断变化而改变位置。而通过观察发现，高光区偏离航带中心的现象一般发生在飞机姿态角在短时间内产生剧烈变化之时（图 2）。因此可解释为飞机姿态发生变化时，传感器的接收视角发生变化，可导致激光扫描镜编码器对天底位置的误判，从而对高光区域的分布产生影响（图 3）。对于城市、山地这些地面坡度变化较大的地形，飞机姿态的影响可以忽略。但对于滩涂这种地势平坦的地区，飞机姿态角的微小变化就会引起强度分布的变化。
图 2 飞机姿态角与强度分布关系视图
Figure 2. Distribution of the Aircraft Attitude Angle and the Intensity Image
图 3 LiDAR系统的光照模型参数视图
Figure 3. Illustration of Reflectance Parameters for LiDAR System
由于LiDAR系统大都具备姿态角自动补偿功能，在飞机姿态角变化较小时，这种高光区偏离的现象会被自动修正。若姿态角在短时间内剧烈波动则会产生补偿迟滞，从而出现高光区偏离航带中心位置的现象。为避免单纯利用入射角进行高光修正可能出现的修正偏移问题，本文利用视点端的飞机姿态角对目标端的激光入射角进行校正，以补偿强度修正中心的位置偏移量。当高光区偏离航迹中心时，将式(7)后一项中目标法向与反射光方向夹角θ校正为飞机自身坐标系Z轴负方向与反射光方向之间的夹角θ′，见图 4。θ′为三维立体角，可以通过坐标转换得到。θ′与θ之间的偏差虽很小，但对于高光区域的分布影响很大，这一观点已得到了后续实验结果的验证。
图 4 坐标变换示意图
Figure 4. Coordinate Transformation Diagram
如图 4所示，设大地坐标系和飞机坐标系分别为$O - XYZ$和$O{\rm{'}} - X{\rm{'}}Y{\rm{'}}Z{\rm{'}}$，O点为地球质心，O′点为飞机重心，P点为目标点，${\theta _R}$、${\theta _P}$、${\theta _H}$分别对应航迹信息记录的飞机侧滚角、俯仰角和偏航角，O′P与Z′O′的夹角即为θ′。通过LiDAR数据的坐标信息和飞机航迹信息可得O′点和P点在大地坐标系的坐标分别为（${\rm{\Delta }}X, {\rm{\Delta }}Y, {\rm{\Delta }}Z$）和（$X, Y, Z$）。为确定角度θ′，需将P点转换至飞机坐标系下，并设P点在飞机坐标系下的坐标为（$X{\rm{'}}, Y{\rm{'}}, Z{\rm{'}}$）。这两个坐标系可通过平移和旋转变换取得一致，坐标间转换关系如下：
$$ \left[ {\begin{array}{*{20}{c}} {X{\rm{'}}}\\ {Y{\rm{'}}}\\ {Z{\rm{'}}} \end{array}} \right] = {\mathit{\boldsymbol{M}}^{\rm{T}}}\left( {\left[ {\begin{array}{*{20}{c}} X\\ Y\\ Z \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{\rm{\Delta }}X}\\ {{\rm{\Delta }}Y}\\ {{\rm{\Delta }}Z} \end{array}} \right]} \right) $$ (9) 式中，M为旋转矩阵，
$$\begin{array}{*{20}{l}} {\mathit{\boldsymbol{M}} = \mathit{\boldsymbol{M}}\left( {{\theta _H}} \right)\mathit{\boldsymbol{M}}\left( {{\theta _P}} \right)\mathit{\boldsymbol{M}}\left( {{\theta _R}} \right) = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}{\theta _H}}&{{\rm{sin}}{\theta _H}}&0\\ { - {\rm{sin}}{\theta _H}}&{{\rm{cos}}{\theta _H}}&0\\ 0&0&1 \end{array}} \right] \cdot }\\ {\;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {{\rm{cos}}{\theta _P}}&0&{ - {\rm{sin}}{\theta _P}}\\ 0&1&0\\ {{\rm{sin}}{\theta _P}}&0&{{\rm{cos}}{\theta _P}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{{\rm{cos}}{\theta _R}}&{{\rm{sin}}{\theta _R}}\\ 0&{ - {\rm{sin}}{\theta _R}}&{{\rm{cos}}{\theta _R}} \end{array}} \right]} \end{array}$$ (10) 在飞机坐标系下，Z′O′方向的单位向量为${\mathit{\boldsymbol{e}}_{Z{\rm{'}}}} = {\left[ {\begin{array}{*{20}{c}} 0&0&{ - 1} \end{array}} \right]^{\rm{T}}}$，$\mathit{\boldsymbol{O}}{\rm{'}}\mathit{\boldsymbol{P}} = {\left[ {\begin{array}{*{20}{c}} {X{\rm{'}}}&{Y{\rm{'}}}&{Z{\rm{'}}} \end{array}} \right]^{\rm{T}}}$，则：
$$\theta {\rm{'}} = {\rm{arccos}}\left( {\left( {\mathit{\boldsymbol{O}}{\rm{'}}\mathit{\boldsymbol{P}} \cdot {\mathit{\boldsymbol{e}}_{Z{\rm{'}}}}} \right)/\left| {\mathit{\boldsymbol{O}}{\rm{'}}\mathit{\boldsymbol{P}}} \right|} \right)$$ (11) 将θ替换为θ′，并结合式(5)和式(7)对强度进行修正，以消除镜面反射对强度的影响，从而只保留漫反射部分的光强度，即：
$$I{{\rm{'}}_d} = {I_{{\rm{diff}}}} = {I_d} - {I_{{\rm{spec}}}} = {I_d} - {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta {\rm{'}}$$ (12) 考虑到LiDAR系统中姿态角自动补偿功能的存在，并非时刻都会出现高光区偏离的现象，因此本文将高光区中心是否与航带中心重合作为角度校正的限制条件。对高光区中心偏离航带中心区域的情况采用校正后的入射角θ′进行高光修正，即按照式(12)进行修正。而对于未发生高光区偏离的区域，则以式(7)为基础进行修正。而界定高光区是否偏离航带中心的步骤如下：
1）强度信息中广泛存在的椒盐噪声容易造成某些点的强度值异常偏高。为避免这些点影响高光区中心位置的判读，本文先通过中值滤波消除椒盐噪声的影响，将每个点i的强度和离该点最近的4个邻点的强度均值作为i点的强度值：
$${\bar I_i} = \left( {\mathop \sum \limits_{j \in Adj\left( i \right)} {I_j} + {I_i}} \right)/5$$ (13) 式中，Adj(i)为点i的邻域点集。
点i的最近邻可以通过对点云构建K维树获得。${\bar I_i}$只是为后续确定高光区中心位置而使用的一个中间值，不需要将${\bar I_i}$赋值给点i。
2）根据.las（LiDAR数据交换标准格式）数据中的全球定位系统时间和扫描方向识别标记，将整条航带的点云数据分割为N条测线。由Phong模型可知，高光区中心的强度最高，因此每条测线中强度最高的点即为该条测线高光区的中心点。虽然步骤1）中已经消除了椒盐噪声的影响，但仍有一些特殊情况影响高光区中心点位置的确定。例如条带边缘处存在船只或建筑等目标物时，覆盖这些目标物的测线上，强度最高的点为目标物上的反射点，而该点并不位于航带中心区域。因此仅比较每条测线上中心区域点的强度值，以确定每条测线高光区的中心点。综合考虑镜面反射实际影响范围和数据运算量，本文设定入射角小于5°的点均属于航带中心区域的点。此外，若特殊地物目标（如船只、建筑和植被等）出现在航带中心区域，且该地物目标的强度很高，则会影响高光区中心点的确定。针对这种特殊情况，可利用高程信息剔除每条测线中心区域的特殊地物目标点。由于滩涂区域地势平坦，每条测线中的地面点应大致呈一条直线，代表地物目标的点高程往往高于地面点。因此本文利用随机采样一致性（random sample consensus，RANSAC）算法将每条测线上中心区域的地物目标点作为直线模型的局外点筛选出来，在后续强度比较过程中不考虑这些点的影响。剔除航带中心区域地物目标点后，比较每条测线中心区域上各点的${\bar I_i}$值，其值最大的点即为该测线的高光区中心点。设第s条测线中心区域的点中${\bar I_i}$值最大的点所对应的入射角为$\theta _{{{\bar I}_{{\rm{max}}}}}^s$，则该点对应这条测线上高光区的中心位置。另外，设该条测线上各点入射角的最小值为$\theta _{{\rm{min}}}^s$，则该入射角所对应的点即为该条测线的航迹中心位置。若$\theta _{{{\bar I}_{{\rm{max}}}}}^s$与$\theta _{{\rm{min}}}^s$之差超过某阈值${\theta _T}$（本文设定${\theta _T}$为0.5°），则说明该条测线上高光区的中心位置偏离航带中心，此时应选用校正后的入射角θ′对整条测线上点的强度值进行高光修正。反之，若$\theta _{{{\bar I}_{{\rm{max}}}}}^s$与$\theta _{{\rm{min}}}^s$之差小于${\theta _T}$，则直接用入射角θ进行高光修正：
$$I{{\rm{'}}_d} = \left\{ {\begin{array}{*{20}{c}} {{I_d} - {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta {\rm{'}}, \left| {\theta _{{{\bar I}_{{\rm{max}}}}}^s - \theta _{{\rm{min}}}^s} \right| \ge {\theta _T}}\\ {{I_d} - {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta , {\rm{\;}}\left| {\theta _{{{\bar I}_{{\rm{max}}}}}^s - \theta _{{\rm{min}}}^s} \right| < {\theta _T}} \end{array}} \right.$$ (14) 由于${I_{{\rm{in}}}}$为标量，当θ=0°时，${I_{{\rm{in}}}} = {I_d}$，因此${I_{{\rm{in}}}}$可以通过计算${I_d}$的最大值获得。文献[9]中已对${k_s}$和n这两个参数对强度的影响进行了分析，${k_s}$影响漫反射光和镜面反射光之间的比率，n影响高光区的范围，因此式（14）中${k_s}$和n可通过数据分析和曲线拟合确定。
经过式（4）和式（14）消除测距和镜面反射的影响后，强度仅与入射角相关。因此将修正后的$I_d^{\rm{'}}$代入式（3），消除入射角的影响，从而获得最终的修正强度${I_c}$：
$${I_c} = \left\{ {\begin{array}{*{20}{c}} {\left( {{I_d} - {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta {\rm{'}}} \right)/{\rm{cos}}\theta , \left| {\theta _{{{\bar I}_{{\rm{max}}}}}^s - \theta _{{\rm{min}}}^s} \right| \ge {\theta _T}}\\ {\left( {{I_d} - {I_{{\rm{in}}}}{k_s}{\rm{co}}{{\rm{s}}^n}2\theta } \right)/{\rm{cos}}\theta , {\rm{\;}}\left| {\theta _{{{\bar I}_{{\rm{max}}}}}^s - \theta _{{\rm{min}}}^s} \right| < {\theta _T}} \end{array}} \right.$$ (15) 本文方法利用MATLAB软件编程实现，利用lastools工具对.las格式的点云数据进行读取和写入，利用Terrasolid软件对.trj格式的航迹数据进行读取和转换，ArcGIS软件进行数据的显示和分析。
-
本文方法的实验环境为Intel(R) Core(TM) i3-4130处理器，3.40 GHz主频，4 GB内存。本文使用的实验数据采集于中国江苏省盐城市（121°08′ E，33°01′ N）的滩涂地形区，设备为Leica ALS60，激光波长为1 064 nm，脉冲频率为110 kHz，扫描频率为41 Hz，飞行高度约为2 200 m，扫描角范围为±32°，数据参数见表 1。
表 1 数据参数列表
Table 1. Date Specification
参数 项目 数据1 数据2 数据
参数点数/个 1 010 677 1 016 797 覆盖面积/km2 2.77 3.17 平均点密度/（个·m-2） 0.69 0.64 强度范围 0~255 0~255 高程范围/m 8.56 ~ 9.91 7.15~16.7 航迹参数 侧滚角范围 -2.06°~-0.08° -2.94°~0.67° 俯仰角范围 0.83°~1.49° 0.80°~1.43° 航偏角范围 -12.25°~-11.72° -14.16°~-12.62° -
以数据1为例，将θ′与θ的分布图分别叠加在强度图像上，并只截取角度小于5°的部分（图 5）显示，可看出θ′在小角度分布区域和强度图像的高光区域基本重合（潮沟中水体部分除外）。本文结合文献[9]中的分析，利用MATLAB软件对式(7)进行曲线拟合，从而确定参数${k_s}$和n。数据1和数据2的最优${k_s}$和n分别为0.95、50和0.70、150。
图 5 数据1强度图像及角度分布图
Figure 5. Angle Distribution and Intensity Images of Data One
由图 6(c)和图 7(c)可看出，利用本文方法对强度进行修正后，高光现象基本消除。且对于同一地物，修正后的高光区强度值与高光区周围的强度值较为接近，体现出较好的同质性。作为对比，本文也利用未校正过的θ，按照式(12)的形式对高光现象进行修正，由图 6(b)和图 7(b)可以看出，由于角度的偏差，高光区的强度并未得到有效修正，反而航带中心位置的强度被过度修正，强度值明显过低（颜色偏暗）。
图 6 数据1的强度图像
Figure 6. Intensity Images of Data One
图 7 数据2的强度图像
Figure 7. Intensity Images of Data Two
对于数据1和数据2中航带中心区域的水体，本文的修正方法并没有完全消除高光现象。这是由于在水面上强度几乎由纯粹的镜面反射机制控制[10]，而Phong模型只有在${k_s} = 1$的情况下才符合这种假设。本文的修正方法是一种全局修正，且${k_s} \ne 1$，这就造成水体的镜面反射不能被完全抑制。但这恰恰提高了水陆强度之间的对比度，为水陆分类提供了便利。
为更清晰直观地评定修正效果，本文还对数据进行剖面采样，构建了强度的剖面散点图。以数据1为例进行显示，即图 6(a)中黄色区域。由图 8(c)可看出，利用本文方法可有效修正高光现象。而直接利用θ进行高光修正时，由于θ的分布与高光区域存在偏差，造成在红色箭头区域存在修正偏移现象，强度明显偏低，如图 8(b)所示。
图 8 数据1强度剖面散点图
Figure 8. Intensity Profiles of Data One
-
为进一步评价修正效果，本文借鉴文献[6]的定量分析方法，引入统计量变异系数（coefficient of variation, CV）和方差均值比（variance-to-mean, VRM）来评估同种地物（wi）回波强度的同质性。本文中结合同一飞行任务期间获得的遥感影像数据，在LiDAR点云中人工选出4个区域，分别代表 3类地物：干沙（区域1）、湿沙（区域2、区域3）、水体（区域4），见图 9。每个区域的原始强度和修正后强度的CV和VRM见表 2。
图 9 同种地物点采样区域分布图
Figure 9. Distribution of Homogeneous Sample Areas
表 2 不同方法的性能对比
Table 2. Capability Contrast of Different Methods
地物编号 点数/个 原始图像 直接利用θ进行修正 本文方法修正 CV VRM CV VRM CV VRM 1 20 548 0.254 2.482 0.244 2.869 0.227 2.445 2 48 615 0.474 10.452 0.548 6.910 0.330 2.960 3 30 813 0.517 9.751 0.468 8.306 0.352 4.213 4 11 504 0.737 18.359 0.733 18.362 0.729 17.846 由表 2可知，直接利用θ进行高光修正后，部分区域强度值的CV和VRM相比于原始强度图像虽有小幅下降，但区域2的CV和区域1、4的VRM非但没有降低，反而升高了。这是由于θ的分布与高光区域存在偏差，造成强度最高的区域得不到应有修正，反而部分区域修正过度(见图 8(b))，从而使本应属于同一地物的激光脚点的强度值差异变得更大，降低了强度的同质性。利用本文方法进行强度修正后，相比于原始强度图像，4个区域强度的CV和VRM均有不同程度的降低。区域3的CV最多可降低32.0%，区域2的VRM可下降71.7%。利用本文方法修正后，水体强度的CV和VRM仍然很高，相较原始强度下降较小，这主要是由于水体表面基本不存在漫反射，本文的全局式修正方法不能完全抑制水体的镜面反射，这也印证了本文之前的分析。本文还利用上述代表不同地物的采样点构建强度频数直方图，对比了原始强度、直接利用θ进行修正后的强度和本文方法修正后强度的频数直方图，如图 10所示，图中横坐标括号内数值分别表示直方图左右强度值。以区域2为例，直接利用θ进行强度修正后的强度频率直方图出现了双峰(见图 10(b))，这与文献[6]中所述的同类物质的强度频率直方图呈高斯分布明显不符。而本文方法修正后的强度频数直方图更接近于高斯分布，也进一步证明本文方法修正后的强度数据具有很好的同质性。
图 10 区域2强度直方图
Figure 10. Intensity Histograms of Region Two
为评估本文强度方法对地物分类的影响，文中还以数据1为例，采用最大似然估计法进行地物分类，并采用3种分类方案：①单纯利用数字表面模型（digital surface model，DSM）数据；②单纯利用强度（Intensity）数据；③联合利用DSM和Intensity数据。方案②和③中用到的强度数据分别采用原始强度数据（original intensity，OI）、利用θ修正后的强度数据（Iθ）和利用本文方法修正后的强度数据（Ip），各种分类方案的分类精度见表 3。比较表 3中各行数据可看出，相比于单纯利用强度或单纯利用DSM进行地物分类的方法，联合DSM和强度方法的不同地物分类精度均有所提升，这说明滩涂地区LiDAR数据的地物分类过程中高程信息和强度信息都是不可或缺的。进一步比较表 3最后一行，可以看出，相比于DSM+OI和${\rm{DSM}} + {I_\theta }$，利用${\rm{DSM}} + {I_p}$进行地物分类的精度更高，可以达到85.9%。说明利用本文方法修正后的强度信息可以有效地辅助高程信息进行地物分类，以提高整体地物的分类精度。
表 3 不同地物分类方案的分类精度/%
Table 3. Results of Classification Accuracy Obtained by Using Different Classification Scenarios/%
分类方案 DSM+OI DSM+Iθ DSM+Ip ① 46.9 46.9 46.9 ② 60.3 64.2 83.6 ③ 61.1 78.2 85.9 本文方法与文献[6, 12]的强度修正方法类似，同样都是利用矩阵变换对光照模型中的角度进行修正，不同点在于上述文献中考虑的是目标端地物的倾角和趋势角，而本文方法考虑的是视点端飞机姿态角，因此本文方法与上述文献中方法的时间复杂度相当。数据1和数据2的方法用时分别为51.64 s和49.89 s，对于百万量级的数据点，方法的时间效率在可接受范围内。
-
本文从现有LiDAR强度修正算法出发，通过对飞机姿态角角度变化和强度信息中高光区分布相关性的分析，针对滩涂地区的机载LiDAR数据，提出了顾及飞机姿态角的强度修正算法，较好地解决了单纯利用入射角进行高光修正而可能出现修正偏移的问题。通过视觉判读和定量分析证明此方法合理可行，能对滩涂LiDAR强度数据进行有效的修正。与直接利用激光入射角进行强度修正的方法相比，本文方法修正后的强度数据具有更好的同质性，可以为滩涂地形LiDAR点云分类、表面分析和特征提取提供更有效的辅助数据。通过比对不同强度修正方法对地物分类结果的影响，证明了本文方法修正后的强度信息可以有效提高地物分类精度，使滩涂地区的地物分类进一步精细化，这对于登陆/抗登陆作战、滩涂演变研究和瞬时水边线提取等工作都具有重要意义。
本文方法是一种全局式修正方法，算法中的参数ks和n需要借助于曲线拟合确定，对于大样本数据运算效率较低。并且本文方法是基于MATLAB软件编程实现的，方法的整体时间效率不高。因此下一步工作将致力于将该方法移植至C++平台，并尝试融入并行算法，以提高整体运算效率。此外，为验证本文强度修正方法对地物分类的影响，文中仅选用了较为经典的最大似然估计算法进行地物分类，后续还将提出更加适用于滩涂地区的地物分类方法，以进一步提高滩涂地区LiDAR数据的地物分类精度。
An Intensity Correction Method for Shoaly Land LiDAR in Consideration of the Attitude Angles of the Aircraft
-
摘要: 在深入分析现有机载激光雷达（light detection and ranging，LiDAR）强度修正算法在滩涂地区适用性的基础上，针对单纯利用激光入射角对滩涂LiDAR强度信息进行高光修正有可能造成强度修正中心位置偏移的问题，利用视点端的飞机姿态角对目标端的激光入射角进行精密校正，结合校正后的入射角和经典的Phong光照模型，提出了一种顾及飞机姿态角的滩涂LiDAR强度修正模型，通过对经典强度修正后的数据进行进一步的高光修正，定量补偿了强度修正中心的位置偏移量，实现了强度信息中高光现象的有效消除。在MATLAB平台下对方法的正确性和有效性进行了实验验证，结果表明，经由该方法修正后的强度信息具有更好的同质性，可显著提高地物分类的准确性。Abstract: Only by using the laser incidence angle to correct the intensity of shoaly land light detection and ranging (LiDAR) data, may cause the offset of the center of the intensity-corrected area. Therefore, an intensity correction method for shoaly land LiDAR in consideration of the attitude angles of the aircraft is proposed, based on the analysis of the existing LiDAR intensity correction method in the shoaly land area. The corrected incidence angle, which is obtained by precisely correcting the laser incidence angle at the target end according to the attitude angle at the viewpoint end, together with the classical Phong illumination model, is combined into our proposed method. After further highlight correction of the intensity data corrected by the classical intensity correction, the offset of the center of the intensity-corrected area is quantitively compensated, and the highlight of the intensity is effectively eliminated. The matrix laboratory platform is employed to validate the correctness and validity of the method. The results indicate that the intensity obtained by this method has better homogeneity, and the overall accuracy of the land-cover classification is improved significantly.
-
Key words:
- LiDAR /
- shoaly land /
- intensity correction /
- attitude angles of the aircraft /
- Phong illumination model
-
图 2 飞机姿态角与强度分布关系视图
Figure 2. Distribution of the Aircraft Attitude Angle and the Intensity Image
图 3 LiDAR系统的光照模型参数视图
Figure 3. Illustration of Reflectance Parameters for LiDAR System
表 1 数据参数列表
Table 1. Date Specification
参数 项目 数据1 数据2 数据
参数点数/个 1 010 677 1 016 797 覆盖面积/km2 2.77 3.17 平均点密度/（个·m-2） 0.69 0.64 强度范围 0~255 0~255 高程范围/m 8.56 ~ 9.91 7.15~16.7 航迹参数 侧滚角范围 -2.06°~-0.08° -2.94°~0.67° 俯仰角范围 0.83°~1.49° 0.80°~1.43° 航偏角范围 -12.25°~-11.72° -14.16°~-12.62° 
下载: 导出CSV
表 2 不同方法的性能对比
Table 2. Capability Contrast of Different Methods
地物编号 点数/个 原始图像 直接利用θ进行修正 本文方法修正 CV VRM CV VRM CV VRM 1 20 548 0.254 2.482 0.244 2.869 0.227 2.445 2 48 615 0.474 10.452 0.548 6.910 0.330 2.960 3 30 813 0.517 9.751 0.468 8.306 0.352 4.213 4 11 504 0.737 18.359 0.733 18.362 0.729 17.846
下载: 导出CSV
表 3 不同地物分类方案的分类精度/%
Table 3. Results of Classification Accuracy Obtained by Using Different Classification Scenarios/%
分类方案 DSM+OI DSM+Iθ DSM+Ip ① 46.9 46.9 46.9 ② 60.3 64.2 83.6 ③ 61.1 78.2 85.9
下载: 导出CSV
-
[1] 李占桥.海洋遥感在潮滩动态预测中的应用[D].青岛: 中国海洋大学, 2008 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1337921 Li Zhanqiao. The Application of Ocean Remote Sensing in Tidal Beach Dynamic Forecasting[D]. Qingdao: Ocean University of China, 2008 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1337921 [2] Schmidt A, Rottensteiner F, Soergel U. Monitoring Concepts for Coastal Areas Using LiDAR Data[C]. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Strasbourg, France, 2013 [3] Zhou Liangyong, Wang Weiwei, Gao Maosheng, et al. Monitoring the Tidal Topography of the North Yellow River Delta with LiDAR Data[J]. Geological Bulletin of China, 2016, 35(10):1661-1668 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=zgqydz201610014 [4] Liu Yongxue, Zhou Minxi, Zhao Saishuai, et al. Automated Extraction of Tidal Creeks from Airborne Laser Altimetry Data[J]. Journal of Hydrology, 2015, 527(32):1006-1020 http://cn.bing.com/academic/profile?id=68bca9058b3a1a38d09b5ecc6601db08&encoded=0&v=paper_preview&mkt=zh-cn [5] Schubert J E, Gallien T W, Majd M S, et al. Terrestrial Laser Scanning of Anthropogenic Beach Berm Erosion and Overtopping[J]. Journal of Coastal Research, 2015, 31(1):47-60 doi: 10.2112/JCOASTRES-D-14-00037.1 [6] Yan W Y, Shaker A. Radiometric Correction and Normalization of Airborne LiDAR Intensity Data for Improving Land-Cover Classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(12):7658-7673 doi: 10.1109/TGRS.2014.2316195 [7] Habib A F, Kersting A P, Shaker A, et al. Geometric Calibration and Radiometric Correction of LiDAR Data and Their Impact on the Quality of Derived Products[J]. Sensors, 2011, 11(9):9069-9097 doi: 10.3390/s110909069 [8] Starek M, Luzum B, Kumar R, et al. Normalizing LiDAR Intensities[R]//GEM Center Report No. Rep_2006-12-001. Gainesville: University of Florida, 2006 [9] Ding Qiong, Chen Wu, King Bruce, et al. Combination of Overlap-Driven Adjustment and Phong Model for LiDAR Intensity Correction[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 75(1):40-47 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=d1ae22e097e2bebd51d574529e7e2d1b [10] Poullain E, Garestier F, Levoy F, et al. Analysis of ALS Intensity Behavior as a Function of the Incidence Angle in Coastal Environments[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(1):313-325 doi: 10.1109/JSTARS.2015.2510337 [11] Tan Kai, Liu Kunbo, Cheng Xiaojun. An Empirical Method in Correcting Specular Highlight Phenomenon in TLS Intensity Data[J]. IEEE Access, 2016, 4:9821-9827 doi: 10.1109/ACCESS.2016.2647559 [12] 伊丕源, 满旺, 童鹏, 等.加入地物倾角分析的LiDAR回波强度校正[J].遥感学报, 2016, 20(4):610-619 http://d.old.wanfangdata.com.cn/Periodical/ygxb201604010 Yi Piyuan, Man Wang, Tong Peng, et al. Calibration Algorithm and Object Tilt Angle Analysis and Calculation for LiDAR Intensity Data[J]. Journal of Remote Sensing, 2016, 20(4):610-619 http://d.old.wanfangdata.com.cn/Periodical/ygxb201604010 [13] Vain A, Kaasalainen S, Pyysalo U, et al. Use of Naturally Available Reference Targets to Calibrate Airborne Laser Scanning Intensity Data[J]. Sensors, 2009, 9(4):2780-2796 doi: 10.3390/s90402780 [14] Ahokas E, Kaasalainen S, Hyyppä J, et al. Calibration of the Optech ALTM 3100 Laser Scanner Intensity Data Using Brightness Targets[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences, 2006, 36:10-16 [15] Kaasalainen S, Hyyppa H, Kukko A, et al. Radiometric Calibration of LiDAR Intensity with Commercially Available Reference Targets[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(2):588-598 doi: 10.1109/TGRS.2008.2003351 [16] 张永军, 吴磊, 林立文, 等.基于LiDAR数据和航空影像的水体自动提取[J].武汉大学学报·信息科学版, 2010, 35(8):936-940 http://ch.whu.edu.cn/CN/abstract/abstract1012.shtml Zhang Yongjun, Wu Lei, Lin Liwen, et al. Automatic Water Body Extraction Based on LiDAR Data and Aerial Images[J]. Geomatics and Information Science of Wuhan University, 2010, 35(8):936-940 http://ch.whu.edu.cn/CN/abstract/abstract1012.shtml -
点击查看大图
图(10) / 表(3)
计量
- 文章访问数: 1082
- HTML全文浏览量: 64
- PDF下载量: 123
- 被引次数: 0
