基于阴影图的航空遥感地面覆盖快速分析方法

谢 鹏，汪超亮，李子扬

（1.中国科学院 定量遥感信息技术重点实验室，北京 100094；2.中国科学院 光电研究院，北京 100094；3.中国科学院大学，北京 100049）

基于阴影图的航空遥感地面覆盖快速分析方法

谢 鹏1,2,3，汪超亮1,2*，李子扬1,2

（1.中国科学院 定量遥感信息技术重点实验室，北京 100094；2.中国科学院 光电研究院，北京 100094；3.中国科学院大学，北京 100049）

航空遥感系统在地形复杂地区侧视成像作业时，有可能会因地形遮挡而出现成像盲区，在任务规划阶段进行载荷成像盲区检测是解决该问题的有效方法。鉴于现有基于地形可视域分析的盲区检测方法存在计算效率低、通用性差等缺点，提出了一种基于阴影图算法的航空遥感地面覆盖快速分析方法。实验结果表明，该方法不仅具有效率高、通用性强的特点，而且可定量分析载荷成像地面覆盖率，为航空遥感作业任务规划的盲区检测与覆盖评估提供了很好的支持。

阴影图算法；载荷成像地面覆盖分析；航迹评估；三维可视化仿真

航空遥感作为一种重要的高分辨率遥感数据获取手段，在资源探查、灾情评估、军事侦察等方面得到广泛的应用。但航空遥感系统在地形复杂地区侧视成像作业时，可能会因地形遮挡而出现成像盲区，因此在航空遥感作业任务规划阶段，需要先进行载荷成像的盲区检测，再采取针对性的措施保证数据获取的完整性，以提高航空遥感作业任务的成功率。

地形视域分析方法是目前航空遥感任务规划阶段载荷成像盲区检测的重要手段，其为复杂地形环境下航空遥感任务地面覆盖分析提供了很好的技术支持。为提高地形视域分析的效率，国内外学者开展了大量的研究工作。现有的视域分析方法大多以LOS（line of sight）算法[1]为基础，该算法利用点与点的几何关系判断地形可视性，通过比较视线上相应位置的高度和对应地形位置的高度判断是否可视。高程内插方法是影响LOS效率的关键所在，常见内插算法有R3算法[2-4]、R2算法[4]、Xdraw[4]算法等。DYNTACS算法[5]采用非均等划分视线段的双线性插值算法进行视域分析，王智杰等对DYNTACS算法进行了改进[6]，根据视线的倾斜方向，只计算与横向格网或与纵向格网的交点，提高了算法效率。应申等基于DEM数据的特征，提出了一种双增量地形可视域[7]分析方法。此外，还有一些研究采用并行计算方法进行可视计算，如基于参考面的可视域算法[8]，通过点与点形成的空间关系判断点与所有目标点是否可视，很大程度上提高了视域的计算效率。上述可视分析算法均需进行地形插值计算或格网的相交计算，计算效率不高，当地形范围和分辨率提高时，会产生过量计算的问题，难以实现地形可视域的快速分析。

现有的地形视域分析方法大多基于CPU实现，在分析地形复杂、覆盖范围较大的区域或高分辨率地形数据时，计算冗余、效率低下，无法直接应用到航空遥感成像地面覆盖快速分析。本文为满足航空遥感成像盲区快速检测的迫切需求，通过研究三维渲染管线和阴影图技术[9]，提出一种基于阴影图的地面覆盖快速分析方法，将大量图形计算任务通过GPU分析执行，充分发挥GPU三维渲染能力，实现载荷在复杂三维场景的覆盖快速分析，大大提高了视域分析的效率和三维可视化的显示效果。

1 航空遥感地面覆盖快速分析

1.1 载荷仿真成像深度图构建

航空遥感载荷可分为被动型和主动型。被动型载荷接收成像视场内的地物反射或发射的电磁波信号；主动型载荷自身发射电磁波，然后接收由地物反射的电磁波信号。不论是哪种类型的遥感载荷，在侧视成像时都有可能受地形遮挡的影响，导致电磁信号被遮挡而无法成像，形成载荷的成像盲区。在计算机图形学中，通过光源模型进行场景阴影渲染时也有类似的过程，光线由光源主动发射，通过与场景地形计算获取到阴影区，最终在屏幕上进行显示。两种类型载荷成像时的地面覆盖计算均可通过光源模型的地形阴影渲染过程进行仿真计算，这里采用的光源模型为聚光灯模型。聚光灯模型作为一种位置性光源，按照一定的位置和方向进行场景的光照阴影处理，在进行三维地形场景渲染过程中，通过载荷成像时的位置和姿态计算得到聚光灯模型的位置和方向参数。三维渲染管线通过聚光灯模型按照摄像机成像模型[9]仿真载荷成像地面覆盖情况，并将成像视场内的三维场景信息投射到二维的成像平面，其对应的瞬时视场内地形覆盖仿真计算借鉴阴影图[10]算法。

利用三维渲染的方法进行载荷的盲区检测，首先需要使用载荷的成像参数构建三维渲染场景的成像深度图，深度图构建过程涉及不同的坐标系变换，其对应的变换公式为：

式中，Pv'为变换后的图像纹理坐标；Mr为归一化矩阵；Mproj为投影矩阵；Mper为视图矩阵；Pv为变换前的世界坐标。Mr矩阵将复合变换结果中所有的分量取值范围限定为[0,1]，使之适合纹理坐标索引的合理范围。文中以面阵载荷为例，通过分析图1所示的载荷成像模型，构建出其成像所涉及的变换矩阵。

假定面阵光学载荷系统的焦距为f，底片的高度为H，可根据载荷成像的几何模型（图1）计算瞬时视场角，其计算方法为：

面阵载荷成像视锥体投影矩阵Mproj描述了载荷成像视锥的大小，可用视场角fovy和成像视锥体宽高比aspect表示， 即

载荷的成像位置和姿态对载荷成像地面覆盖的影响可通过视图矩阵近似表达。假定载荷成像位置为eye，载荷成像中心位置为center，以及垂直于成像相机向上的向量up，3个参数唯一确定了成像相机的位置和姿态。在构建视图矩阵之前通过这3个矢量构建相机的成像坐标系，即

然后在相机成像坐标系的基础上构建出载荷成像的视图矩阵，即

通过上述投影变换矩阵和视图变换矩阵，可近似模拟面阵光学载荷的成像模型，同时结合式（1）、（2），渲染得到载荷在当前成像状态下的深度图（图2）。

图1 载荷成像三维几何模型

图2 模拟成像深度图

1.2 载荷成像可视区域三维渲染

构建出载荷成像瞬时视场内三维场景的深度图后，再根据深度信息判断场景区域是否可见。具体方法为：首先通过将当前视点下的像素转换到光源坐标系下，计算每个像素与光源的距离；然后将距离值和阴影图中对应的深度值进行比较，确定当前像素是否处于阴影中；最后根据比较的结果对载荷成像盲区和载荷成像区域分别进行着色处理，从而实现对载荷成像盲区的检测。

三维渲染管线对载荷成像区域进行三维场景渲染着色时，需区分成像可视区域、成像盲区和成像背景区域，文中采用3种不同的颜色进行着色处理，如绿色表示成像可视区域、红色表示成像盲区，蓝色表示成像背景区域。通过着色渲染后，不但能区分地形区域可见与否，而且能区分载荷的成像范围和背景范围，便于后续的定量分析。载荷成像盲区检测渲染所对应的伪代码如下：

算法名称：载荷成像盲区检测渲染算法

输入：三维场景的顶点，深度纹理，纹理坐标gl_ TexCoord；

Fragment Shader:

{

Vec4 shadowcolor，visiblecolor，outscopecolor； //定义3种不同着色颜色；

//深度比较，当前片元是否可见，结果保存在compareResult中；

float compareResult = shadow2DProj( osgShadow_ shadowTexture, gl_TexCoord[1] )；

//归一化坐标；

Vec3 outXYZ=vec3(gl_TexCoord[1].x/gl_TexCoord[1]. wgl_TexCoord[1].y/gl_TexCoord[1].w，gl_TexCoord[1.z]/ gl_TexCoord[1].w)；

//判断位置进行着色；

If((outXYZ.x＞1)||(outXYZ.x＜0)||(outXYZ. y＞1)||(outXYZ.y＜0) ) //表明不在成像范围；

{

gl_FragColor=outscopecolor；

}

Else if((outXYZ.x＞0)&&(outXYZ.x＜1)&(outXYZ. y＞0)&(outXYZ.y＜1)) //表明在成像范围；

{

//判断是否可视，并着色；

gl_FragColor=mix(shadowColor,visibleColor,compar eResult)；

}

}

1.3 载荷成像覆盖参考面构建

虽然通过三维渲染管线的渲染区分了载荷的地面覆盖区域内的成像范围和成像盲区，但由于渲染管线作为一个状态机，保持了高度封装的特性，未提供访问每一块渲染区域的接口，无法直接获得渲染区域的坐标信息。因此本文在场景中增加了一个额外的“虚拟”相机，用于记录渲染管线的渲染结果，为下一步的定量分析提供二维栅格数据。此相机采用平行投影的投影方式，垂直于地面向下投影，如图3所示。采用平行投影的原因有两方面：在平行投影方式下，成像视锥是一个平行的长方体，可最大范围覆盖成像区域，不会产生成像盲区；在平行投影方式下，视口图像不会发生扭曲变形。

图3 “虚拟相机”的平行投影方式

1.4 载荷成像覆盖率计算

覆盖率是指载荷在瞬时视场内可观测到的地面面积占载荷成像投影于参考面所涵盖范围面积的百分比。覆盖率计算的可以有两种定义方式：①将载荷成像区域的地形作为参考面，载荷当前位置下投影于地形参考面的成像外边界范围作为覆盖率计算的总面积，如图4a所示黄色折线围起来的绿色区域和红色区域；②采用成像区域平均海拔高度平面作为参考面，将载荷当前位置下投影于海拔平均高度参考面上的成像范围作为覆盖率计算的总面积，如图4b所示黄色折线圈起来的绿色区域。在实际应用中，根据需求选择不同的覆盖率计算方法。

载荷成像覆盖率计算是通过分析§1．2三维渲染后的二维栅格图像获取的。假定栅格图像的宽度与高度分别为width和height，共有width×height个像元。设二维栅格图像中红色像素数目为a（表示载荷成像盲区的不可见像素数目），绿色像素数目为b（表示载荷成像可视区域的可见像素数目），蓝色的像素数目为c（表示不在成像范围的像素数目），a+b+c=width×height；将载荷在成像区域平均海拔高度平面上成像的可见像素数目记为d。根据覆盖率的两种定义，其对应的计算公式分别为：

式中，Φa表示第一种参考面定义下的载荷成像地面覆盖率；Φb表示第二种参考面定义下的载荷成像地面覆盖率。

图4 参考面的三维渲染结果

2 实验与结果分析

2.1 实验环境

实验计算机使用Windows 7操作系统，硬件配置为CPU i3 3．4GHz、显卡NVIDIA GeForce GT620 M、内存4 GB，开发环境OpenSceneGraph 3．0．1，开发工具Visual Studio 2010。实验采用内蒙古包头地区的DEM数据：经度范围108．775 6～109．203 8，纬度范围40．618 4～40．766 0。

2.2 实验结果

实验采用面阵CCD载荷参数：焦距34．51 mm、底片的元件长7．1 um。成像位置为（42．708 2N，109．100 7E，3 000 m），其视场角计算公式为：

面阵CCD载荷45°侧视角成像的效果如图5所示。

图5 侧视成像效果图

通过阴影图算法的渲染分析，得到载荷在当前成像位置下的仿真计算结果，如图6所示。其中，a图像为视域分析的结果，红色区域为载荷成像盲区，绿色区域为载荷成像可视区域，b图像为左侧对应方框的局部放大图。

图6 视域分析结果

通过本文方法对载荷成像地形区域的三维渲染，得到载荷地面覆盖二维栅格图像，如图7所示。其中，a图为实际地形参考面的成像渲染结果，b图为采用平均海拔高度为1 270 m的参考面的渲染结果。二维栅格图像相对应的像素统计信息如表1所示。

图7 载荷成像渲染结果

表1 像素数目统计表

按照第一种参考面的定义，计算得到当前载荷瞬时成像覆盖率为：

按照第二种参考面的定义，计算得到当前载荷瞬时成像覆盖率为：

2.3 算法效率对比实验

为了验证本文方法的计算效率，将本文提出的方法分别与传统的R3视域分析算法、商业软件ArcGIS10．0的可视性分析模块的可视域分析效率进行了比较。实验随机抽取不同视点，采用大小为325×224、650×507、1 446×520的高程数据进行了比对，3种算法的处理时间如表2所示。

表 2 不同分辨率的DEM数据下算法的处理时间/s

由表2可知，随着地形数据精细程度的增加，R3算法和ArcGIS分析算法的计算时间会越来越长，本文算法在计算效率上有明显优势，且对DEM分辨率不敏感，处理时间比较稳定。

由于在实际应用中需考虑盲区分析的实时性，因此本次实验同时也分析了深度图分辨率对算法效率的影响，分别分析了512×512、1 024×1 024、2 048×2 048、4 096×4 096四种分辨率下的深度图算法效率。

从图8可知，随着深度图分辨率的增加，算法的效率呈明显降低的趋势，当深度图分辨率达到4 096×4 096或更高时，很难进行实时的覆盖分析。所以在进行载荷盲区实时检测时，需根据需求和计算机配置条件选择合适分辨率的深度图，在效率与精度方面进行折衷考虑。由于实时进行场景三维渲染分析的帧率至少为24 Hz，文中实验计算机配置条件下选择1 024×1 024的深度图分辨率较为合适。

图8 深度图分辨率对计算效率影响

3 结 语

本文借鉴了阴影图算法，将其成功应用于航空遥感成像覆盖快速分析，实现了成像过程的盲区分析检测，在保持一定精度的基础上，实时显示当前载荷位置下的成像盲区以及快速计算成像覆盖率，提高了航空载荷盲区检测分析的效率。在未来的工作中，可从以下3个方面进行深入研究：①在深度图分辨率较低时，渲染结果会出现锯齿形边缘现象，可进行渲染显示的抗锯齿优化；②结合载荷的飞行航迹，估算载荷在整条航迹上的成像覆盖情况；③该方法可扩展应用到建筑规划、景区观测点选址、手机基站布设等领域，可根据具体应用领域的特点进行实时分析算法优化。

[1] Chen C, Zhang L, Ma J, et al． Adaptive Multi-resolution Labeling in Virtual Landscapes [J]． International Journal of Geographical Information Science, 2010, 24(6)∶ 949-964

[2] Shapira A． Visibility and Terrain Labeling [J]． Master's Thesis, Rensselaer Polytechnic Institute, 1994, 8(1)∶ 15-17

[3] Franklin W R, Ray C K, Mehta S． Geometric Algorithms for Sitting of Air Defense Missile Batteries Research Project for Battelle, Columbus Division[J]． Delivery Order, 1994, 27(5)∶41-44

[4] 缪贲术,屈军亮． 基于Xdraw 算法的SAR 航线规划[J]． 舰船电子工程,2011,31(10)∶ 52-55

[5] De Floriani L, Marzano P, Puppo E． Line-of-sight Communication on Terrain Models [J]． International Journal of Geographical Information Systems, 1994, 8(4)∶ 329-342

[6] 王智杰,印欣,朱良家． 地形通视性算法改进[J]． 情报指挥控制系统与仿真技术, 2004 (2)∶ 54-57

[7] 应申,李霖,梅洋,等．增量法地形可视计算与分析[J]．测绘学报, 2007,36(2)∶192-197

[8] Zhi Y, Wu L, Sui Z, et al． An Improved Algorithm for Computing Viewshed Based on Reference Planes [C]．Geoinformatics, 2011 19th International Conference on IEEE, 2011

[9] 苏国中,郑顺义,张剑清,等． OpenGL 模拟摄影测量方法研究[J]．中国图像图形学报,2006,11(4)∶ 540-544

[10] Williams L． Casting Curved Shadows on Curved Surfaces[C]．ACM Siggraph Computer Graphics, ACM, 1978

6 结 语

在本检定实例中，电子罗盘显示的定向角度值为109°30'16"，示值误差小于最小允许误差的15'，远远大于本校准方法的扩展不确定度，对于GPS电子罗盘，是一种可靠的校准方法。

参 考 文 献

[1] 田湘,范胜林,刘建业．基于GPS的短基线姿态系统研究[J]．传感器与微系统,2007,26(10)∶29-31

[2] 王建斌．GPS方位角及其在导线测量中的应用研究[J]．北京测绘,2005(3)∶39-42

[3] 熊建明,王宇飞．GPS短边方位测量若干问题的探讨[J]．四川测绘,2002,25(3)∶117-118

[4] 熊建明．GPS短边方位测量的精度分析[J]．测绘通报,2000(11)∶18-20

[5] 王解先,刘慧芹,唐立军．不同站心地平坐标系下的坐标归算[J]．工程勘察,2005(5)∶58-59

[6] 孔祥元,郭际明,刘宗泉．大地测量学基础[M]．武汉∶武汉大学出版社,2001

[7] 曹金国,廖望东,李建伟,等．GPS短边方位角测量精度定量分析研究[J]．全球定位系统,2013,38(1)∶18-19

第一作者简介：彭友志，工程师，主要从事GPS和全站仪等测绘仪器的计量和相关科研工作。

P237

B

1672-4623（2016）01-0071-05

10.3969/j.issn.1672-4623.2016.01.021

谢鹏，硕士，主要从事遥感信息处理的研究。

2014-10-21。

项目来源：国家高技术研究发展计划资助项目（2013AA122105）。（*为通讯作者）