CN117347971B - 一种激光雷达建筑物边界估测方法及计算机可读介质 - Google Patents
一种激光雷达建筑物边界估测方法及计算机可读介质 Download PDFInfo
- Publication number
- CN117347971B CN117347971B CN202311146876.2A CN202311146876A CN117347971B CN 117347971 B CN117347971 B CN 117347971B CN 202311146876 A CN202311146876 A CN 202311146876A CN 117347971 B CN117347971 B CN 117347971B
- Authority
- CN
- China
- Prior art keywords
- building
- satellite
- laser
- laser radar
- borne
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000005259 measurement Methods 0.000 claims abstract description 13
- 230000006870 function Effects 0.000 claims description 56
- 230000000295 complement effect Effects 0.000 claims description 16
- 230000003287 optical effect Effects 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 230000003595 spectral effect Effects 0.000 claims description 9
- 238000002310 reflectometry Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 7
- 238000001228 spectrum Methods 0.000 claims description 6
- 238000013507 mapping Methods 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 3
- 238000005316 response function Methods 0.000 claims description 3
- 230000009131 signaling function Effects 0.000 claims description 3
- GHMWZRWCBLXYBX-UHFFFAOYSA-M sodium;4-chlorobenzoate Chemical compound [Na+].[O-]C(=O)C1=CC=C(Cl)C=C1 GHMWZRWCBLXYBX-UHFFFAOYSA-M 0.000 claims description 3
- 230000002123 temporal effect Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 11
- 238000010586 diagram Methods 0.000 description 4
- 238000000605 extraction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Abstract
本发明公开一种激光雷达建筑物边界估测方法及计算机可读介质，属于激光遥感技术领域。本发明首先根据扩展光斑效应的基本原理，建立了建筑物目标的激光回波信号理论模型；然后通过统计飞越城市或其他包含建筑物目标的区域时的星载单光子激光雷达数据点云密度，计算建筑物目标的实测点云密度曲线；最后结合实测点云密度曲线与回波信号理论模型，基于最小二乘方法在复频域拟合估计建筑物的边界位置，提取对应的光子点云。本发明能准确确定星载单光子激光雷达数据中建筑物边界位置，有效消除星载激光雷达的扩展光斑效应，所提取的建筑物边界点云可作为与其他遥感数据源配准和融合的基本特征点，在城市三维遥感领域具有很大的应用价值和潜力。
Description
技术领域
本发明属于激光遥感技术领域，尤其涉及一种激光雷达建筑物边界估测方法及计算机可读介质。
背景技术
星载激光雷达作为一种可用于对地观测(测高)的主动遥感设备，是通过卫星平台获取地表绝对地理坐标精度最高的遥感设备，具有快速获取全球范围地表高精度高程控制点的能力。新一代光子计数(单光子)体制的星载激光雷达具有微脉冲、高重频、多波束等显著优点，已经成为国内外星载激光雷达(激光测高)技术发展的主要趋势。目前，美国激光测高卫星ICESat-2(Ice,Cloud and land Elevation Satellite-2)上搭载的ATLAS(Advanced Topographic Laser Altimeter System)是全球唯一在轨运行的星载单光子激光雷达。得益于高达10kHz的激光重复频率，ICESat-2/ATLAS沿卫星飞行方向的相邻激光测量点间距仅为0.7m，使得ICESat-2/ATLAS可以直接获取近乎连续的沿轨地表高程轮廓，这是以往的星载测高手段(例如微波高度计、线性体制激光测高仪)所无法实现的。
目前，ICESat-2卫星获取的光子点云数据已被广泛应用于地形测绘、林业调查、海洋环境和极地环境监测等研究领域。除上述应用外，ICESat-2在城市遥感方面同样具有广阔的应用前景。从ICESat-2光子点云数据中提取的建筑物边界，即获取地面激光高程轮廓与建筑物顶面边界的交点，可作为与其他遥感数据源(例如光学、SAR图像)配准和融合的基本特征点。
准确提取ICESat-2光子点云数据中的建筑物边界面临的主要问题是：由于卫星飞行高度约为500km，ICESat-2星载单光子激光雷达的地面光斑直径达到了约11m，远远超过了沿轨方向上相邻两光斑中心0.7m的间距，而ICESat-2接收到的信号实际上是光斑范围内的所有地表目标对发射激光脉冲的调制结果，即光斑内的能量分布与目标高程分布、反射率分布等地面目标响应之间的空间卷积，可以称为“扩展光斑效应”。该效应会造成ICESat-2光子点云数据中相邻地表目标之间的边界模糊现象。在城市地区，扩展光斑效应导致ICESat-2光子点云数据中的地面光子点云与建筑物光子点云在建筑物屋顶边界处产生了重叠现象，所呈现建筑物目标的测量宽度也明显大于其实际宽度，难以确定建筑物边界的准确位置。
实际上，“扩展光斑效应”在所有星载激光雷达都存在，本发明中以ICESat-2/ATLAS的光子点云数据为例，阐述针对建筑物目标，有效剔除该效应，从星载单光子激光雷达的光子点云数据中准确提取建筑物边界的方法。
发明内容
为了解决上述技术问题，本发明提出了一种激光雷达建筑物边界估测方法及计算机可读介质。
本发明方法的技术方案为一种激光雷达建筑物边界估测方法，其特征在于，包括：
建立星载激光雷达接收的激光回波信号功率分布模型，建立描述建筑物几何和反射特征的目标模型，结合互补误差函数构建建筑物目标的激光回波信号理论模型；
计算每个网格单元的点云密度，结合选定区域每个网格单元的点云密度计算激光雷达飞越建筑物时光子点云数据的实际测量值；
在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云。
本发明的具体步骤如下：
步骤1：建立星载激光雷达接收的激光回波信号功率分布模型，建立描述建筑物几何和反射特征的目标模型，结合互补误差函数构建建筑物目标的激光回波信号理论模型；
步骤2：将星载单光子激光雷达飞越城市或其他包含建筑物目标的区域时采集的坐标为沿轨位置和高程的光子点云数据，按照高程间隔、沿轨距离间隔划分为二维网格，统计落在各网格单元内的光子点云个数从而得到每个网格单元的点云密度，通过目标识别算法选取建筑物目标的光子点云沿轨、高程所在的选定区域，对选定区域内网格的每个单元的点云密度在高程方向上进行求和，计算激光雷达飞越建筑物时光子点云数据的实际测量值；
步骤3：在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云；
作为优选，步骤1所述建立星载激光雷达接收的激光回波信号功率分布模型，具体定义如下：
激光接收信号函数表示为地面目标响应函数与发射激光脉冲函数的卷积；
发射激光脉冲在空间上服从二维高斯分布，时间分布近似为高斯函数；
所述星载激光雷达接收的激光回波信号功率分布模型，具体如下：
其中，Pr(x,y,t)表示下文构建的星载激光雷达平面坐标系下的地表光斑中心的平面坐标(x,y)处时间t的接收光学系统接收得到的激光功率，表示卷积运算，Pt(t)是发射激光脉冲的功率分布，表示为：
其中，t表示激光脉冲的飞行时间，E0表示发射激光脉冲的总能量，σt表示发射脉冲宽度；
将星载激光雷达沿卫星飞行方向定义为沿轨方向，将星载激光雷达垂直于卫星飞行方向定义为垂轨方向，通过沿轨方向、垂轨方向构建星载激光雷达平面坐标系下；
其中，(x,y)表示星载激光雷达平面坐标系下的地表光斑中心的平面坐标，x表示星载激光雷达平面坐标系下的地表光斑中心的沿轨方向的坐标，y表示星载激光雷达平面坐标系下的地表光斑中心的垂轨方向的坐标；
(u,v)表示星载激光雷达平面坐标系下的光斑内微面元的平面坐标，u表示星载激光雷达平面坐标系下的光斑内微面元的沿轨方向的坐标，v表示星载激光雷达平面坐标系下的光斑内微面元的垂轨方向的坐标；h表示垂直于星载激光雷达平面坐标系所在的平面并指向天顶的垂直高程，κ0为总衰减系数，表示地面的法向量，/>表示激光指向向量，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，ρ(u,v,h)表示(u,v,h)处的目标反射率，l表示卫星上的激光参考点与地表光斑内微面元的位置(u,v,h)的几何距离，如下式所示：
其中，c为大气中的光速，(xs,ys,hs)表示星载激光雷达平面坐标系下星载激光雷达的激光参考点坐标，xs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的沿轨方向的坐标，ys表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂轨方向的坐标，hs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂直于星载激光雷达平面坐标系的平面并指向天顶的垂直高程；
步骤1所述建立描述建筑物几何和反射特征的目标模型，具体定义如下：
假设建筑物目标表面具有一致的反射特征，对于近天底入射的激光脉冲，建筑物的反射表面为建筑物顶面，即假设接收光学系统接收的来自建筑物顶面上每一处微面的反射率相同；
所述建筑物顶面几何和反射特征的目标模型，具体定义如下：
γ(u,v,t)＝γ0δ(h-Ht(u,v))(u,v)∈Duv
其中，δ(·)表示狄拉克函数，Duv表示建筑物在地面上的投影范围，γ(u,v,t)对应为S1.1的星载激光雷达接收的激光回波信号功率分布模型中的
对于被激光光斑照射的一个建筑物目标，激光指向矢量建筑物屋顶表面法向量/>几何距离l、总衰减系数κ0均为常数，使用衰减及反射常系数γ0综合表示建筑物顶面的反射特征。为了描述建筑物目标的空间分布情况，使用三维曲面方程h-Ht(u,v)＝0表示建筑物顶面的几何形状；
假设与激光轨迹相交的两条建筑物边界线即Duv的边界分别为规则的直线段l1和规则的直线段l2，(x1,y1)、(x2,y2)分别为前述设定的坐标系下l1、l2与激光轨迹的交点，x1和x2分别表示前述设定的平面坐标系下l1、l2与激光轨迹的交点的沿轨方向的坐标，y1和y2分别表示星载激光雷达平面坐标系下l1、l2与激光轨迹的交点的垂轨方向的坐标；
在星载激光雷达平面坐标系中y1＝y2＝0，星载激光雷达平面坐标系下l1、l2的直线方程分别表示为：
y＝k1(x-x1)
y＝k2(x-x2)
其中，k1和k2分别表示前述设定的平面坐标系下l1、l2的直线方程的斜率；
步骤1所述结合互补误差函数构建建筑物目标的激光回波信号理论模型，具体定义如下：
在时间轴t上积分得到建筑物目标反射到星载激光雷达接收光学系统的激光能量分布er(x,y)，表示为：
卫星沿固定轨道飞行，光斑中心坐标(x,y)的两个元素之间存在映射关系，er(x,y)简化表示为er(x)；
通过引入互补误差函数erfc(·)代替垂轨方向的高斯函数积分，进一步得到建筑物目标的激光回波信号理论模型，具体如下：
其中，E0表示发射激光脉冲的总能量，x表示星载激光雷达平面坐标系下地表光斑中心的沿轨方向的坐标，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，表示卷积运算，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标；
作为优选，步骤2所述统计飞越城市或其他包含建筑物目标的区域时的星载单光子激光雷达数据点云密度，具体如下：
每个网格单元的点云密度计算公式如下：
其中，dmea(m,n)表示网格中第m行第n列的单元的点云密度，Lh表示网格划分的高程间隔，Lx表示网格划分的沿轨距离间隔，(m,n)表示网格中第m行第n列的单元，对应网格中心位置的高程和沿轨距离分别为mLh和nLx，Wh和Wx分别为统计点云密度时使用的窗口在高程和沿轨距离方向上的半宽，hpk和xpk分别表示整段数据中第k个点云所对应的高程和沿轨距离值，ε(·)表示单位阶跃函数，K表示整段数据的总光子点云数；
步骤2所述计算激光雷达飞越建筑物时光子点云数据的实际测量值，具体如下：
其中，Dtarg表示自动或人工选取的建筑物目标在沿轨、高程方向上所包含网格单元的集合；建筑物目标激光雷达沿轨实测的累计点云密度曲线smea(x)与建筑物目标激光雷达理论上沿轨接收的能量分布er(x)呈线性关系，即：
smea(x)＝C0·er(x)
其中，C0表示建筑物目标未知常系数，smea(x)表示计算得到的激光雷达飞越建筑物时光子点云数据的实际测量值；
作为优选，步骤3所述在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，具体如下：
将激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)与建筑物目标的激光回波信号理论模型er(x)在复频域中使用非线性最小二乘法进行拟合和参数估计；
激光回波理论模型er(x)经过离散傅里叶变换得到的频谱为Er(ω)，即其中/>表示离散傅里叶变换，具体公式如下：
其中，ω表示角频率，j表示复数虚部，E0表示发射激光脉冲的总能量，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标，复频域中基于最小二乘法得到的代价函数f(x1,x2,k1,k2,C)如下式所示：
其中，星载激光雷达平面坐标系下建筑物两个边界在星载激光雷达沿轨方向的准确位置x1和x2、建筑物两个边界直线方程的斜率常数k1和k2以及建筑物目标常未知系数C为待估计的未知参数，发射激光脉冲的总能量E0、衰减及反射常系数γ0与建筑物目标常未知系数C0合并为未知常系数C来估计，即C＝C0E0γ0；
Smea(ω)是激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)经过离散傅里叶变换得到的频谱函数，即ωi表示序列Smea(ω)中第i个频谱分量所对应的角频率，Smea(ωi)表示序列Smea(ω)中第i个频谱分量，I表示序列Smea(ω)的频谱分量总数的集合。
步骤3所述通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云，具体如下：
使用Levenberg-Marquardt法，估计复频域中基于最小二乘法得到的代价函数的极小值点，准确得到建筑物边界在星载激光雷达沿轨方向上的位置参数x1*、x2*的最佳估计结果。光子点云数据包含每个光子点云的沿轨位置和高程，代入位置参数x1*、x2*到光子点云数据中，可以得到对应的建筑高度h1*和h2*。
本发明还提供了一种计算机可读介质，所述计算机可读介质存储电子设备执行的计算机程序，当所述计算机程序在电子设备上运行时，执行所述激光雷达建筑物边界估测方法的步骤。
与现有技术相比，本发明所提供的方法能准确确定星载单光子激光雷达数据中建筑物边界位置，有效消除星载激光雷达的扩展光斑效应，所提取的建筑物边界点云可作为与其他遥感数据源配准和融合的基本特征点，在城市三维遥感领域具有很大的应用价值和潜力。同时，本发明也为有效处理和消除星载激光雷达数据中的扩展光斑效应问题提供了一种新思路。
附图说明
图1：本发明实施例的方法流程图；
图2：本发明实施例的激光地面光斑轨迹及建筑物边界示意图；
图3：本发明实施例的星载单光子激光雷达获取的光子点云实际测量值示意图；
图4：本发明实施例的星载单光子激光雷达数据的建筑物边界提取结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。
实施例以美国ICESat-2卫星搭载的单光子激光雷达ATLAS所获取的光子点云数据作为输入样例。ICESat-2/ATLAS激光脉冲重复频率为10kHz，地面光斑直径约11m，以沿轨约0.7m的间隔获取地表及其覆盖物的光子点云。ICESat-2星载单光子激光雷达ATL03标准数据产品提供了沿卫星飞行方向所获取的光子点云数据，如图3(a)中黑色散点所示，横轴为沿轨距离，单位为米，纵轴为高程，单位为米；同时，ATL03标准数据产品也提供了每个光子点云对应的获取时刻、经纬度坐标(基于WGS84椭球)。具体实施样例选择为ICESat-2/ATLAS在新西兰惠灵顿大区哈特市一条激光轨迹所包含的光子点云数据，该轨迹飞越一幢建筑物，该建筑物具有较规则的边界和相对一致的表面反射率。实施例中所使用参数取值为：Lh＝0.15m，Lx＝0.2m，Wh＝0.3m，Wx＝6m；由于光斑直径为11m，因此激光束截面参数σB＝2.75m。
下面结合图1至图4介绍本发明实施例的具体实施方式为一种激光雷达数据中建筑物边界估测方法，具体如下：
如图1所示，为本发明实施例的方法流程图。
步骤1：建立星载激光雷达接收的激光回波信号功率分布模型，建立描述建筑物几何和反射特征的目标模型，结合互补误差函数构建建筑物目标的激光回波信号理论模型；
步骤1所述建立星载激光雷达接收的激光回波信号功率分布模型，具体定义如下：
激光接收信号函数表示为地面目标响应函数与发射激光脉冲函数的卷积；
发射激光脉冲在空间上服从二维高斯分布，时间分布近似为高斯函数；
所述星载激光雷达接收的激光回波信号功率分布模型，具体如下：
其中，Pr(x,y,t)表示下文构建的星载激光雷达平面坐标系下的地表光斑中心的平面坐标(x,y)处时间t的接收光学系统接收得到的激光功率，表示卷积运算，Pt(t)是发射激光脉冲的功率分布，表示为：
其中，t表示激光脉冲的飞行时间，E0表示发射激光脉冲的总能量，σt表示发射脉冲宽度；
将星载激光雷达沿卫星飞行方向定义为沿轨方向，将星载激光雷达垂直于卫星飞行方向定义为垂轨方向，通过沿轨方向、垂轨方向构建星载激光雷达平面坐标系下；
其中，(x,y)表示星载激光雷达平面坐标系下的地表光斑中心的平面坐标，x表示星载激光雷达平面坐标系下的地表光斑中心的沿轨方向的坐标，y表示星载激光雷达平面坐标系下的地表光斑中心的垂轨方向的坐标；
(u,v)表示星载激光雷达平面坐标系下的光斑内微面元的平面坐标，u表示星载激光雷达平面坐标系下的光斑内微面元的沿轨方向的坐标，v表示星载激光雷达平面坐标系下的光斑内微面元的垂轨方向的坐标；h表示垂直于星载激光雷达平面坐标系所在的平面并指向天顶的垂直高程，κ0为总衰减系数，表示地面的法向量，/>表示激光指向向量，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，ρ(u,v,h)表示(u,v,h)处的目标反射率，l表示卫星上的激光参考点与地表光斑内微面元的位置(u,v,h)的几何距离，如下式所示：
其中，c为大气中的光速，(xs,ys,hs)表示星载激光雷达平面坐标系下星载激光雷达的激光参考点坐标，xs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的沿轨方向的坐标，ys表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂轨方向的坐标，hs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂直于星载激光雷达平面坐标系的平面并指向天顶的垂直高程；
步骤1所述建立描述建筑物几何和反射特征的目标模型，具体定义如下：
假设建筑物目标表面具有一致的反射特征，对于近天底入射的激光脉冲，建筑物的反射表面为建筑物顶面，即假设接收光学系统接收的来自建筑物顶面上每一处微面的反射率相同；
所述建筑物顶面几何和反射特征的目标模型，具体定义如下：
γ(u,v,t)＝γ0δ(h-Ht(u,v))(u,v)∈Duv
其中，δ(·)表示狄拉克函数，Duv表示建筑物在地面上的投影范围，γ(u,v,t)对应为S1.1的星载激光雷达接收的激光回波信号功率分布模型中的
对于被激光光斑照射的一个建筑物目标，激光指向矢量建筑物屋顶表面法向量/>几何距离l、总衰减系数κ0均为常数，使用衰减及反射常系数γ0综合表示建筑物顶面的反射特征。为了描述建筑物目标的空间分布情况，使用三维曲面方程h-Ht(u,v)＝0表示建筑物顶面的几何形状；
如图2所示，假设与激光轨迹相交的两条建筑物边界线即Duv的边界分别为规则的直线段l1和规则的直线段l2，(x1,y1)、(x2,y2)分别为前述设定的坐标系下l1、l2与激光轨迹的交点，x1和x2分别表示前述设定的平面坐标系下l1、l2与激光轨迹的交点的沿轨方向的坐标，y1和y2分别表示星载激光雷达平面坐标系下l1、l2与激光轨迹的交点的垂轨方向的坐标；
在星载激光雷达平面坐标系中y1＝y2＝0，星载激光雷达平面坐标系下l1、l2的直线方程分别表示为：
y＝k1(x-x1)
y＝k2(x-x2)
其中，k1和k2分别表示前述设定的平面坐标系下l1、l2的直线方程的斜率；
步骤1所述结合互补误差函数构建建筑物目标的激光回波信号理论模型，具体定义如下：
在时间轴t上积分得到建筑物目标反射到星载激光雷达接收光学系统的激光能量分布er(x,y)，表示为：
卫星沿固定轨道飞行，光斑中心坐标(x,y)的两个元素之间存在映射关系，er(x,y)简化表示为er(x)；
通过引入互补误差函数erfc(·)代替垂轨方向的高斯函数积分，进一步得到建筑物目标的激光回波信号理论模型，具体如下：
其中，E0表示发射激光脉冲的总能量，x表示星载激光雷达平面坐标系下地表光斑中心的沿轨方向的坐标，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，表示卷积运算，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标；
步骤2：将星载单光子激光雷达飞越城市或其他包含建筑物目标的区域时采集的坐标为沿轨位置和高程的光子点云数据，按照高程间隔、沿轨距离间隔划分为二维网格，统计落在各网格单元内的光子点云个数从而得到每个网格单元的点云密度，通过目标识别算法选取建筑物目标的光子点云沿轨、高程所在的选定区域，对选定区域内网格的每个单元的点云密度在高程方向上进行求和，计算激光雷达飞越建筑物时光子点云数据的实际测量值；
步骤2所述统计飞越城市或其他包含建筑物目标的区域时的星载单光子激光雷达数据点云密度，具体如下：
将星载单光子激光雷达飞越城市或其他包含建筑物目标的区域时采集的坐标为沿轨位置和高程的光子点云数据，按照高程间隔Lh、沿轨距离间隔Lx划分为二维网格，统计落在各网格单元内的光子点云个数，从而得到每个网格单元的点云密度；
每个网格单元的点云密度计算公式如下：
其中，dmea(m,n)表示网格中第m行第n列的单元的点云密度，Lh表示网格划分的高程间隔，Lx表示网格划分的沿轨距离间隔，(m,n)表示网格中第m行第n列的单元，对应网格中心位置的高程和沿轨距离分别为mLh和nLx，Wh和Wx分别为统计点云密度时使用的窗口在高程和沿轨距离方向上的半宽，hpk和xpk分别表示整段数据中第k个点云所对应的高程和沿轨距离值，ε(·)表示单位阶跃函数，K表示整段数据的总光子点云数；
步骤2所述计算建筑物目标的实测点云密度曲线，具体如下：
通过常见的目标识别算法或直接使用人工框选方式，选取建筑物目标的光子点云沿轨、高程所在的选定区域，对选定区域内网格的点云密度dmea(m,n)在高程方向上进行求和，计算得出沿轨方向建筑物目标的高程累计点云密度曲线smea(x)：
其中，Dtarg表示自动或人工选取的建筑物目标在沿轨、高程方向上所包含网格单元的集合；建筑物目标激光雷达沿轨实测的累计点云密度曲线smea(x)与建筑物目标激光雷达理论上沿轨接收的能量分布er(x)呈线性关系，即：
smea(x)＝C0·er(x)
其中，C0表示建筑物目标未知常系数，smea(x)表示计算得到的激光雷达飞越建筑物时光子点云数据的实际测量值；
此步骤对应图3(a)到(b)的过程，即由离散光子点云按照网格统计得到点云密度dmea(m,n)，图3(a)是ICESat-2获取的光子点云数据，图3(b)是点云密度分布图，点云密度由式(7)计算得到，其中参数在实施例中选用Lh＝0.15m，Lx＝0.2m，Wh＝0.3m，Wx＝6m。图3(b)右边灰度条表示点云密度结果的数值范围，其单位为pts/m2，每平方米点个数。
步骤3：在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云；
步骤3所述在复频域构建实测点云密度曲线与回波信号理论模型之间的代价函数，具体如下：
将激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)与建筑物目标的激光回波信号理论模型er(x)在复频域中使用非线性最小二乘法进行拟合和参数估计；
激光回波理论模型er(x)经过离散傅里叶变换得到的频谱为Er(ω)，即其中/>表示离散傅里叶变换，具体公式如下：
其中，ω表示角频率，j表示复数虚部，E0表示发射激光脉冲的总能量，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标，复频域中基于最小二乘法得到的代价函数f(x1,x2,k1,k2,C)如下式所示：
其中，星载激光雷达平面坐标系下建筑物两个边界在星载激光雷达沿轨方向的准确位置x1和x2、建筑物两个边界直线方程的斜率常数k1和k2以及建筑物目标常未知系数C为待估计的未知参数，发射激光脉冲的总能量E0、衰减及反射常系数γ0与建筑物目标常未知系数C0合并为未知常系数C来估计，即C＝C0E0γ0；
Smea(ω)是激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)经过离散傅里叶变换得到的频谱函数，即ωi表示序列Smea(ω)中第i个频谱分量所对应的角频率，Smea(ωi)表示序列Smea(ω)中第i个频谱分量，I表示序列Smea(ω)的频谱分量总数的集合。
步骤3所述通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云，具体如下：
使用Levenberg-Marquardt法，估计复频域中基于最小二乘法得到的代价函数的极小值点，准确得到建筑物边界在星载激光雷达沿轨方向上的位置参数x1*、x2*的最佳估计结果。光子点云数据包含每个光子点云的沿轨位置和高程，代入位置参数x1*、x2*到光子点云数据中，可以得到对应的建筑高度h1*和h2*。
为了评估该发明方法的可行性和准确性，本实施例选择ICESat-2样例数据所在区域当地高精度机载激光雷达点云数据作为验证数据。本实施例中建筑物边界的提取结果，以目标建筑物在ICESat-2沿轨轨迹上的宽度(即两个边界点云之间距离|x2-x1|)作为评估标准，其与机载激光雷达点云数据计算得到的建筑物宽度对比，相对误差为3％，表明了本发明方法的可行性和准确性。
图4详细展示了实施例中采用本发明方法得到的一例建筑物边界提取结果，该结果与卫星光学影像数据、机载点云数据、ICESat-2星载单光子激光雷达ATL03数据产品中光子点云数据进行了对比，图4(a)为高分辨率卫星光学影像，右下角的箭头指向正北方向；图4(b)为机载高精度、高密度点云数据，图中的点云灰度与右侧灰度条对应表明了点云高程信息；图4(a)、图4(b)中的水平白色虚线表示ICESat-2飞越建筑物时的地面激光轨迹。图4(c)为ICESat-2获取的光子点云数据，浅色圆圈所圈的点是建筑物目标对应的光子点云(位于人工框选的虚线方框内)，图4(d)为图4(c)中光子点云计算的点云密度分布图，点云密度dmea(m,n)由S2.1的式(7)计算得到，右边灰度条标明了点云密度数值范围，其单位为pts/m2(每平方米点个数)。
图4(e)中不平滑曲线为建筑物目标的累计点云密度曲线smea(x)，由图4(d)中位于虚线方框内的点云密度dmea(m,n)代入S2.2的式(8)计算得到；经过最小二乘拟合，图4(e)不平滑曲线对应的频域最佳拟合曲线式(10)，再经过傅里叶逆变换后，得到图4(e)中的平滑曲线，其中拟合得到沿轨位置x1、x2在图4(e)中使用垂直方向虚线表示。图4(a)至图4(e)的垂直虚线表示使用本发明方法提取的建筑目标准确边界位置，图4(a)至图4(e)中的垂直实线表示ICESat-2实测的原始光子点云对应建筑物边界位置，通过与图4(a)卫星光学影像、图4(b)机载高精度高密度点云数据中的建筑物屋顶边界与ICESat-2沿轨轨迹的交点相对比，可以看出本发明方法的准确性更高。垂直实线对应的ICESat-2实测原始光子点云所获取的建筑物两侧边界，由于光斑扩展效应影响，显然远超出实际建筑物边界，难以直接确定建筑物边界的准确位置。
因此，本发明所提供的方法能准确确定星载单光子激光雷达数据中建筑物边界，有效消除星载激光雷达的扩展光斑效应，所提取的建筑物边界点云可以在未来城市区域遥感数据配准和融合中发挥重要作用。同时，本发明也为有效解决星载激光雷达数据中的扩展光斑效应问题提供了一种新思路。
本发明的具体实施例还提供了一种计算机可读介质。
所述计算机可读介质为服务器工作站；
所述服务器工作站存储电子设备执行的计算机程序，当所述计算机程序在电子设备上运行时，使得所述电子设备执行本发明实施例的激光雷达建筑物边界估测方法的步骤。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。
Claims (2)
1.一种激光雷达建筑物边界估测方法，其特征在于，包括以下步骤：
步骤1：建立星载激光雷达接收的激光回波信号功率分布模型，建立描述建筑物几何和反射特征的目标模型，结合互补误差函数构建建筑物目标的激光回波信号理论模型；
步骤2：将星载单光子激光雷达飞越城市或其他包含建筑物目标的区域时采集的坐标为沿轨位置和高程的光子点云数据，按照高程间隔、沿轨距离间隔划分为二维网格，统计落在各网格单元内的光子点云个数从而得到每个网格单元的点云密度，通过目标识别算法选取建筑物目标的光子点云沿轨、高程所在的选定区域，对选定区域内网格的每个单元的点云密度在高程方向上进行求和，计算激光雷达飞越建筑物时光子点云数据的实际测量值；
步骤3：在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云；
步骤1所述建立星载激光雷达接收的激光回波信号功率分布模型，具体定义如下：
激光接收信号函数表示为地面目标响应函数与发射激光脉冲函数的卷积；
发射激光脉冲在空间上服从二维高斯分布，时间分布近似为高斯函数；
所述星载激光雷达接收的激光回波信号功率分布模型，具体如下：
其中，Pr(x,y,t)表示构建的星载激光雷达平面坐标系下的地表光斑中心的平面坐标(x,y)处时间t的接收光学系统接收得到的激光功率，表示卷积运算，Pt(t)是发射激光脉冲的功率分布，表示为：
其中，t表示激光脉冲的飞行时间，E0表示发射激光脉冲的总能量，σt表示发射脉冲宽度；
将星载激光雷达沿卫星飞行方向定义为沿轨方向，将星载激光雷达垂直于卫星飞行方向定义为垂轨方向，通过沿轨方向、垂轨方向构建星载激光雷达平面坐标系下；
其中，(x,y)表示星载激光雷达平面坐标系下的地表光斑中心的平面坐标，x表示星载激光雷达平面坐标系下的地表光斑中心的沿轨方向的坐标，y表示星载激光雷达平面坐标系下的地表光斑中心的垂轨方向的坐标；
(u,v)表示星载激光雷达平面坐标系下的光斑内微面元的平面坐标，u表示星载激光雷达平面坐标系下的光斑内微面元的沿轨方向的坐标，v表示星载激光雷达平面坐标系下的光斑内微面元的垂轨方向的坐标；h表示垂直于星载激光雷达平面坐标系所在的平面并指向天顶的垂直高程，κ0为总衰减系数，表示地面的法向量，/>表示激光指向向量，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，ρ(u,v,h)表示(u,v,h)处的目标反射率，l表示卫星上的激光参考点与地表光斑内微面元的位置(u,v,h)的几何距离，如下式所示：
其中，c为大气中的光速，(xs,ys,hs)表示星载激光雷达平面坐标系下星载激光雷达的激光参考点坐标，xs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的沿轨方向的坐标，ys表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂轨方向的坐标，hs表示星载激光雷达平面坐标系下星载激光雷达的激光参考点的垂直于星载激光雷达平面坐标系的平面并指向天顶的垂直高程；
步骤1所述建立描述建筑物几何和反射特征的目标模型，具体定义如下：
假设建筑物目标表面具有一致的反射特征，对于近天底入射的激光脉冲，建筑物的反射表面为建筑物顶面，即假设接收光学系统接收的来自建筑物顶面上每一处微面的反射率相同；
所述建筑物顶面几何和反射特征的目标模型，具体定义如下：
γ(u,v,t)＝γ0δ(h-Ht(u,v)) (u,v)∈Duv
其中，δ(·)表示狄拉克函数，Duv表示建筑物在地面上的投影范围，γ(u,v,t)对应为S1.1的星载激光雷达接收的激光回波信号功率分布模型中的
对于被激光光斑照射的一个建筑物目标，激光指向矢量建筑物屋顶表面法向量/>几何距离l、总衰减系数κ0均为常数，使用衰减及反射常系数γ0综合表示建筑物顶面的反射特征；为了描述建筑物目标的空间分布情况，使用三维曲面方程h-Ht(u,v)＝0表示建筑物顶面的几何形状；
假设与激光轨迹相交的两条建筑物边界线即Duv的边界分别为规则的直线段l1和规则的直线段l2，(x1,y1)、(x2,y2)分别为前述设定的坐标系下l1、l2与激光轨迹的交点，x1和x2分别表示前述设定的平面坐标系下l1、l2与激光轨迹的交点的沿轨方向的坐标，y1和y2分别表示星载激光雷达平面坐标系下l1、l2与激光轨迹的交点的垂轨方向的坐标；
在星载激光雷达平面坐标系中y1＝y2＝0，星载激光雷达平面坐标系下l1、l2的直线方程分别表示为：
y＝k1(x-x1)
y＝k2(x-x2)
其中，k1和k2分别表示前述设定的平面坐标系下l1、l2的直线方程的斜率；
步骤1所述结合互补误差函数构建建筑物目标的激光回波信号理论模型，具体定义如下：
在时间轴t上积分得到建筑物目标反射到星载激光雷达接收光学系统的激光能量分布er(x,y)，表示为：
卫星沿固定轨道飞行，光斑中心坐标(x,y)的两个元素之间存在映射关系，er(x,y)简化表示为er(x)；
通过引入互补误差函数erfc(·)代替垂轨方向的高斯函数积分，进一步得到建筑物目标的激光回波信号理论模型，具体如下：
其中，E0表示发射激光脉冲的总能量，x表示星载激光雷达平面坐标系下地表光斑中心的沿轨方向的坐标，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，表示卷积运算，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标；
步骤2统计飞越城市或其他包含建筑物目标的区域时的星载单光子激光雷达数据点云密度，具体如下：
每个网格单元的点云密度计算公式如下：
其中，dmea(m,n)表示网格中第m行第n列的单元的点云密度，Lh表示网格划分的高程间隔，Lx表示网格划分的沿轨距离间隔，(m,n)表示网格中第m行第n列的单元，对应网格中心位置的高程和沿轨距离分别为mLh和nLx，Wh和Wx分别为统计点云密度时使用的窗口在高程和沿轨距离方向上的半宽，hpk和xpk分别表示整段数据中第k个点云所对应的高程和沿轨距离值，ε(·)表示单位阶跃函数，K表示整段数据的总光子点云数；
步骤2所述计算激光雷达飞越建筑物时光子点云数据的实际测量值，具体如下：
其中，Dtarg表示自动或人工选取的建筑物目标在沿轨、高程方向上所包含网格单元的集合；建筑物目标激光雷达沿轨实测的累计点云密度曲线smea(x)与建筑物目标激光雷达理论上沿轨接收的能量分布er(x)呈线性关系，即：
smea(x)＝C0·er(x)
其中，C0表示建筑物目标未知常系数，smea(x)表示计算得到的激光雷达飞越建筑物时光子点云数据的实际测量值；
步骤3所述在复频域构建激光雷达飞越建筑物时光子点云数据的实际测量值与回波信号理论模型之间的代价函数，具体如下：
将激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)与建筑物目标的激光回波信号理论模型er(x)在复频域中使用非线性最小二乘法进行拟合和参数估计；
激光回波理论模型er(x)经过离散傅里叶变换得到的频谱为Er(ω)，即其中/>表示离散傅里叶变换，具体公式如下：
其中，ω表示角频率，j表示复数虚部，E0表示发射激光脉冲的总能量，γ0表示衰减及反射常系数，σB表示激光束截面参数，等于激光半发散角和卫星飞行高度的乘积，erfc(·)表示互补误差函数，x1和x2分别表示星载激光雷达平面坐标系下对应建筑物两个边界的准确位置的沿轨方向的坐标，复频域中基于最小二乘法得到的代价函数f(x1,x2,k1,k2,C)如下式所示：
其中，星载激光雷达平面坐标系下建筑物两个边界在星载激光雷达沿轨方向的准确位置x1和x2、建筑物两个边界直线方程的斜率常数k1和k2以及建筑物目标常未知系数C为待估计的未知参数，发射激光脉冲的总能量E0、衰减及反射常系数γ0与建筑物目标常未知系数C0合并为未知常系数C来估计，即C＝C0E0γ0；
Smea(ω)是激光雷达飞越建筑物时光子点云数据的实际测量值smea(x)经过离散傅里叶变换得到的频谱函数，即ωi表示序列Smea(ω)中第i个频谱分量所对应的角频率，Smea(ωi)表示序列Smea(ω)中第i个频谱分量，I表示序列Smea(ω)的频谱分量总数的集合；
步骤3所述通过最小二乘迭代算法得到代价函数最小时的建筑物的边界位置提取对应的光子点云，具体如下：
使用Levenberg-Marquardt法，估计复频域中基于最小二乘法得到的代价函数的极小值点，准确得到建筑物边界在星载激光雷达沿轨方向上的位置参数x1*、x2*的最佳估计结果；光子点云数据包含每个光子点云的沿轨位置和高程，代入位置参数x1*、x2*到光子点云数据中，可以得到对应的建筑高度h1*和h2*。
2.根据权利要求1所述的激光雷达建筑物边界估测方法，其特征在于：
一种计算机可读介质，其特征在于，其存储电子设备执行的计算机程序，当所述计算机程序在电子设备上运行时，使得所述电子设备执行如权利要求1所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311146876.2A CN117347971B (zh) | 2023-09-05 | 一种激光雷达建筑物边界估测方法及计算机可读介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311146876.2A CN117347971B (zh) | 2023-09-05 | 一种激光雷达建筑物边界估测方法及计算机可读介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117347971A CN117347971A (zh) | 2024-01-05 |
| CN117347971B true CN117347971B (zh) | 2024-06-11 |
Family
ID=
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009143986A1 (en) * | 2008-05-27 | 2009-12-03 | The Provost, Fellows And Scholars Of The College Of The Holy And Undivided Trinity Of Queen Elizabeth Near Dublin | Automated building outline detection |
| CN106097311A (zh) * | 2016-05-31 | 2016-11-09 | 中国科学院遥感与数字地球研究所 | 机载激光雷达数据的建筑物三维重建方法 |
| CN111060899A (zh) * | 2019-12-31 | 2020-04-24 | 武汉大学 | 星地一体化激光雷达回波波形仿真方法及系统 |
| CN114325661A (zh) * | 2021-10-25 | 2022-04-12 | 武汉大学 | 基于激光角反射器点云信号的单光子激光雷达在轨检校方法 |
| WO2023015623A1 (zh) * | 2021-08-13 | 2023-02-16 | 复旦大学 | 一种多旋翼无人机载合成孔径雷达分段孔径成像及定位方法 |
| CN115932883A (zh) * | 2022-09-19 | 2023-04-07 | 国网河南省电力公司电力科学研究院 | 基于激光雷达的导线舞动边界识别方法 |
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009143986A1 (en) * | 2008-05-27 | 2009-12-03 | The Provost, Fellows And Scholars Of The College Of The Holy And Undivided Trinity Of Queen Elizabeth Near Dublin | Automated building outline detection |
| CN106097311A (zh) * | 2016-05-31 | 2016-11-09 | 中国科学院遥感与数字地球研究所 | 机载激光雷达数据的建筑物三维重建方法 |
| CN111060899A (zh) * | 2019-12-31 | 2020-04-24 | 武汉大学 | 星地一体化激光雷达回波波形仿真方法及系统 |
| WO2023015623A1 (zh) * | 2021-08-13 | 2023-02-16 | 复旦大学 | 一种多旋翼无人机载合成孔径雷达分段孔径成像及定位方法 |
| CN114325661A (zh) * | 2021-10-25 | 2022-04-12 | 武汉大学 | 基于激光角反射器点云信号的单光子激光雷达在轨检校方法 |
| CN115932883A (zh) * | 2022-09-19 | 2023-04-07 | 国网河南省电力公司电力科学研究院 | 基于激光雷达的导线舞动边界识别方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于自然地表的星载光子计数激光雷达在轨标定;赵朴凡 等;《红外与激光工程》;20201130;第49卷(第11期);全文 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110441758B (zh) | 一种星载线阵多波束测高激光雷达的在轨几何定标方法 | |
| Blacknell et al. | Geometric accuracy in airborne SAR images | |
| CN108627834A (zh) | 一种基于地基InSAR的地铁路基结构监测方法及装置 | |
| CN108414998B (zh) | 一种卫星激光测高仪回波波形模拟仿真方法及设备 | |
| Brunt et al. | Assessment of NASA airborne laser altimetry data using ground-based GPS data near Summit Station, Greenland | |
| Goulden et al. | The forward propagation of integrated system component errors within airborne lidar data | |
| US9097792B2 (en) | System and method for atmospheric correction of information | |
| CN102866393B (zh) | 一种基于pos与dem数据的sar多普勒参数估计方法 | |
| CN112098958B (zh) | 一种基于数字地图和气象水文信息的雷达杂波预测方法 | |
| CN113325422B (zh) | 天基测雨雷达目标定位及降雨信息三维处理方法和系统 | |
| CN102176003A (zh) | 一种机载激光雷达航测参数的优化设计方法 | |
| Wang et al. | UAV-Based P-Band SAR Tomography with Long Baseline: A Multi-Master Approach | |
| CN117347971B (zh) | 一种激光雷达建筑物边界估测方法及计算机可读介质 | |
| Spore et al. | Collection, processing, and accuracy of mobile terrestrial lidar survey data in the coastal environment | |
| Weiler | Bias correction using ground echoes for the airborne demonstrator of the wind lidar on the ADM-Aeolus mission | |
| CN116466361A (zh) | 一种基于无人机平台的机场净空测量装置及使用方法 | |
| CN117347971A (zh) | 一种激光雷达建筑物边界估测方法及计算机可读介质 | |
| Mitishita et al. | Study of stability analysis of the interior orientation parameters from the small-format digital camera using on-the-job calibration | |
| CN116500648A (zh) | 一种地基激光雷达目标区风廓线反演方法 | |
| de Haag et al. | Application of laser range scanner based terrain referenced navigation systems for aircraft guidance | |
| Li et al. | Within-footprint roughness measurements using ICESat/GLAS waveform and LVIS elevation | |
| Zhao et al. | Modeling and correcting building boundary in ICESat-2 spaceborne laser altimeter data considering the extended laser spot effect | |
| Salehi-Dorcheabedi et al. | Improving LiDAR height precision in urban environment: Low-cost GNSS ranging prototype for post-mission airborne laser scanning enhancement | |
| Zhou et al. | An improved method of AGM for high precision geolocation of SAR images | |
| Volkmer et al. | Consideration of the cloud motion for aircraft-based stereographically derived cloud geometry and cloud top heights |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |