CN106291505A - 一种非植被覆盖区机载LiDAR数据回波强度值校正方法 - Google Patents

Abstract

本发明属于遥感测绘技术领域，具体涉及一种非植被覆盖区机载LiDAR数据回波强度值校正方法。本发明包括以下步骤：基于航高、激光点高程值进行初次校正回波强度值，将激光传输距离对回波强度值的影响均一化；基于初次校正后的LiDAR数据，以平面距离值、高程值、回波强度值为参考，通过设定阈值判定条件，逐点计算激光点的趋势角度值；结合趋势角和入射角，计算出反射角，重新对原始LiDAR数据的回波强度值进行校正计算。本发明解决了现有技术难以实现非植被类地物、尤其是倾斜地物的机载LiDAR数据回波强度值校正的技术问题，能够使机载LiDAR数据回波强度值更好地反映地物介质属性，提高其在地物分类识别方面的应用价值。

Description

一种非植被覆盖区机载LiDAR数据回波强度值校正方法

技术领域

本发明属于遥感测绘领域，涉及一种机载LiDAR数据回波强度值的校正方法。

背景技术

机载激光雷达(Airborne Light Detection And Ranging，LIDAR)是近十几年来快速发展的一种新型的测量技术，能够直接获取真实地面的高精度三维地形信息，可广泛用于大面积地形测绘、快速生成DEM等数字产品。激光点云数据则是激光测距仪对地发射激光，并通过计算发射激光的角度、脉冲信号返回时间等来计算获取的地面点三维坐标。除空间坐标信息外，激光雷达点云数据中记录的内容还包括回波强度、回波次序、扫描角度等信息。其中，回波强度值是经地物反射后返回的激光脉冲信号由接收器感应产生电压信号，经过放大处理和计算，转化为数字。

此外，激光点云密度是机载激光雷达点云数据质量的重要指标，如依据国家测绘地理信息局2011年颁布的《机载激光雷达数据获取技术规范》，要达到1:2000的成图比例尺，激光点云密度必须大于1点/平方米。

近年来随着数据采集技术的提高，激光回波强度信息精度大大提高，带有全波形记录仪的机载激光雷达设备也开始在市场上应用，因此越来越多的学者开始重视对回波强度信息的使用。国内外对融合激光回波强度信息的滤波分类算法也进行了一些研究。

理论上，每种物质对激光信号的反射特性是不一样的，根据这一特性，能够区分识别很多地物。地面介质表面的反射系数决定了激光回波能量的多少。地面介质对激光的反射系数取决于激光的波长、介质材料以及介质表面的明暗黑白程度。反射介质的表面越亮，反射率就越高。表1列出了德国Riegl公司测试的部分地物对900nm波长激光反射率。从表1中可以看出，即使同种材质的物体，其水分含量等不同，也会导致反射率差异。

表1部分介质900nm激光的反射率值

材质	反射率
		白纸	接近100％
形状规则的木料(干松木)	94％
		雪	80-90％
泡沫	88％
		白石块	85％
石灰石、粘土	接近75％
		有印记的新闻纸	69％
绵纸	60％
		落叶树	典型值60％
松类针类常青树	典型值30％
		碳酸盐类沙(干)	57％
碳酸盐类沙(湿)	41％
		海岸沙滩、沙漠裸露地	典型值50％
粗糙木料	25％
		光滑混凝土	24％
带小卵石沥青	17％
		火山岩	8％
黑色氯丁(二烯)橡胶	5％
		黑色橡皮轮胎	2％

激光强度信号虽然能够在一定程度上提供地物反射特性的信息，但实际上并不能很好地用来重构地面物体的反射性质，主要原因是激光回波强度不仅与反射介质的特性有关，还同激光的入射角度、激光脉冲作用距离产生的大气对激光的吸收等因素相关。这一缺点不仅制约了激光强度数据的精度，而且还使得根据强度数据进行地面物体分类变得困难。因此，如何校正激光回波强度值便成为关键问题。

要实现精确的激光回波强度值校正，理论上要依据最基本的雷达方程式，结合大气的散射、透过率等多项参数，才可以实现。但实际情况下，激光雷达设备的很多基本参数难以获取，因此雷达方程的应用较为困难，而且航空飞行过程中的精确大气参数获取同样很困难，所以大大限制了回波强度值的应用。

国外激光雷达硬件设备研制较早，在数据处理方面也开展了大量研究。其中，LiDAR回波强度值校正处理方法主要分为两类：基于理论的辐射校正模型和归一化校正模型。基于理论模型的辐射校正是基于雷达方程和激光传输理论，假定地面为朗伯体反射、大气条件已知(散射、透过率等)、激光发射功率恒定、接收功率与回波强度记录值线性相关，在满足以上条件下通过雷达方程计算公式进行回波强度值的校正。归一化校正模型，假定记录的回波强度值与地面反射率呈比例关系，与飞行高度呈线性关系或者平方反比，然后通过计算对距离值变化引起的差异进行归一化处理，通常计算公式为：I_{(R)normalized}＝I·R²/(Rs²·cosα)，其中I_{(R)normalized}为标准化后的回波强度值，I为原始回波强度测量值，R为激光点与传感器之间距离值，Rs为参考距离值，α为激光光束入射角，此外还有一种回波强度值计算公式：I_{(R)normalized}＝I·R²/Rs²，该方法同样假设地面为朗伯体，且地物平坦，忽略地物倾角对激光脉冲回波强度的影响，且通常只针对单次回波和首次回波激光点有效。然而，实际上地物倾角对激光回波强度值的影响不可忽略。

植被的回波强度校正难以解决，主要原因是植被自身的树叶形态对激光回波产生影响，但植被树叶形态的准确计算是十分困难的，其回波强度值与激光光斑面积、叶面积、叶倾角、树叶成分等都有关系，在有全波形记录仪数据的情况下可以开展一定工作。相对而言，非植被类地物的形态较为稳定，倾角容易计算(如房顶)，所以针对非植被类地物的回波强度校正，尤其是倾斜地物趋势角的影响分析和校正算法是可以研究的。

发明内容

本发明需要解决的技术问题：现有技术难以实现非植被类地物、尤其是倾斜地物的机载LiDAR数据的回波强度值校正。

本发明的技术方案如下所述：

一种非植被覆盖区机载LiDAR数据回波强度值校正方法，包括以下步骤：步骤1回波强度初步校正；步骤2确定激光点趋势角；步骤3回波强度最终校正。

步骤1中，首先采用依据下式计算高程平均值：

Z_平均＝(Z_最高+Z_最低)/2

式中，

Z_最高为高程最大值；

Z_最低为高程最小值；

Z_平均为高程平均值；

然后依据下式进行回波强度初步校正：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)，其中，R_j＝H-Z_j，R_A＝H-Z_平均；

式中，

H为机载LiDAR数据获取时的飞行航高；

Z_j为激光点P_j高程值，j＝1,…,n；

R_j为激光点P_j对地传输距离；

R_A为整个测区的n个激光点平均对地传输距离；

I_j原始值为激光点P_j回波强度原始值；

I_{j初步校正值}为激光点P_j回波强度初步校正值。

步骤2包括以下步骤：

步骤2.1分析平面距离

采用下式计算临近激光点平面距离：

$L_j=\sqrt{{(X_{\left(P_j\right)}-X_{\left(P_{j-1}\right)})}^2+{(Y_{\left(P_j\right)}-Y_{\left(P_{j-1}\right)})}^2}$

$R_j=\sqrt{{(X_{\left(P_{j+1}\right)}-X_{\left(P_j\right)})}^2+{(Y_{\left(P_{j+1}\right)}-Y_{\left(P_j\right)})}^2}$

式中，

L_j为激光点P_j与相临近的P_j-1的平面距离；

R_j为激光点P_j与相临近的P_j+1的平面距离；

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；

设定平面距离阈值，然后将L_j、R_j与平面距离阈值进行比较：

若L_j、R_j均小于平面距离阈值，进行步骤2.2.1；

若只有R_j小于平面距离阈值，进行步骤2.2.2；

若只有L_j小于平面距离阈值，进行步骤2.2.3；

若L_j、R_j均大于平面距离阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；

步骤2.2设定高差阈值，分析高程值

步骤2.2.1若激光点P_j的高程值同时大于其临近激光点P_j-1和P_j+1的高程值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，分别计算激光点P_j与激光点P_j-1的高程差、激光点P_j与激光点P_j+1的高程差，若上述两个高程差均大于高程差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个高程差均小于高程差阈值，则进行步骤2.3.1；若只有激光点P_j与激光点P_j+1的高程差小于高程差阈值，则进行步骤2.3.2；若只有激光点P_j与激光点P_j-1的高程差小于高程差阈值，则进行步骤2.3.3；

步骤2.2.2计算激光点P_j与激光点P_j+1的高程差，若高程差大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.2；

步骤2.2.3计算激光点P_j与激光点P_j-1的高程差，若高程差均大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.3；

步骤2.3设定回波强度差阈值，分析回波强度值

2.3.1计算激光点P_j与激光点P_j-1的回波强度差、激光点P_j与激光点P_j+1的回波强度差：若上述两个回波强度差都大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个回波强度差都小于回波强度差阈值，进行步骤2.4.1；若只有激光点P_j与激光点P_j+1的回波强度差小于回波强度差阈值，则进行步骤2.4.2；若只有激光点P_j与激光点P_j-1的回波强度差小于回波强度差阈值，则进行步骤2.4.3；

2.3.2计算激光点P_j与激光点P_j+1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.2；

2.3.3计算激光点P_j与激光点P_j-1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.3；

步骤2.4趋势角度计算

2.4.1采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\β}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-Z_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{(\mathrm{pj}-1)})}^2}}$

2.4.2采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\β}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-Z_{\left(\mathrm{Pj}\right)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{\left(\mathrm{Pj}\right)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{\left(\mathrm{pj}\right)})}^2}}$

2.4.3采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\β}}_j=\arctan \frac{\sqrt{{(Z_{\left(\mathrm{Pj}\right)}-Z_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{\left(\mathrm{Pj}\right)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{\left(\mathrm{Pj}\right)}-Y_{(\mathrm{pj}-1)})}^2}}$

式中，

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；；

Z_(Pj)为激光点P_j高程值；

Z_(Pj+1)为激光点P_j+1的高程值；

Z_(Pj-1)为激光点P_j-1的高程值；

步骤2.5趋势角数值正负判定

采用以下判据判定激光点P_j的趋势角β_j正负：

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

其中，α_(Pj-1)、α_(Pj)、α_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的扫描角，Z_(Pj-1)、Z_(Pj)、Z_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的高程值。

步骤3中，采用下式对初始机载LiDAR数据逐点计算回波强度最终校正值：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)·cos(|α_j|+β_j)

式中，I_{j初步校正值}为激光点P_j回波强度最终校正值。

作为优选方案：

步骤2.1中，所述平面距离阈值依据机载LiDAR数据的平均点密度取值，若激光点密度为n个/米²，则平面点间距为1/n米，距离阈值≧1/n米；

步骤2.2中，所述高程差阈值采用下式计算：

h_{高程差阈值}≧平面距离阈值/tan(90°-γ)

式中，γ为地物最大倾角；

步骤2.3中，回波强度差阈值依据地物类型设定，取1至15。

本发明的有益效果为：

(1)本发明的一种非植被覆盖区机载LiDAR数据回波强度值校正方法，基于机载LiDAR数据自身记录的信息，充分考虑了激光传输距离、扫描角、地物趋势角的影响，结合航高即可实现非植被类地物激光点的回波强度值的校正处理，可以有效的去除距离、反射角、地物倾角等因素造成的同类地物回波强度值差异，尤其是地物自身倾角带来的影响，从而提高激光回波强度值的准确性和利用价值，在地物分类及目标探测领域发挥更大作用；

(2)现有技术中采用采用传统的雷达方程式进行校正，由于参数难以获取导致校正难以实现，而国外已有的基于激光传输距离与扫描角度的归一化校正方法，则没有考虑到地物倾角对激光回波强度的影响。本发明提出的方法可以快速实现LiDAR数据回波强度值的校正处理，有效去除地物倾斜等因素的影响，校正后同类地物的回波强度值趋于一致，能够有效的用于进行地物分类、目标识别等工作；

(3)充分结合机载LiDAR数据的获取方式，考虑到在激光摆扫方向上的地物倾角影响大，与摆扫方向垂直的地物倾角可以忽略，因此选择P_j及相邻的P_j-1和P_j+1进行运算，既可以提高校正效果，又可以减少运算量，

(4)步骤2.1中，首先进行平面距离值分析，然后通过平面距离值的判定有效避免存储顺序相邻、实际空间距离不相邻的激光点之间的运算错误，如P_j和P_j+1分别位于航带的最边缘两侧，再如P_j和P_j+1之间存在水体(1.06um激光基本无回波信号)等地物，导致P_j和P_j+1相隔距离远且分属不同地物；

(5)步骤2.2可有效避免地物高程突变带来的计算误差，如P_j位于楼房顶部边缘，P_j+1紧邻地面，此时高差为楼房高度，且分属不同地物，若直接用于趋势角计算无疑是错误的，通过高程阈值判定，有效去除这一错误；

(6)步骤2.3中的回波强度阈值判定是在步骤1的初次校正后的数据基础上进行的，经过步骤1的初次校正后可以在一定程度上有效去除传输距离差异带来的影响，相同地物的回波强度值比较接近，因此步骤2.3的回波强度阈值判定可以有效去除部分小型地物的干扰，如房顶的太阳能板等小型设施；

(7)步骤2.4中充分利用激光扫描角方向与趋势角方向之间的关系，充分利用激光点的扫描角与高程值的变化相结合，制定了判定规则，确定地物趋势角与激光扫描角之间的加减关系，从计算出正确的反射角数值。

附图说明

图1为机载LiDAR的激光脉冲扫描原理示意图；

图2为依据原始回波强度值显示的实验区激光点云数据灰度图像；

图3为本发明的一种非植被覆盖区机载LiDAR数据回波强度值校正方法流程图；

图4(a)、(b)依次为三角形、拱形倾斜房顶激光点回波强度值分布示意图；

图5为激光扫描点平面分布示意图；

图6(a)、(b)依次为激光点房顶分布左视图、右视图；

图7为反射角计算(扫描角与地物趋势角加减关系)原理示意图；

图8为依据校正后回波强度值显示的实验区激光点云数据灰度图像；

图9(a)、(b)依次为经过初次校正后的三角形、拱形房顶激光点回波强度值分布显示图。

具体实施方式

下面结合附图和实施例对本发明的一种非植被覆盖区机载LiDAR数据回波强度值校正方法进行详细说明。

机载激光雷达的扫描作业过程如图1所示。机载LiDAR采用摆扫方式作业，当激光光斑扫描的地物水平时，如图1右侧所示，则反射角θ等于扫描角α；但当地物在扫描方向上具有一定的趋势角β，如图1左侧，则反射角θ等于扫描角α减去趋势角β。如背景技术中所述，传统的归一化校正方法，都是假定地面为朗伯体，且激光反射角等于入射角(即扫描角)，忽略地物倾角对激光脉冲回波强度的影响，没有提出计算地物趋势角的方法。本发明的主要创新点正在于此，通过实验分析倾斜地物的趋势角对激光回波强度的影响不可忽略，并提出一种计算地物趋势角的回波强度值校正方法。

实验数据基于ALTM Gemini机载激光雷达系统获得，该设备激光波长1064nm，波束角0.3mrad，具有4次回波记录能力。实验区位于石家庄机场区域，采用运五飞机为搭载平台，飞行海拔高度1100米，激光脉冲频率100kHz，共采集两条航线，1号航线摆扫角度为±25°、摆扫频率35Hz，平均点密度2个/米²，2号航线为±15°、摆扫频率50Hz，平均点密度3个/米²。图2所示为所采集的两条航线的LiDAR数据依据原始回波强度值显示的效果图。如前所述，由于受到激光传输距离、扫描角、地物反射等因素影响，可以看出图2中航带两侧边缘地物的回波强度值偏低，两条航带重叠区域的道路等地物回波强度值也不一致。

如图3所示，本发明的方法包括以下步骤：

步骤1机载激光雷达回波强度值的初次校正

设定机载LiDAR数据共有n个激光点数据，每个激光点的记录信息包括：激光点投影坐标(X，Y)、激光点高程Z、回波强度、扫描角等参数。

对n个激光点的高程值进行统计，确定高程最大值和高程最小值，然后依据下式计算高程平均值：

Z_平均＝(Z_最高+Z_最低)/2

式中，

Z_最高为高程最大值；

Z_最低为高程最小值；

Z_平均为高程平均值。

采用下式进行回波强度初步校正：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)，其中，R_j＝H-Z_j，R_A＝H-Z_平均。

式中，

H为机载LiDAR数据获取时的飞行航高；

Z_j为激光点P_j高程值，j＝1,…,n；

R_j为激光点P_j对地传输距离；

R_A为整个测区的n个激光点平均对地传输距离；

I_j原始值为激光点P_j回波强度原始值；

I_{j初步校正值}为激光点P_j回波强度初步校正值。

这一步骤主要是参考已有方法将激光传输距离对回波强度值的影响均一化，从而为下一步的回波强度值差值判定提供较好的数据基础。图4所示为三角形房顶(图4a)、拱形房顶(图4b)初步校正后得回波强度值分布显示，图中可以看出由于房顶两侧的反射角不同，造成回波强度值存在明显差异，这种现象经初步校正后并没有明显改观。

步骤2确定激光点趋势角

以步骤1校正后的机载LiDAR为基础，逐点计算每个激光点的趋势角。

步骤2.1分析平面距离

采用下式计算激光点P_j与临近激光点P_j-1、P_j+1的平面距离：

$L_j=\sqrt{{(X_{\left(P_j\right)}-X_{\left(P_{j-1}\right)})}^2+{(Y_{\left(P_j\right)}-Y_{\left(P_{j-1}\right)})}^2}$

$R_j=\sqrt{{(X_{\left(P_{j+1}\right)}-X_{\left(P_j\right)})}^2+{(Y_{\left(P_{j+1}\right)}-Y_{\left(P_j\right)})}^2}$

式中，

L_j为激光点P_j与相临近的P_j-1的平面距离；

R_j为激光点P_j与相临近的P_j+1的平面距离；

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；

设定平面距离阈值，然后将L_j、R_j与平面距离阈值进行比较：

若L_j、R_j均小于平面距离阈值，进行步骤2.2.1；

若只有R_j小于平面距离阈值，进行步骤2.2.2；

若只有L_j小于平面距离阈值，进行步骤2.2.3；

若L_j、R_j均大于平面距离阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；

本实施例中，两条航带数据的平均激光点密度分别为2个/米²和3个/米²，平均点间距为0.5米和0.3米，平面距离阈值大于等于平均点间距，取1米。

步骤2.2设定高差阈值，分析高程值

步骤2.2.1若激光点P_j的高程值同时大于其临近激光点P_j-1和P_j+1的高程值，则将激光点P_j的趋势角β_j取0，直接进行步骤3。否则，分别计算激光点P_j与激光点P_j-1的高程差、激光点P_j与激光点P_j+1的高程差，若上述两个高程差均大于高程差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个高程差均小于高程差阈值，则进行步骤2.3.1；若只有激光点P_j与激光点P_j+1的高程差小于高程差阈值，则进行步骤2.3.2；若只有激光点P_j与激光点P_j-1的高程差小于高程差阈值，则进行步骤2.3.3。

步骤2.2.2计算激光点P_j与激光点P_j+1的高程差，若高程差大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.2。

步骤2.2.3计算激光点P_j与激光点P_j-1的高程差，若高程差大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.3。

本实施例中，所述高程差阈值采用下式计算：

h_{高程差阈值}≧平面距离阈值/tan(90°-γ)

式中，

γ为地物最大倾角，其取值可以根据测区地物实际状况设置，一般不超过60°。

本实施例中，平面距离阈值为2.1步骤中所设置1米，经计算h_{高程差阈值}≧0.58米，最终取值0.8米。

步骤2.3设定回波强度差阈值，分析回波强度值

2.3.1计算激光点P_j与激光点P_j-1的回波强度差、激光点P_j与激光点P_j+1的回波强度差：若上述两个回波强度差均大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个回波强度差都小于回波强度差阈值，进行步骤2.4.1；若只有激光点P_j与激光点P_j+1的回波强度差小于回波强度差阈值，则进行步骤2.4.2；若只有激光点P_j与激光点P_j-1的回波强度差小于回波强度差阈值，则进行步骤2.4.3；

2.3.2计算激光点P_j与激光点P_j+1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.2；

2.3.3计算激光点P_j与激光点P_j-1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.3；

所述回波强度差阈值根据测区地物复杂程度具体设置，本实施例中，回波强度差阈值可以取1-15，优选值为15。

步骤2.4趋势角数值计算

2.4.1采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-Z_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{(\mathrm{pj}-1)})}^2}}$

2.4.2采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-Z_{\left(\mathrm{Pj}\right)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{\left(\mathrm{Pj}\right)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{\left(\mathrm{pj}\right)})}^2}}$

2.4.3采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{\left(\mathrm{Pj}\right)}-Z_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{\left(\mathrm{Pj}\right)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{\left(\mathrm{Pj}\right)}-Y_{(\mathrm{pj}-1)})}^2}}$

式中，

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；；

Z_(Pj)为激光点P_j高程值；

Z_(Pj+1)为激光点P_j+1的高程值；

Z_(Pj-1)为激光点P_j-1的高程值；

步骤2.5趋势角数值正负判定

机载LiDAR扫描激光点的分布如图5所示。由于激光雷达采用摆扫方式发射激光脉冲，加之飞机处于飞行运动中，因此激光点的分布呈图5所示呈“之”字形。因此结合图1分析，与扫描方向平行的地物趋势角影响较大，与扫描方向垂直的倾斜地物影响可以忽略，选取其中的3个激光点1、2、3。

如果其扫描顺序为从1至3，扫描角度值依次增大，地物倾斜如图6(a)、(b)所示包括两种情况，若3个激光点高程依次增加(图6a)，则如图7中激光束2所示情况，反射角等于扫描角绝对值减去趋势角，则趋势角为负，若3个激光点高程依次减小(图6b)，则如图7中激光束1所示情况，反射角等于扫描角绝对值加上趋势角，则趋势角为正；

如果其扫描顺序为从1至3，扫描角度值依次减小，地物倾斜如图6(a)、(b)所示包括两种情况，若3个激光点高程依次增加(图6a)，则如图7中激光束1所示情况，反射角等于扫描角绝对值加上趋势角，则趋势角为正，若3个激光点高程依次减小(图6b)，则如图7中激光束2所示情况，反射角等于扫描角绝对值减去趋势角，则趋势角为负。

因此，激光点P_j的趋势角β_j正负判定规则如下：

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

其中，α_(Pj-1)、α_(Pj)、α_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的扫描角，由于部分机载激光雷达系统通常以扫描器中轴法线为参考(图1、图7所示)将中轴法线两边扫描角以正、负数值表示加以区别，因此须转换为绝对值进行运算，Z_(Pj-1)、Z_(Pj)、Z_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的高程值。

步骤3回波强度最终校正

采用下式对原始机载LiDAR数据逐点计算回波强度最终校正值：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)·cos(|α_j|+β_j)

式中，I_{j初步校正值}为激光点P_j回波强度最终校正值。

图8为校正后的航带数据依据回波强度值显示的灰度图，对比图2可以看出整个跑道及周边混凝土路面回波强度值比较均一，总体效果明显改善。

图9是加入趋势角计算校正后的三角形、拱形房顶回波强度值显示效果，对比图4其显示效果得到明显改善。进一步分别对三角形房顶82个激光点、拱形房顶3470个激光点的标准差进行计算分析为：

计算结果证明，经过趋势角度计算校正后，倾斜地物回波强度值的均方差更小，能够更好的反映出倾斜地物的介质属性，从而提高其分类识别的精确度。

Claims (5)

1.一种非植被覆盖区机载LiDAR数据回波强度值校正方法，其特征在于：包括以下步骤：

步骤1回波强度初步校正；

步骤2确定激光点趋势角；

步骤3回波强度最终校正。

2.根据权利要求1所述的一种非植被覆盖区机载LiDAR数据回波强度值校正方法，其特征在于：步骤1中，首先采用依据下式计算高程平均值：

Z_平均＝(Z_最高+Z_最低)/2

式中，

Z_最高为高程最大值；

Z_最低为高程最小值；

Z_平均为高程平均值；

然后依据下式进行回波强度初步校正：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)，其中，R_j＝H-Z_j，R_A＝H-Z_平均；

式中，

H为机载LiDAR数据获取时的飞行航高；

Z_j为激光点P_j高程值，j＝1,…,n；

R_j为激光点P_j对地传输距离；

R_A为整个测区的n个激光点平均对地传输距离；

I_j原始值为激光点P_j回波强度原始值；

I_{j初步校正值}为激光点P_j回波强度初步校正值。

3.根据权利要求2所述的一种非植被覆盖区机载LiDAR数据回波强度值校正方法，其特征在于：步骤2包括以下步骤：

步骤2.1分析平面距离

采用下式计算临近激光点平面距离：

$L_j=\sqrt{{(X_{\left(P_j\right)}-X_{\left(P_{j-1}\right)})}^2+{(Y_{\left(P_j\right)}-Y_{\left(P_{j-1}\right)})}^2}$

$R_j=\sqrt{{(X_{\left(P_{j+1}\right)}-X_{\left(P_j\right)})}^2+{(Y_{\left(P_{j+1}\right)}-Y_{\left(P_j\right)})}^2}$

式中，

L_j为激光点P_j与相临近的P_j-1的平面距离；

R_j为激光点P_j与相临近的P_j+1的平面距离；

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；

设定平面距离阈值，然后将L_j、R_j与平面距离阈值进行比较：

若L_j、R_j均小于平面距离阈值，进行步骤2.2.1；

若只有R_j小于平面距离阈值，进行步骤2.2.2；

若只有L_j小于平面距离阈值，进行步骤2.2.3；

若L_j、R_j均大于平面距离阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；

步骤2.2设定高差阈值，分析高程值

步骤2.2.1若激光点P_j的高程值同时大于其临近激光点P_j-1和P_j+1的高程值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，分别计算激光点P_j与激光点P_j-1的高程差、激光点P_j与激光点P_j+1的高程差，若上述两个高程差均大于高程差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个高程差均小于高程差阈值，则进行步骤2.3.1；若只有激光点P_j与激光点P_j+1的高程差小于高程差阈值，则进行步骤2.3.2；若只有激光点P_j与激光点P_j-1的高程差小于高程差阈值，则进行步骤2.3.3；

步骤2.2.2计算激光点P_j与激光点P_j+1的高程差，若高程差大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.2；

步骤2.2.3计算激光点P_j与激光点P_j-1的高程差，若高程差均大于高程差阈值，则将激光点P_j的趋势角β_j取0直接进行步骤3，否则进行步骤2.3.3；

步骤2.3设定回波强度差阈值，分析回波强度值

2.3.1计算激光点P_j与激光点P_j-1的回波强度差、激光点P_j与激光点P_j+1的回波强度差：若上述两个回波强度差都大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；若上述两个回波强度差都小于回波强度差阈值，进行步骤2.4.1；若只有激光点P_j与激光点P_j+1的回波强度差小于回波强度差阈值，则进行步骤2.4.2；若只有激光点P_j与激光点P_j-1的回波强度差小于回波强度差阈值，则进行步骤2.4.3；

2.3.2计算激光点P_j与激光点P_j+1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.2；

2.3.3计算激光点P_j与激光点P_j-1的回波强度差，若回波强度差大于回波强度差阈值，则将激光点P_j的趋势角β_j取0，直接进行步骤3；否则，进行步骤2.4.3；

步骤2.4趋势角度计算

2.4.1采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-X_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{(\mathrm{Pj}-1)})}^2}}$

2.4.2采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{(\mathrm{Pj}+1)}-X_{\left(\mathrm{Pj}\right)})}^2}}{\sqrt{{(X_{(\mathrm{Pj}+1)}-X_{\left(\mathrm{Pj}\right)})}^2+{(Y_{(\mathrm{Pj}+1)}-Y_{\left(\mathrm{Pj}\right)})}^2}}$

2.4.3采用下式计算激光点P_j的趋势角β_j，并进行步骤2.5：

${\mathrm{\&beta;}}_j=\arctan \frac{\sqrt{{(Z_{\left(\mathrm{Pj}\right)}-X_{(\mathrm{Pj}-1)})}^2}}{\sqrt{{(X_{\left(\mathrm{Pj}\right)}-X_{(\mathrm{Pj}-1)})}^2+{(Y_{\left(\mathrm{Pj}\right)}-Y_{(\mathrm{Pj}-1)})}^2}}$

式中，

(X_(Pj)、Y_(Pj))为激光点P_j的投影坐标；

(X_(Pj+1)、Y_(Pj+1))为激光点P_j+1的投影坐标；

(X_(Pj-1)、Y_(Pj-1))为激光点P_j-1的投影坐标；；

Z_(Pj)为激光点P_j高程值；

Z_(Pj+1)为激光点P_j+1的高程值；

Z_(Pj-1)为激光点P_j-1的高程值；

步骤2.5趋势角数值正负判定

采用以下判据判定激光点P_j的趋势角β_j正负：

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|<|α_(Pj)|<|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)<Z_(Pj)<Z_(Pj+1)，则β_j取负值；

若|α_(Pj-1)|>|α_(Pj)|>|α_(Pj+1)|、Z_(Pj-1)>Z_(Pj)>Z_(Pj+1)，则β_j取正值；

其中，α_(Pj-1)、α_(Pj)、α_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的扫描角，Z_(Pj-1)、Z_(Pj)、Z_(Pj+1)依次表示激光点P_j-1、P_j、P_j+1的高程值。

4.根据权利要求3所述的一种非植被覆盖区机载LiDAR数据回波强度值校正方法，其特征在于：步骤3中，采用下式对初始机载LiDAR数据逐点计算回波强度最终校正值：

I_{j初步校正值}＝I_j原始值·(R_j ²/R_A ²)·cos(|α_j|+β_j)

式中，I_{j初步校正值}为激光点P_j回波强度最终校正值。

5.根据权利要求2或3或4所述的一种非植被覆盖区机载LiDAR数据回波强度值校正方法，其特征在于：

步骤2.1中，所述平面距离阈值依据机载LiDAR数据的平均点密度取值，若激光点密度为n个/米²，则平面点间距为1/n米，距离阈值≧1/n米；

步骤2.2中，所述高程差阈值采用下式计算：

h_{高程差阈值}≧平面距离阈值/tan(90°-γ)

式中，γ为地物最大倾角；

步骤2.3中，回波强度差阈值依据地物类型设定，取1至15。