PDERL：基于DEM的快速精确通视域分析算法介绍（一）

本文主要介绍我们团队研发的一种新的通视域算法——“临近度-方向-俯仰度”空间参考线算法（PDERL）。该算法已经在国际地理信息期刊《Earth Science Informatics》开放发表，欢迎各位对通视分析研究有兴趣的小伙伴全文围观：

出版社原文： PDERL: an accurate and fast algorithm with a novel perspective on solving the old viewshed analysis problem | SpringerLink

research gate : research gate 链接

Springer Nature： PDERL: an accurate and fast algorithm with a novel perspective on solving the old viewshed analysis problem（这个比较全，能直接下载到supplements）

github算法源码： https://github.com/blct-w/Pedrl-Algorithm

CSDN下载源码： 国内点击这个下载比较靠谱，原谅我没积分了，要点积分...

通视域算法大部分是基于规则格网DEM的，分为精确算法和近似算法，其中R3算法为经典的精确算法，该算法计算了视线经过的每条格网线对通视的影响，结果非常可靠，但由于其巨大的计算量难以运用于实际；而XDraw和参考面算法为经典的近似算法，其计算速度较快，自提出以来影响深远，常被用于实际分析运算，但由于计算时并未考虑所有网格线对通视的影响，这类算法及以它们为基础的快速算法不可避免的会出现误差，且我们经研究发现，它们的误差比率虽然不高，但误差点的分布却不均匀，极有可能聚集在一起，导致区域性通视判定错误，以致影响结果的应用，这个问题至今没有得到实质性改善。长久以来，传统DEM的采集精度与分辨率均较差，相较而言算法精度上的影响实际上并不大，故众多学者将研究重点放在了对算法速度的提升上，时至今日依然有学者在XDraw等算法基础上进行并行化研究，但随着现代传感器技术和存储能力的极大提升，DEM的采集精度与分辨率均已得到了很大的提升，粗略的近似计算已无法体现DEM数据质量提升所带来的优势，而精确算法的速度又无法达到实用的要求，通视域计算中速度与精度之间的矛盾，成为当前急需解决的关键问题。我们针对这个问题建立了“邻近度-方向-俯仰度”空间概念坐标系（proximity-direction-elevation 坐标系，简称PDE空间坐标系），在这个新的坐标系的基础上提出了PDE空间参考线算法(PDERL)，该算法采取分区计算的形式，在各区内分别将格网点转换到PDE空间，在构造和更新参考线的过程中判定每个格网点的通视情况，计算过程中每条DEM格网线对通视的影响均以函数的形式从中心向外围传递，无需内插，计算结果精确，运算速度快。经过大量测试，证实PDERL精度与R3完全相同，速度又能达到参考面算法的52%，XDraw算法的61%，是一种可以运用于大面积通视域分析的精确算法。

风险提示：我们团队研究的该算法属于职务成果，已被大学提交国家专利局进行专利申请。发布此工程的主要目的是为了方便国内做地理信息工程的程序员们读懂算法的核心思想，复现算法，应用于研究、学习等。但是计划将该算法用于商业目的的小伙伴们请注意：我们是算法发明人，不是专利申请人（申请人是大学组织），无法对您进行使用授权，如果您计划将其应用于商业项目，我们强烈建议您首先联系大学或发邮件给我代为咨询：blct_w@foxmail.com。

此外，如果您成功应用或改进了这个算法，我们会感到非常荣幸，如果方便请告知我们论文或系统的名称，非常感谢！

文章目录

前言

一、整体思路

二、PDE空间定义及性质

1.PDE空间定义

2.PDE空间的线性不变性

三、PDERL算法

1.计算分区

2.地形线在PDE空间的表达

3.俯仰度差的定义及通视判定函数

4.参考线的构造更新与通视判定法

总结

前言

该算法是目前你能找到的直接基于DEM规则格网通视分析中唯一精确的快速算法，对通视分析比较熟悉的小伙伴应该最熟悉的通视域算法是R2,R3,XDraw和参考面算法，以及基于这几个经典算法的大量并行算法。R3算法精度高，但计算量巨大，其它的几个属于近似算法，虽然速度比较快，但精度普遍不高，并且会出现错误点聚集的问题，而我们研制的PDERL算法精度和R3一样，是一种精确算法，但速度却与XDraw和参考面算法接近，非常实用。

写这篇文章是为了方便国内做地理信息工程的程序员们以中文的方式读懂算法的核心思想，便于你们在自己的程序中实现算法。本文只介绍算法核心思想和推导，学有余力的伙伴们欢迎直接阅读原文，了解一下通视域算法的前世今生。

一、整体思路

目标是否通视决定于目标与遮挡物之间的相对位置关系，在地形上如果遮挡物对目标形成遮挡必然同时满足如下三个遮挡条件：

(ⅰ)遮挡物在目标的前面；

(ⅱ)遮挡物和目标位于观察者的相同方位；

(ⅲ)遮挡物视线仰角大于目标视线仰角。

传统通视域算法均在地理空间中进行推导，地理空间坐标对目标的位置、大小等信息的表达比较简单，而对上述三个条件的表达较为复杂，不利于通视计算。本文提出PDE空间坐标系，直接对上述三个条件进行量化，该空间对目标的描述从传统的几何位置坐标转换为目标相对于观测者的邻近度（标记为 p）、方向（标记为 d）、俯仰度（标记为 e ）三个概念维度坐标。在PDE空间的同一方向上，e值越大，目标视线越高，p 值越大，目标距离观测者越近，如果两目标的 d 和 e 都相同，那么它们的视线将重合。

假设DEM上相邻两格网点连线能够表达这两点之间的地形，那么某行（列）DEM点连接而成的连续折线即可表达该行（列）对应地理空间的地形，如图2所示，本文称其为横向（纵向）地形线，本文以观测者的位置为界，将计算区域中的横向地形线和纵向地形线各分为两个区域。在同一区中，DEM上的通视域问题可以抽象为讨论距离观测者近的地形线（在PDE空间邻近度大）对距离观测者远的地形线（在PDE空间邻近度小）的遮挡情况，由于观测者的视线是以观测者为中心向四周发射的状态，在地理空间中求解这种遮挡关系极为复杂，但在PDE空间中较为简单。只需各区域分别转换入对应的PDE空间中计算，四个区域分别计算了行/列地形线对格网点的遮挡效果，某点仅有同时不被行/列地形线遮挡才能与观察者通视，如图2中点，仅在区域Ⅱ和区域Ⅲ中同时被判定为通视，才能判定为与观测者通视。根据以上原理完成四个区域计算后，再合并计算结果即可得到区域通视。

在计算邻近度较大的地形线对邻近度较小地形线的遮挡时，前者称为参考线，后者称为考察线。在PDE空间中，沿邻近度将参考线平移到考察线所在邻近度平面，两地形线交点所在的方向即是考察线被遮挡的状态发生改变的方向。以各交点方向值作为分界点将两地形线分为若干段，考察线上轴斜率较小的区段会被遮挡，位于这些区段上的格网点也会被遮挡，取两地形线中轴斜率较大的区段构成新的地形线，而这条新的地形线又将在计算邻近度更小的地形线时作为参考线。基于这个原理，在各区域中从距离观测者最近的两条地形线开始，按照上述步骤逐条对地形线进行计算，直至完成，各区域的计算思路如图1。任意格网点在各区域计算中都不被遮挡则与观测者通视。

图1 某区域的计算思路（图中OSP空间即PDE空间，懒得改了）

图2 计算分区

本算法建立在PDE空间坐标系的基础上，本文将在详细介绍PDE空间坐标系定义及性质的基础上推导PDE空间参考线通视域分析算法。

二、PDE空间定义及性质

1.PDE空间定义

为便于推导，定义Ⅰ号站心坐标系（Topocentric Coordinate System）：以观察者视点为原点，以DEM的纵向网格线偏北方向为Y轴（纵轴），以DEM的横向网格线偏东方向为X轴（横轴），以观察者所在的参考椭球面的法向量方向为Z轴（竖轴）建立右手空间直角坐标系。图3 站心坐标系中A点站心坐标为 (x,y,z)，与原点水平距离为 l ，视线方向角为 ，仰角为 。

本文将定义三个概念：邻近度、方向、俯仰度。

临近度（记为 p ）定义为：站心坐标系中考察点的横轴 x 值的倒数，定义式为公式(1)。邻近度的意义是：在方向相同时，邻近度大的对象（距离观测者近），有可能会对遮挡邻近度小的对象（距离观测者远），反之不然，邻近度是“遮挡条件(ⅰ)”的判定值。

 (公式1)

方向（记为 d ）定义为：考察点在站心坐标系中方向角 的正切值，定义式为公式(2)。如果两点 d  值不同，代表位于观测点的方向不同，方向不重叠的对象不构成遮挡关系。d 值可以作为“遮挡条件(ⅱ)”的判定值。

(公式2)

俯仰度（标记为 e ）定义为：考察点的相对高程 z  与邻近度 p  之积，定义式为公式(3)。根据该式，e  在地理空间中的意义为：(0,0,0)、 (x,y,z)  、 (x,0,z)  三点形成面的俯仰程度，故本文称 e  为“俯仰度”。斜率（标记为 s  ）是考察点在站心坐标系中仰角  的正切值，斜率越大意味着视线仰角越大，而俯仰度与斜率符合关系式(4)，则在方向 d  相同时，轴斜率越大，考察点的视线仰角就越大，轴斜率可作为“遮挡条件(ⅲ)”的判定值。

(公式3)

(公式4)

                                                                                          
PDE坐标系定义为：以方向(d)为横轴，以邻近度( p )为纵轴，以俯仰度( e)为竖轴，构成空间直角坐标系，位于地理空间的考察点转换到该坐标系形成“方向-轴斜率-邻近度”的空间关系。                                                                                                                   

2.PDE空间的线性不变性

对于任意站心坐标系中任意不在ZOY平面的直线，一定有合适的常数c,d,m,n,v,w 使该直线能够表达为如下方程：

(公式5)

将方程(5)各式左右同除以 x  得到方程(6)，而方程(6)正是直线转换到OSP空间之后的表达式。结合相关数学原理，已知直线的方向、轴斜率、邻近度之间互为线性关系，本文称此性质为OSP空间的线性不变性，该性质为本文的通视域算法奠定了理论基础。

(公式6)

三、PDERL算法

1.计算分区

为了便于计算，本文对计算区域进行分区。如图4所示，以观测点为界，把纵向地形线分为Ⅰ区、Ⅱ区，把横向地形线分为Ⅲ区、Ⅳ区。以Ⅰ号站心坐标系Z轴为旋转轴，按顺时针方向分别旋转180°、270°、90°所得三个坐标系分别定义为Ⅱ、Ⅲ、Ⅳ号站心坐标系。以观测点为中心，在Ⅰ、Ⅱ、Ⅲ、Ⅳ区分别建立Ⅰ、Ⅱ、Ⅲ、Ⅳ号站心坐标系。按此方法划分区域后，任意横/纵地形线从相应的站心坐标系转换到PDE空间后均位于邻近度为正常数的平面上。

2.地形线在PDE空间的表达

设 u 为DEM的纵向格网宽度的倒数，以Ⅰ区为例。地形线上任一小段均是由两临近网格点构造的直线段，根据PDE空间性质，一定有合适的常数 a,b 使得线段在PDE空间符合公式(7)，且地形线转换到PDE空间后位于某邻近度为常数的PDE平面上，也是一条连续折线。本文将转换到PDE空间后的线段记为  ，其中 分别为线段方向定义域的左右边界，a,b,称为 L 的参数。如果已知两相邻格网点的站心坐标分别为A 、B ，则常数 a,b  的值可以表达为(8)式。实际上b的值不必每次都求，只需记住第一个点的a,b值和每一段末尾的d值就可以完整表达一条PDE空间的地形线。

(公式7)

(公式8)

3.俯仰度差的定义及通视判定函数

如图5所示，设  参数已知，邻近度为  ，空间某点C的PDE坐标为 ，且有 ， 。记  与PDE平面 交于点D，交点坐标设为 。定义C点相对于 L 的轴斜率差为公式(9)，对应于图5中红色有向线段，根据轴斜率的定义易知， 则C点通视，反之不通视。如果用  标记C点的通视性，取值０代表不通视，取值1代表通视，则可定义C的通视性判别函数为公式(10)。

(公式9)

(公式10)

图5  俯仰度差的定义

4.参考线的构造更新与通视判定法

更新过程以Ⅰ区为例，如图6所示以 列格网点构造地形线，再转换到PDE空间形成初始参考线 ，以 列格网点构造地形线，再转换到PDE空间构成待考察折线 ，将 沿P轴平移到 所在的平面上，以 与 的交点为界将两地形线分为若干段，合并 与中轴斜率较高的区段构成新的参考线 ，而 又将与 列格网点构造的地形线通过上述步骤继续更新，直至区域计算结束。

参考线更新的关键是寻找两地形线之间交点所在的方向值。由于地形线是折线，交点应当分段求解。假设 为参考线， 为考察线，以 和 中所有折线拐点的方向值为分段点，将两地形线分割为若干段，设第 i  段末尾方向为  ，该段中 和 上的线段的斜率分别为  和  ，此段末尾俯仰度差记为  ，由于直线特点，有如下关系：

(公式11)

公式(11)称为轴斜率差的累积公式，通过该式只需求出 上各格网点的方向d  和参数 a  即可累积计算出 各格网点相对于 的轴斜率差，再通过(10)式判定格网点通视性，记录到一个与整体计算区域同维的通视矩阵中。如果 与 正负相反，则 和  在第 i 段存在交叉点，交叉点方向记为  ，由公式(12)求出即可。

(公式12)

某格网点必须同时在Ⅰ/Ⅱ区和Ⅲ/Ⅳ区中判定为通视，才能最终确定该点为通视。本文用与通视域DEM同维的0-1矩阵记录通视域分析结果，每个格网点的通视计算结果保存在列号相同的对应元素中，取值0代表不通视，1代表通视。Ⅰ、Ⅱ、Ⅲ、Ⅳ区的计算结果分别记录在与计算区域同维的通视矩阵 中，则最终的通视矩阵应为：

(公式13)

总结

PDERL算法是我们团队研发的新的通视域算法，该算法同三阶算法R3一样是一种精确算法，但是计算速度却和二阶算法XDraw相近，非常实用。以上主要介绍了PDERL的主要推导过程，算法实现代码可以在CSDN中下载得到。由于时间有限，本文仅介绍了该算法的推导和原理，后续如果有需要将讲解这个算法与XDraw和参考面算法的对比，特别是XDraw等传统基于LOS算法的错误点团聚问题，也可以讲解一下算法代码中细节的实现。