CN105550639A - 对地观测激光测高卫星高程控制点自动提取方法和数据处理方法 - Google Patents

Abstract

对地观测激光测高卫星高程控制点自动提取方法和数据处理方法。该激光高程控制点提取方法包括：采用对地观测激光测距值有效性评估方法，保留判定为无云的足印影像块，剔除掉判定为薄云或厚云的足印影像块的激光高程数据；以激光足印点的反射率ε<1作为筛选参数，筛选所述保留的足印影像块的激光足印点，其中，选择回波波形中仅有一个波峰、峰值大于阈值且波形拟合后的标准偏差σ≤3.2ns的激光足印点，其中，用于确定所述阈值的参数包括激光器的发射能量和接收口径。根据本发明的对地观测激光测高卫星高程控制点自动提取方法和数据处理方法，可以减少云层对激光测距的影响，保证激光测距的精度，有效提高激光高程参考数据的准确性。

Description

对地观测激光测高卫星高程控制点自动提取方法和数据处理方法

技术领域

本发明涉及卫星数据处理方法，更具体而言涉及一种激光足印影像的云雾量检测方法及基于回波波形处理和云雾量检测的高程控制点自动提取方法和对地观测激光测高卫星数据处理方法。

背景技术

随着经济的发展，三维地理信息已经在数字地球、城市规划、环境保护等各方面得到广泛应用。而遥感卫星技术的快速发展，使航天摄影测量已经成为继航空摄影测量后又一个快速获取三维地理信息的手段，特别是近年来随着三线阵立体相机等技术的发展，三维遥感信息的获取技术取得了长足的进步。

例如，资源三号是我国首颗高精度民用立体测图卫星，它的成功发射和有效应用打破了我国长期以来对国外高精度卫星影像的依赖，产生了巨大的社会和经济效益。但由于光学立体卫星立体测图方式以及卫星姿轨测量精度、相机畸变等方面的原因，导致其在无地面控制点的条件下，其高程测量精度还难以满足高精度测绘的需求。

地面控制点GCP(GroundControlPoints)是卫星遥感影像几何纠正和地理定位时重要的基准数据源。在遥感影像几何纠正处理过程中，为达到一定的纠正精度，一定数量的地面控制点是必不可少的，通过控制点的物方坐标和相应的像点坐标构建影像成像模型并求解模型参数或对已有成像模型进行优化补偿求解补偿参数，最终建立成像过程中物方与像方的正确转换关系。

传统作业过程中的控制点一般采用全野外测量的模式，需要经过“收集已有控制资料、测区踏勘、选点埋石、外业测量、内业整理”等一系列复杂的工序，虽然随着GPS-RTK等先进测量技术的发展，外业测量的工作量和复杂度大大减少，但必要的控制点野外测量工作仍不可避免。

此外，由于卫星遥感影像覆盖范围大(以资源三号卫星为例，覆盖范围为50公里×50公里)，要想获得均匀分布的地面控制点，通常需要在几百公里乃至上千平方公里的范围内施测，其外业工作量之大、耗费人力和物力之巨是不言而喻的。此外，在地震、洪水、泥石流等自然灾害多发区域或原始森林、沼泽、沙漠等人迹罕至地区，测量人员常常无法进入实地进行测量。

即使控制点是通过人工选取遥感影像和地形图同名点来获得，也存在效率低下、精度难以保障等问题。而对同一地区的不同时相或不同传感器影像进行纠正，又会出现重复选点的情况。

地面控制点一般分为：平高控制点、平面控制点和高程控制点。中国专利申请201310143369.3公开了一种多源异构遥感影像控制点自动采集方法，能够从多源异构影像自动提取遥感控制影像点(控制点影像片)，提高了控制数据获取的效率和精度，但其实质上仍为一种平面控制点的自动采集方法。

在地面控制点数据不足的情况下，卫星影像的区域网平差技术可以作为精确几何定位的一种重要手段。例如，中国发明专利申请201510191096.9公开了一种基于星载激光测高数据的卫星影像立体区域网平差方法。在该专利申请所公开的方案中，采用星载激光测高数据作为广义高程控制点数据库的控制数据。

开展全球测绘，是为保护我国领土主权、相关利益和地理信息安全的必要支撑。卫星遥感技术在获取境外地理信息方面具有独特的优势。目前，国家测绘地理信息局正围绕“一带一路”、“高铁走出去”等重大战略需求，逐步开展全球高精度立体测图。而进行全球高程控制点的测量与获取，则是开展全球测绘的重要技术保障。因此，为顺利开展全球测图相关的工作，必须优先开展全球高精度控制点的获取，建立相应的控制点数据库。利用激光测高卫星进行全球高程控制点获取，是在目前卫星遥感技术条件下，一种有效的技术手段。

星载激光测高是一种地面点高程测定技术，它以卫星为平台，搭载激光测高仪从空间、时间对地球进行观测，高精度、实时地测定卫星到被测物体之间的距离，并通过数据处理与分析，获得地球表面的地貌、植被覆盖情况、海面形态等信息。国外星载激光测高技术发展迅速，世界主要发达国家都在开展星载激光测高仪的研究。例如，美国于2003年发射了ICESat卫星，其上搭载的地球科学激光测高仪***(GeoscienceLaserAltimeterSystem，GLAS)是全球第一个用于连续观测地球的星载激光测距***，其主要科学目的是测量极地冰盖高程及其变化，云层和气溶胶的分布特征等。ICESat卫星在轨运行期间，获取了大量高精度的高程数据，其激光足印平面精度达10m量级，高程精度约15cm。

国内近年来十分重视星载激光测高研究。自主研发的探月工程系列卫星“嫦娥一号”和“嫦娥二号”卫星均搭载了高精度的激光测高仪，获取了月球表面三维立体影像。除了“嫦娥”系列卫星上搭载激光测高仪外，还没有用于对地观测的激光测高***，但已经列入后续对地观测激光测高卫星发展规划，如1:1万比例尺的高精度民用立体测图卫星、陆地生态***碳监测卫星等，均将搭载对地观测的激光测高仪。

对地观测激光测高卫星在激光束往返大气的过程中，会与大气分子及气溶胶产生散射、折射等现象，进而导致激光能量出现衰减、回波变形、测距精度下降等问题，其中云特别是厚云的影响尤为显著，为了保证激光测距的精度，需要尽量减少云层对激光测距的影响，保证获得的激光测距值有效。

发明内容

根据本发明的实施例的一方面，提供了一种激光足印影像的云雾量检测方法，其包括：利用大量的含云雾影像进行样本训练，获得云雾的优选灰度阈值及纹理特征值；计算整幅足印影像的灰度直方图，基于所述优选灰度阈值初步判断是否存在云雾，以及云雾的含量；计算足印影像的子块的灰度均值，当所述子块的灰度均值高于第一阈值时判定含云雾，当所述子块的灰度均值低于第二阈值时判定不含云雾；利用灰度共生矩阵计算灰度均值处于所述第一阈值和第二阈值之间的影像子块的纹理特征，与从所述样本训练获得的所述云雾纹理特征值进行对比；统计被判定为云雾的像素总数及其在所述整幅足印影像中的数量占比，确定云雾量。

根据本发明实施例的激光足印影像云雾量检测方法，可选地，在计算整幅足印影像的灰度直方图时，对影像进行直方图均衡化增强处理，均衡化的公式如下：

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\&Sigma;}}{P}_r\left(r_i\right)=\underset{i=0}{\overset{k}{\&Sigma;}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

其中，s_k为原图灰度值为k的像素经过变换后的新灰度值；P_r(r_i)是灰度值为i的像素频率；L是图像的灰度等级。

根据本发明实施例的激光足印影像云雾量检测方法，可选地，所述第一阈值时235，所述第二阈值是80。

根据本发明实施例的激光足印影像云雾量检测方法，可选地，利用灰度共生矩阵计算所述影像子块的纹理特征时所考虑的纹理特征包括：角二阶矩，同质性，对比度和相关性中的一个或多个。

根据本发明的实施例的另一方面，提供了一种对地观测激光测距值有效性评估方法，其包括：对于激光测高卫星下传的经解码后的原始数据进行***几何改正及初步的波形处理，根据几何对应关系，获得激光光斑对应的足印影像，基于硬件参数和地理坐标或者足印相机与激光测高仪的初步标定来实现激光和足印影像的基本配准；确定所获得的与激光光斑对应的足印影像的中心点位置，并确定足印影像中激光光斑对应的面积；将确定的影像进行裁剪或存入内存中形成足印影像块；采用前述的云雾量检测方法，对所述足印影像块进行云雾量检测，获得足印影像块的云雾量值；基于所述云雾量值判断所述足印影像块的激光测距值的有效性或可用性。

根据本发明的实施例的又一方面，提供了一种激光高程控制点提取方法，其包括：采用前述的对地观测激光测距值有效性评估方法，保留判定为无云的足印影像块，剔除掉判定为薄云或厚云的足印影像块的激光高程数据；以激光足印点的反射率ε<1作为筛选参数，筛选所述保留的足印影像块的激光足印点，其中，选择回波波形中仅有一个波峰、峰值大于阈值且波形拟合后的标准偏差σ≤3.2ns的激光足印点，其中，用于确定所述阈值的参数包括激光器的发射能量和接收口径。

根据本发明实施例的激光高程控制点提取方法，可选地，采用多个高斯函数叠加的方式对所述回波波形进行波形拟合，拟合公式如下式所示：

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}$

上式中t为时间，A_m,t_m,σ_m分别为第m个高斯函数的幅值、均值和标准差，ε为波形噪声值。

根据本发明的实施例的再一方面，提供了一种对地观测激光测高卫星数据处理方法，其包括：对激光的发射和接收波形进行平滑、去噪，提取波形特征参数，确定发射和接收波形重心对应的时刻，根据激光传输时间间隔，计算激光单向传输初始距离值；根据卫星姿轨参数、激光发射时间、初始距离值及激光几何定位模型，计算激光足印点的粗略三维坐标；根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，重采样生产出足印影像几何粗纠正产品；采用前述的云雾量检测方法，对足印影像进行云雾量检测，对检测到的云雾量值未超过预定标准的足印影像的激光足印点进行下一步处理。

根据本发明实施例的对地观测激光测高卫星数据处理方法，可选地，还包括：采用大气延迟改正模型，根据大气参数及所述激光足印点粗略三维坐标获得大气延迟改正值；根据几何检校获得的激光测距***性误差值、所述初始距离值、所述大气延迟改正值，计算精确的距离值；根据卫星姿轨参数、激光发射时间、所述精确的距离值及激光几何定位模型，计算激光足印点三维坐标；采用潮汐改正模型，计算潮汐改正值，对所述计算激光足印点三维坐标进行修正，获得激光足印点精确三维坐标。

根据本发明实施例的对地观测激光测高卫星数据处理方法，可选地，还包括：根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，采用地形数据生产足印影像正射纠正产品；判断云雾量值未超过预定标准的足印影像为无云，则对回波波形特征参数进行特征约束条件判断；如果满足所述特征约束条件，则用所述激光足印点精确三维坐标和所述足印影像正射纠正产品来组合形成激光高程控制点。

根据本发明实施例的对地观测激光足印影像的云雾量检测方法、基于云雾量检测和回波波形处理的对地观测激光高程控制点的提取方法和激光测高卫星数据处理方法，可以减少云层对激光测距的影响，保证激光测距的精度，有效提高激光高程参考数据的准确性。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍，显而易见地，下面描述中的附图仅仅涉及本发明的一些实施例，而非对本发明的限制。

图1是根据本发明实施例的激光足印影像的云雾量检测方法的流程图；

图2示意性地示出了根据本发明实施例的激光高程控制点精细提取流程；

图3A是激光测高卫星严密几何定位模型示意图，图3B是激光发射方向与本体坐标坐标系的夹角示意图；

图4A和图4B分别示意性地示出了激光测高仪的发射波形和接收波形；

图5示意性地示出了根据本发明一个实施例的星载激光测高数据的处理流程；

图6示意性地示出了激光测高卫星数据产品的分级体系；

图7示出了基本产品的示意性生产流程；

图8示出了标准产品的示意性生产流程；

图9示出了专题产品/高级产品的示意性生产流程。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例的附图，对本发明实施例的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于所描述的本发明的实施例，本领域普通技术人员在无需创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

除非另作定义，此处使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明专利申请说明书以及权利要求书中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。同样，“一个”或者“一”等类似词语也不表示数量限制，而是表示存在至少一个。

我国计划于2018年发射的1:1万高精度民用立体测图卫星，其将搭载一套激光测高***。该激光测高***作为辅助载荷，用于与CCD相机拍摄的光学影像联合平差，达到提升高分辨率卫星无控定位精度的目的。该套激光测高***还将配备足印相机，能够获取脉冲到达地面形成的足印影像，从而可以确定足印的真实位置，克服激光测高***平面定位误差大的问题。例如，激光测高***采用波长为1064nm的激光脉冲，该波长处于大气窗口波段，对大气有较好的穿透性，然而在面临云遮挡或大雾的情况时，激光无法穿透云雾到达地面，因而获得的足印位置可能是云雾表面的位置，这样的数据达不到观测要求，要予以剔除。

可选地，可以通过将获得的测高数据与已有的数字高程模型(DEM)对比，将差距超过一定阈值的测量数据判断认为测量的足印点位于云雾上，从而进行剔除。但是这种方法对于有一定局限性，一方面对DEM的精度有一定要求，另一方面在激光能部分穿透云雾到达地面，且能接收到地面回波信号时，该算法易失效，不能准确地反映激光测距值的可靠性。

根据本发明实施例的激光测距数据提取方法基于云雾量的检测结果对激光测距数据进行筛选。利用足印影像，采用遥感影像分类的思想进行云雾检测判别，还可以进一步结合激光回波波形特征，若足印路径上存在云及大雾时，且回波波形峰值较小、展宽较大的则认为测量数据存在较大误差或不可用。

以下具体说明利用足印影像进行云雾量检测的方法。遥感影像的云雾量检测本质上就是影像的分类，将云雾看做诸多地物的一种。一般而言，云雾量检测依赖的是云雾的两个特征：光谱特征和纹理特征。光谱特征反映云团在不同波段的辐射反射情况，包括云的辐射特性和云及地物的反射光谱特性。纹理特征可以被认为是，在局部范围内，图像灰度等级间的空间分布及其空间的相互作用。

由于激光测高***的足印影像不包含多光谱信息，根据本发明，考虑利用足印全色影像的云雾检测方法，可选地，采用基于灰度直方图的阈值法和基于灰度共生矩阵的检测法来检测云雾。进一步，为提升云雾检测的准确率，可以综合这两种方法进行云雾判别。

灰度直方图反映的是一副图像中各灰度级像素出现的频率，以灰度级为横坐标，以灰度级的频率为纵坐标，它是图像的一个重要特征，反映了图像灰度分布的情况。当整张足印影像的灰度直方图存在明显的二峰性时，则表明这个图像的较亮区域和较暗区域可以较好的分离，以这两个峰值作为参考设定阈值，可以得到较好的二值图像处理效果。

经过一定量的遥感影像样本训练获取云雾的灰度值特征的最佳阈值，在整幅遥感影像的灰度直方图中，阈值处于一定频率时即可判别有云。利用灰度直方图可以从总体上可判断含云量。将影像分割为若干子图块，对每个子图块统计其灰度均值，根据均值的大小，初步实现云雾判定。

由于云在图像上的灰度分布均匀、跳变程度小，纹理较粗且模糊，与部分地物的特征相似，比如沙漠与层云的纹理特征极为接近。为了增加云雾检测的准确度，可以对影像进行直方图均衡化增强处理，突出隐含的纹理细节。均衡化的公式如下：

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\&Sigma;}}{P}_r\left(r_i\right)=\underset{r=0}{\overset{k}{\&Sigma;}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

上式中，s_k为原图灰度值为k的像素经过变换后的新灰度值；P_r(r_i)是灰度值为i的像素频率；L即图像的灰度等级。

另一方面，纹理可以被认为是，在局部范围内，图像灰度等级间的空间分布及其空间的相互关系。灰度共生矩阵是对与图像保持某个距离的两个像素分别具有的特定灰度的统计而得到，它能反映出灰度的分布特性，同时也可以反映具有同样灰度或接近同样灰度的像素间的位置分布特性，是一种基于图像灰度联合概率矩阵的方法，即关于图像灰度变化的二阶统计特征。由于云雾检测是基于足印影像灰度值的检测，因此云检测方案可使用灰度共生矩阵作为其特征提取的依据。

灰度共生矩阵是从图像灰度值为i的像元(x,y)出发，统计与其距离为d，灰度值为j的像元(x+a,y+b)同时出现的频度P(i,j,d,θ)，数学表达式为：

P(i,j,d,θ)＝{[(x,y),(x+a,y+b)|f(x,y)＝i；f(x+a,y+b)＝j]}

其中，θ为灰度共生矩阵的生成方向。简单地说，d与θ决定了两个像素的位置和方向，在两者一定的条件下可得出灰度共生矩阵P：

$P=\left(\begin{array}{ccc}P(0,0)& ...& P(0,L-1)\\ ...& ...& ...\\ P(L-1,0)& ...& P(L-1,L-1)\end{array}\right)$

其中每一个元素P(i,j)代表灰度i与灰度j在给定位置和方向出现的次数，L表示灰度级数。只需要求出4个方向的值即可代表所有纹理的信息。

可使用以下4个特征量f₁，f₂，f₃，f₄作为纹理特征：

(1)角二阶矩f₁

$f_1=ASM=\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}{P}^2(i,j)$

反映了图像分布的均匀性和纹理的粗细程度：当大多数元素分布靠近主对角线时，说明图像灰度分布较均匀，此时纹理比较粗，所以f₁的值会较大；反之f₁较小时，说明纹理比较细。

(2)同质性f₂

$f_2=IDM=\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}\frac{P(i,j)}{1+{(i-j)}^2}$

逆差矩反映了图像纹理的同质性，如果图像局部纹理缺少变化非常均匀，则该值会比较大。

(3)对比度f₃

$f_3=\underset{n=1}{\overset{L}{\&Sigma;}}{n}^2\left(\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}P\right(i,j\left)\right)$

其中|i-j|＝n，图像对比度反映了纹理的清晰度。图像中，如果像素对的灰度差别越大，那么纹理的对比越强烈，视觉效果越明显。反之则表明纹理不明显。往往粗糙的纹理对比度比较小，所以对比度也可以反映纹理的粗细程度。

(4)相关性f₄

$f_4=COR=\frac{\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}(i-j)P(i,j)-{\mathrm{\&mu;}}_x{\mathrm{\&mu;}}_y}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

其中，

${\mathrm{\&mu;}}_x=\underset{i=1}{\overset{L}{\&Sigma;}}i\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j),{\mathrm{\&mu;}}_y=\underset{i=1}{\overset{L}{\&Sigma;}}j\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j)$

${\mathrm{\&sigma;}}_x=\underset{i=1}{\overset{L}{\&Sigma;}}{(i-{\mathrm{\&mu;}}_x)}^2\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j),{\mathrm{\&sigma;}}_y=\underset{i=1}{\overset{L}{\&Sigma;}}{(j-{\mathrm{\&mu;}}_y)}^2\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j)$

相关性反映纹理的方向，它衡量了灰度共生矩阵元素在行或列方向的相似度。若某种灰度在某一方向延伸长度越长，则相关性值在该方向就越大。当一幅图像中相似的纹理区域有某种方向性时，其相关性值较大。

利用大量的含云雾影像作为训练样本，获取云雾的纹理特征，作为云雾量判别的阈值，优选地采用角二阶矩作为纹理特征。

于是，将灰度直方图与灰度共生矩阵两种方法结合，进行足印影像的云雾检测，步骤如下：

利用大量的含云雾影像进行样本训练，获得云雾的优选灰度阈值及纹理特征值；

计算整幅足印影像的灰度直方图，初步判断是否存在云雾，以及云雾的含量，可选地，对足印影像进行灰度值均衡化增强；

计算足印影像的子块的灰度均值，根据该均值判定是否包含云雾，例如，灰度高于235的代表含云雾，低于80的不含云雾，对处于两阈值之间的影像子块进行进一步的检测；

利用灰度共生矩阵计算所述处于两阈值之间的影像子块的纹理特征，与样本训练获取的云雾纹理特征值进行对比；

统计被判定为云雾的像素总数及其在足印影像中的数量占比(百分比)，确定云雾量。

图1示意性地示出了上述激光足印影像云雾量检测方法的步骤。

基于上述的激光足印影像云雾量检测方法，根据本发明的实施例，提出了一种基于足印影像云检测的激光测距值有效性评估方法。

该激光测距值有效性评估方法的基本流程如下：

对于激光测高卫星下传的经解码后的原始数据进行***几何改正及初步的波形处理，根据几何对应关系，获得激光光斑对应的足印影像，基于硬件参数和地理坐标实现激光和足印影像的基本配准。如果已做了足印相机与激光测高仪的初步标定，也可将标定结果直接用于激光和足印影像的配准。

由于足印相机获得的影像远大于激光的实际地面足印大小，因此需要做一定裁剪。确定上述获得的与激光光斑对应的足印影像的中心点位置，并确定足印影像中激光光斑对应的面积S＝πR²，其中R为地面激光足印的半径，考虑到激光指向精度，可将R扩大一倍，以保证足印影像对激光足印点的全覆盖。将确定的影像进行裁剪或存入内存中形成足印影像块。

采用前述的云雾量检测算法，对足印影像块进行云雾量检测，获得足印影像块的云雾量值，并基于该云雾量值对云雾状况进行判断，例如，若云雾量<10％则判定为无云，云雾量>10％且<30％则判定为薄云，若云雾量>30％则判定为厚云。

之后可以根据云雾量值判定足印影像块的激光测距值的有效性或者可用性。例如，将判定为无云的足印影像块的激光测距值标记为优，将判定为薄云的足印影像块的激光测距值标记为良，将判定为厚云的足印影像块的激光测距值标记为差。进一步，标记为优的激光测距值完全可用，标记为良的激光测距值可视为待用，标记为差的完全不可用，不参与后续产品的生产。

激光测高卫星的一个重要应用就是能够获得全球范围的高精度的高程控制点数据，然而激光测高卫星的地面足印点的最终高程精度除受激光器本身测距精度影响外，还与卫星位置和姿态测量精度、大气影响、地面起伏大小等因素有关，因此必须对经过精确几何标定、大气改正、潮汐改正等处理后的激光高程参考点结合相关辅助参数进行筛选，保证最终保留可靠的激光高程参考数据作为专题产品。结合前述的基于足印影像云检测的激光测距值有效性评估方法，根据本发明的实施例，进一步提出了一种多参数约束的激光高程控制点提取方法，以实现激光高程控制点专题产品的有效生产。

该多参数约束的激光高程控制点提取方法的基本步骤如下：

保留判定为无云的足印影像块的激光测距值，剔除掉判定为薄云或厚云的足印影像块的激光高程数据；

由于接收脉冲能量易受大气后向散射、太阳杂光、曝光过度等影响，因此确定激光足印点的反射率ε<1(其中)作为筛选参数，筛选上述激光高程数据中的激光足印点；

并且选择回波波形中仅有一个波峰，即nPeaks＝1，峰值E大于一定阈值V(具体阈值需根据激光器的发射能量和接收口径等参数确定)，且波形拟合后的标准差σ≤3.2ns的激光足印点；

将筛选后得到的激光足印点作为激光高程控制点，联合其足印影像，构建激光高程控制点数据库。

图2示意性地示出了根据本发明实施例的激光高程控制点精细提取流程。在图2所示的流程中，进一步结合了地形参考数据和足印影像正射纠正产品来提高激光高程控制点提取的精确度。

如下说明对地观测激光测高卫星回波波形的界定和拟合方法。

为了更准确地反映地物的垂直分布，揭示地物的几何和物理属性，对地观测的激光测高卫星一般以很小的等采样间隔(如1GHz)，对发射激光脉冲与目标作用形成的后向散射能量按时间先后顺序进行采样记录，得到一个随时间变化的回波信号，这个回波信号被称为回波波形数据。回波波形包括振幅、脉宽、波峰数、后向散射截面等特征信息，基于回波波形数据进行波形特征分析，也可以构成激光测高数据处理中的一个的重要步骤。

激光测高卫星的基本原理是：通过卫星发射激光束经地面反射后由卫星接收，计算激光发射和接收的时间间隔t，光的传播速度为c，则激光单程传输距离p＝c*t/2，然后结合卫星上搭载的GPS和星敏感器获得的卫星位置和姿态信息，即可获得激光足印点的三维地面坐标。其严密几何模型如图2所示，其中P_laser为激光发射的参考点，P_GPS为GPS天线相位中心，O_Body为卫星质心或本体坐标系的原点，P_Ground为激光地面足印点。

卫星的原始轨姿参数测量均在国际天球参考***ICRS(InternationalCelestialReferenceSystem)下进行定义，而国际天球参考框架ICRF(InternationalCelestialReferenceFrame)是ICRS的具体实现，目前主要采用J2000坐标系。通常描述的地面三维坐标是在国际地球参考系ITRS(InternationalTerrestrialReferenceSystem)下定义的，如与ITRF2000符合度最好的WGS84坐标系，其中ITRF是ITRS的具体实现。国际天球参考框架ICRF到国际地球参考框架ITRF的坐标转换矩阵与地球自转、岁差、章动等有关，其中P(t)为岁差矩阵，N(t)为章动矩阵，R(t)为地球自转矩阵，W(t)为极移矩阵。

卫星本体坐标系定义为：卫星质心为原点，X轴指向卫星飞行方向，Z轴指向天顶方向，Y轴垂直于卫星轨道平面，与X、Z轴构成右手坐标系，如图3A所示。激光发射时指向方向与卫星本体坐标系存在一定的夹角，假设激光指向与本体坐标系Z轴的负向夹角为θ，在XOY平面上的投影与X轴正向夹角为α，如图3B所示。

若激光的测距值为ρ，激光参考点在本体坐标系中的坐标为：则在激光足印点在卫星本体坐标系下的坐标为：

$\stackrel{\&OverBar;}{\mathrm{\&rho;}}=\left(\begin{array}{c}Lx\\ Ly\\ Lz\end{array}\right)+{\left(\begin{array}{c}\mathrm{\&rho;}sin\mathrm{\&theta;}cos\mathrm{\&alpha;}\\ \mathrm{\&rho;}sin\mathrm{\&theta;}sin\mathrm{\&alpha;}\\ -\mathrm{\&rho;}cos\mathrm{\&theta;}\end{array}\right)}_{Body}$

GPS的相位中心一般与卫星质心不完全重合，两者之间存在一定的偏差，假定GPS的相位中心在卫星本体坐标系下的坐标为：卫星星敏感器本体坐标系与卫星本体坐标之间的旋转矩阵为星敏感器测得的是星敏本体系到J2000坐标系下的旋转矩阵

因此，激光测高卫星的严密几何定位公式为：

${\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{ITRF}=R_{J2000}^{ITRF}\&lsqb;{\left(\begin{array}{c}{X}_{GPS}\\ Y_{GPD}\\ Z_{GPS}\end{array}\right)}_{J2000}+R_{Star}^{J2000}{R}_{Body}^{Star}\&lsqb;\left(\begin{array}{c}Lx\\ Ly\\ Lz\end{array}\right)+\left(\begin{array}{c}(\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\sin \mathrm{\&theta;}\cos \mathrm{\&alpha;}\\ (\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\sin \mathrm{\&theta;}\sin \mathrm{\&alpha;}\\ -(\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\cos \mathrm{\&theta;}\end{array}\right)-\left(\begin{array}{c}Dx\\ Dy\\ Dz\end{array}\right)\&rsqb;\&rsqb;$

其中Δρ为因大气折射及硬件误差而引起的测距改正值。

地球参考框架坐标系ITRF下的空间直角坐标(X,Y,Z)与经纬度(B,L)及大地高H之间的对应关系如下式所示。

$\left\{\begin{array}{c}X=(N+H)\cos B\cos L\\ Y=(N+H)\cos B\sin L\\ Z=\&lsqb;N(1-e^2)+H\&rsqb;\sin B\end{array}\right.$

上式中，N为该点卯酉圈半径，e为参考椭球第一偏心率。N的计算公式为:

$N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}$

其中，a为参考椭球长半轴。

由于激光测高仪发射波形可近似为高斯脉冲，激光经地表面反射后回波波形可近似看做一次或多次高斯脉冲的叠加。图4A和图4B分别示意性地示出了激光测高仪的发射波形和接收波形。

对于多回波波形数据，采用多个高斯函数叠加的方式进行拟合，拟合公式如下式所示。

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}---\left(2.14\right)$

上式中t为时间，A_m,t_m,σ_m分别为第m个高斯函数的幅值、均值和标准差，ε为波形噪声值。

回波波形的高斯函数模型的3个参数可用向量形式表示为

c_m＝[A_m,t_m,σ_m](2.15)

则未知数可用向量表示为X,且X的维数为3p+1。

X＝[ε,c₁,c₂,…,c_p]^T

设定回波波形的采样数目为N，以GLAS为例，N＝544，则每个采样点的实际回波值表示为：R＝[r₁,r₂,...,r_N]^T(2.16)

采用式(2.14)计算出的每个采样点对应的波形模型的值为

W＝[w₁,w₂,...,w_N]^T(2.17)

定义采样值和波形模型的计算值误差l为

F(C)＝W-R(2.18)

为了确定模型的最优化时的未知参数向量C，即F(C)＝0；

对(2.14)的波形模型，求未知参数向量C偏导得：

$A=\frac{\&part;W}{\&part;C}=\left[\begin{array}{cccc}\frac{\&part;w_1}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_1}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_1}{\&part;c_p}\\ \frac{\&part;w_2}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_2}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_2}{\&part;c_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{\&part;w_N}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_N}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_N}{\&part;c_p}\end{array}\right]---\left(2.19\right)$

对于波形模型用一阶泰勒公式展开得：

$w_n^{j+1}=w_n^j+\frac{\&part;w_n^j}{\&part;\mathrm{\&epsiv;}}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}+\underset{m=1}{\overset{p}{\&Sigma;}}\frac{\&part;w_n^j}{\&part;c_m}{\mathrm{\&delta;}}_{c_m}---\left(2.20\right)$

其中：

$\frac{\&part;w_n}{\&part;\mathrm{\&epsiv;}}=1$

$\frac{\&part;w_n}{\&part;c_m}=\&lsqb;\frac{\&part;w_n}{\&part;L},\frac{\&part;w_n}{\&part;t_m},\frac{\&part;w_n}{\&part;{\mathrm{\&sigma;}}_m}\&rsqb;$

$\frac{\&part;w_n}{\&part;A_m}=\exp \&lsqb;-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;$

$\frac{\&part;w_n}{\&part;t_m}=\frac{A_m(t-t_m)}{{{\mathrm{\&sigma;}}_m}^2}\exp \&lsqb;-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;=\frac{A_m(t-t_m)}{{{\mathrm{\&sigma;}}_m}^2}\frac{\&part;w_n}{\&part;A_m}$

$\frac{\&part;w_n}{\&part;{\mathrm{\&sigma;}}_m}=\frac{A_i{(t-t_i)}^2}{{{\mathrm{\&sigma;}}_m}^3}\exp \&lsqb;-\frac{{(t-t_i)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;=\frac{(t-t_m)}{{\mathrm{\&sigma;}}_m}\frac{\&part;w_n}{\&part;t_m}---\left(2.21\right)$

对于式(2.18)，

有F(C)＝W(C)+Aδx-R,

其中，W(C)是根据待估参数C及模型(2.14)计算出的回波近似值。

约束条件为：F(C)＝0，则有

Aδx＝L

L＝R-W(C)(2.22)

则δx＝[A^TA]^-1A^TL，其中，考虑权重和先验值，公式可以写为：

δx＝[A^TPA+V₀]^-1A^TPL(2.23)

其中：P为权重矩阵，V₀为先验值矩阵,P_ij＝p_iδ_ij，[V₀]_jk＝p_ckδ_jk，p_i是第i次采样值的权重，p_ck为第k个向量c_k的先验值权重。

其中，δ_ij为克罗内克δ函数值：

${\mathrm{\&delta;}}_{ij}=\left\{\begin{array}{c}0,i\&NotEqual;j\\ 1,i=j\end{array}\right.$

计算出δx后，可以得到新的参数：

Cⁱ⁺¹＝Cⁱ+δx(2.24)

判断|δx|<d，其中d为设定的阈值。如果成立，则停止迭代，输出参数C；如果不成立，则采用新的参数C，回到(2.18)进行下一次迭代运算。

结合上述足印影像云雾检测和回波波形数据处理，对地观测激光测高卫星数据处理流程可描述为：

1)对激光的发射和接收波形进行平滑、去噪，提取波形特征参数，确定发射和接收波形重心对应的时刻，根据激光传输时间间隔，计算激光单向传输初始距离值；

2)根据卫星姿轨参数、激光发射时间、初始距离值及激光几何定位模型，计算激光足印点的粗略三维坐标；

3)根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，重采样生产出足印影像几何粗纠正产品；

其中，上述1～3步获得的波形特征参数、激光初始距离值激光足印点粗略三维坐标、足印影像几何粗纠正产品均为基本产品，需要处理所有的激光数据点，

4)对足印影像采用云雾检测模型进行处理，获得足印影像云雾检测产品；

5)如果足印影像云雾检测产品中云雾量超标(如比例超20％)，则回到第4步进行下一个激光点的处理；若云雾量不超标则进入第6步；

6)采用大气延迟改正模型，根据大气参数及第2步中的激光足印点粗略三维坐标获得大气延迟改正值；

7)根据几何检校获得的激光测距***性误差值、第1步中的初始距离值、第6步中的大气延迟改正值，计算精确的距离值；

8)根据卫星姿轨参数、激光发射时间、精确的距离值及激光几何定位模型，计算激光足印点三维坐标；

9)采用潮汐改正模型，计算潮汐改正值，对第8步中的结果进行修正，获得激光足印点精确三维坐标；

10)根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，采用地形数据生产足印影像正射纠正产品；

其中，上述4～10步获得的足印影像云雾检测产品、大气延迟改正值、精确距离值、潮汐改正值、激光足印点精确的三维坐标、足印影像正射纠正产品均为标准产品，

11)对第5步中云雾量不超标的足印影像云雾检测产品进行判断，如果无云，则进入第12步；否则，进入第17步；

12)对回波波形特征参数进行特征约束条件判断(例如采用前述的回波波形拟合)，如满足，则进入第13步；否则，进入第16步；

13)利用第9步和第10步的激光足印点精确三维坐标和足印影像正射纠正产品，组合形成激光高程控制点；

14)利用同区域多期测量获得的激光高程参考点进行高程变化分析，获得极地冰盖变化监测产品；

15)利用多轨测量获得的激光高程参考点构建数据库，形成激光高程控制点库产品；

16)对激光回波波形特征参数进行林业参数提取分析，获得林业专题产品；

17)对云量进行分析，结合第9、10步的结果，获得其他类专题产品。

其中，第14到17步产生的极地冰盖变化监测产品、激光高程控制点库产品、林业专题产品、其他类专题产品均属于通过进一步处理生成的专题类高级产品。

图5示意性地示出了根据本发明一个实施例的星载激光测高数据的处理流程。例如，根据卫星姿轨参数以及相机基本参数，建立经***几何校正后的足印影像，根据卫星姿轨参数、激光指向参数、激光初始测距值确定激光光斑在足印影像里的大致位置。

图6示意性地示出了激光测高卫星数据产品的分级体系。具体说明如下：

Level_0(0级)原始数据产品：指激光测高卫星下传的经解码后的原始数据，包括激光足印影像、卫星姿态测量数据、卫星轨道测量数据、激光发射波形数据、激光回波波形数据以及相关硬件参数等。原始数据产品一般不对外提供，仅供地面数据处理***使用。

Level_1(1级)基本产品：指针对原始数据产品经过***几何改正及初步的波形处理后的产品，包括波形的特征参数、激光粗略距离参数、足印影像几何粗纠正产品、激光足印粗略三维坐标等，其中将经事后精密处理后的姿轨参数归为基础产品。基本产品一般不对外提供，仅供地面数据处理***使用。

基本产品的示意性生产流程例如图7所示。如图7所示的生产流程，得到的产品为：波形特征参数、激光初始距离值、激光足印点初略三维坐标、足印影像几何粗纠正产品。在此生产阶段，按产品的生产要求，只要是格式完整、数据没有缺失，都应该生产，只是在后面加入质检，将数据质量较差的标记出来，不予发布。

Level_2(2级)标准产品：指经过精确几何标定、大气改正、潮汐改正等处理后的产品，包括大气改正参数、潮汐改正参数、足印影像云检测产品、足印影像正射纠正产品、激光足印精确三维坐标产品。标准产品可对用户发布。

标准产品的示意性生产流程例如图8所示。如图8所示的生产流程，得到的产品为：大气延迟改正参数、潮汐改正参数、精确距离值、激光足印三维坐标、足印影像云检测产品、足印影像正射纠正产品。

Level_3(3级)专题产品/高级产品：指在标准产品的基础上，针对用户特定需求开发的激光测高卫星高级产品，包括激光高程控制点库产品、林业专题产品、极地冰盖监测产品、其他类专题产品等。专题产品可对用户发布。

专题产品/高级产品的示意性生产流程例如图9所示。如图9所示的生产流程，得到的产品为：极地冰盖监测产品、激光高程控制点库产品、林业专题产品、其他类专题产品。

其中，前述的激光测距值有效性评估方法可以用于1级和2级数据(产品)的后续处理；前述的多参数约束的激光高程控制点提取方法可以用于2级数据(产品)的后续处理。

开展利用激光测高卫星进行全球高程控制点获取的研究，能在现阶段为我国提供一种获得境外陆地高程控制点的技术手段；为我国全球测绘提供高程控制点数据保障；另外还能用于大范围的获取和改善我国境内地面控制点的高程控制精度，为其它遥感设备提供检校数据；最后还能进一步提升国内遥感卫星数据的定位精度，拓宽国内卫星数据的应用领域。

根据本发明实施例的对地观测激光足印影像的云雾量检测方法、基于云雾量检测和回波波形处理的对地观测激光高程控制点提取方法和激光测高卫星数据处理方法，可以减少云层对激光测距的影响，保证激光测距的精度，有效提高激光高程参考数据的准确性。

以上所述仅是本发明的示范性实施方式，而非用于限制本发明的保护范围，本发明的保护范围由所附的权利要求确定。

Claims (10)

1.一种激光足印影像的云雾量检测方法，其特征在于，包括：

利用大量的含云雾影像进行样本训练，获得云雾的优选灰度阈值及纹理特征值；

计算整幅足印影像的灰度直方图，基于所述优选灰度阈值初步判断是否存在云雾，以及云雾的含量；

计算足印影像的子块的灰度均值，当所述子块的灰度均值高于第一阈值时判定含云雾，当所述子块的灰度均值低于第二阈值时判定不含云雾；

利用灰度共生矩阵计算灰度均值处于所述第一阈值和第二阈值之间的影像子块的纹理特征，与从所述样本训练获得的所述云雾纹理特征值进行对比；

统计被判定为云雾的像素总数及其在所述整幅足印影像中的数量占比，确定云雾量。

2.根据权利要求1所述的云雾量检测方法，其特征在于，在计算整幅足印影像的灰度直方图时，对影像进行直方图均衡化增强处理，均衡化的公式如下：

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\&Sigma;}}{P}_r\left(r_i\right)=\underset{i=0}{\overset{k}{\&Sigma;}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

其中，s_k为原图灰度值为k的像素经过变换后的新灰度值；P_r(r_i)是灰度值为i的像素频率；L是图像的灰度等级。

3.根据权利要求1所述的云雾量检测方法，其特征在于，所述第一阈值时235，所述第二阈值是80。

4.根据权利要求1所述的云雾量检测方法，其特征在于，利用灰度共生矩阵计算所述影像子块的纹理特征时所考虑的纹理特征包括：角二阶矩，同质性，对比度和相关性中的一个或多个。

5.一种对地观测激光测距值有效性评估方法，其特征在于，包括：

对于激光测高卫星下传的经解码后的原始数据进行***几何改正及初步的波形处理，根据几何对应关系，获得激光光斑对应的足印影像，基于硬件参数和地理坐标或者足印相机与激光测高仪的初步标定来实现激光和足印影像的基本配准；

确定所获得的与激光光斑对应的足印影像的中心点位置，并确定足印影像中激光光斑对应的面积；

将确定的影像进行裁剪或存入内存中形成足印影像块；

采用根据权利要求1-4中任一项所述的云雾量检测方法，对所述足印影像块进行云雾量检测，获得足印影像块的云雾量值；

基于所述云雾量值判断所述足印影像块的激光测距值的有效性或可用性。

6.一种激光高程控制点提取方法，其特征在于，包括：

采用根据权利要求5所述的对地观测激光测距值有效性评估方法，保留判定为无云的足印影像块，剔除掉判定为薄云或厚云的足印影像块的激光高程数据；

以激光足印点的反射率ε<1作为筛选参数，筛选所述保留的足印影像块的激光足印点，其中，

选择回波波形中仅有一个波峰、峰值大于阈值且波形拟合后的标准偏差σ≤3.2ns的激光足印点，其中，用于确定所述阈值的参数包括激光器的发射能量和接收口径。

7.根据权利要求6所述的激光高程控制点提取方法，其特征在于，

采用多个高斯函数叠加的方式对所述回波波形进行波形拟合，拟合公式如下式所示：

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}$

上式中t为时间，A_m,t_m,σ_m分别为第m个高斯函数的幅值、均值和标准差，ε为波形噪声值。

8.一种对地观测激光测高卫星数据处理方法，其特征在于，包括：

对激光的发射和接收波形进行平滑、去噪，提取波形特征参数，确定发射和接收波形重心对应的时刻，根据激光传输时间间隔，计算激光单向传输初始距离值；

根据卫星姿轨参数、激光发射时间、初始距离值及激光几何定位模型，计算激光足印点的粗略三维坐标；

根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，重采样生产出足印影像几何粗纠正产品；

采用根据权利要求1-4中任一项所述的云雾量检测方法，对足印影像进行云雾量检测，对检测到的云雾量值未超过预定标准的足印影像的激光足印点进行下一步处理。

9.根据权利要求8所述的对地观测激光测高卫星数据处理方法，其特征在于，还包括：

采用大气延迟改正模型，根据大气参数及所述激光足印点粗略三维坐标获得大气延迟改正值；

根据几何检校获得的激光测距***性误差值、所述初始距离值、所述大气延迟改正值，计算精确的距离值；

根据卫星姿轨参数、激光发射时间、所述精确的距离值及激光几何定位模型，计算激光足印点三维坐标；

采用潮汐改正模型，计算潮汐改正值，对所述计算激光足印点三维坐标进行修正，获得激光足印点精确三维坐标。

10.根据权利要求9所述的对地观测激光测高卫星数据处理方法，其特征在于，还包括：

根据卫星姿轨参数、足印影像、足印相机参数及足印影像几何定位模型，采用地形数据生产足印影像正射纠正产品；

如果判断所述云雾量值未超过预定标准的足印影像为无云，则对回波波形特征参数进行特征约束条件判断；

如果满足所述特征约束条件，则用所述激光足印点精确三维坐标和所述足印影像正射纠正产品来组合形成激光高程控制点。