CN108319717B - 一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法 - Google Patents

Abstract

本发明涉及一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法。若设定一个栅格尺寸的值t，并获得一个矢量点所在地图的最大点位中误差，可根据本专利提出的公式计算该矢量点转为栅格点后，该矢量点落入该栅格点的概率。然后可计算该地图上m个矢量点转换为栅格点后，至少n(n≤m)个矢量点落入对应的栅格点的事件的概率；随后给t赋一系列不同值，重复上述步骤，得到t与概率的对应关系表；最后根据用户指定的概率，依据该对应关系表确定t的最佳范围，即最佳栅格尺寸范围。本发明能够科学地计算出最佳栅格尺寸范围。

Description

一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法

技术领域

本发明涉及地理信息系统领域，为一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法。

背景技术

目前在地理信息系统研究和应用领域，对于矢量点转栅格点的操作,还没有计算最佳栅格尺寸(栅格边长)范围的算法。而该操作作为本领域中绝大多数空间分析的第一步，可能会导致空间位置的不可靠性，其理由如下:

一般认为，一个矢量点A(x_A,y_A)的坐标的概率密度函数满足二元正态分布。那么，若令A(x_A,y_A)的均值为A₀(x₀,y₀)，根据二元正态分布的性质，A处于A₀周围或与A₀重合。

根据矢量转为栅格操作的算法，当A转化为栅格点G时，A₀应该处于G之中。但是A可能处于G之中(如图1所示)，也可能处于G之外(如图2所示，图1和图2中坐标系的原点均在G的中心，G的边长均为t，矢量点A(x_A,y_A)的2维坐标均值对应点为A₀(x₀,y₀)。)。图2对应的事件会导致矢量点转栅格点操作中的空间位置不可靠性。这种不可靠性会随着空间分析本身以及空间分析结果的研究和使用而逐级放大。因此，若盲目地选择栅格尺寸会使得所有后继众多空间分析(如查询分析,复合叠置分析,缓冲区分析等)的结果出现误差甚至错误,导致相应研究结果和空间决策的错误，造成人员，财力和物力的不可挽回的损失。

图1对应的事件的概率p可以由下式计算(那么图2对应的事件的概率为1-p)：

式中,ρ_A分别为点A(x_A,y_A)的横坐标x_A的中误差，纵坐标y_A中误差和横纵坐标间的相关系数。D₁为积分区域，对应于A和A₀均在G内的情况(如图1所示)；D₂为积分区域，对应于A在G外但A₀在G内的情况(如图2所示)；一般认为，x_A,y_A,x₀,y₀相互独立。然而，上式在实践中无法使用，因为ρ_A常常未知。

发明内容

有鉴于此，本发明的目的是提出一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，能够科学地确定最佳栅格尺寸范围。

本发明采用以下技术方案实现：一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，包括以下步骤：

步骤S1：设定一个栅格尺寸的值t和矢量点A(x_A,y_A)最大的点位中误差σ_{A_max}，可根据下式计算点A转为栅格点后，点A落入该栅格点的概率p：

上式中，(x₀,y₀)为(x_A,y_A)的均值。

步骤S2：计算m个矢量点转换为栅格点后，至少n个矢量点落入对应的栅格点的事件的概率p_m,n；

步骤S3：给t赋不同值，并重复步骤S1与步骤S2，计算当t为不同值时对应的p_m,n，进而得到t与p_m,n的对应关系表；

步骤S4：根据用户指定的概率p₁，在步骤S3得到的t与p_m,n的对应关系表中选出大于p₁的p_m,n所对应的t的范围，并将该范围作为t的最佳范围，即最佳栅格尺寸范围。

进一步地，步骤S1中公式的定积分采用蒙特卡洛法计算。

进一步地，步骤S1中所述σ_{A_max}根据《GB/T 14912-20051:5001:10001:2000外业数字测图技术规程》查得。

进一步地，所述步骤S2中，p_m,n的计算采用下式：

与现有技术相比，本发明有以下有益效果：

本发明能够通过上述方法科学地确定最佳栅格尺寸范围，使得栅格数据的研究者或使用者从空间分析的第一步起即可准确地掌握空间分析的可靠性(通过概率量化)，并可以定量地分析后继的空间分析(如缓冲区分析，复合叠置分析)的可靠性，并为决策者提供可靠的基础数据，从而避免社会和经济损失。

附图说明

图1为背景技术中A和A₀均在G内的坐标示意图。

图2为背景技术中A在G外而A₀在G内的坐标示意图。

具体实施方式

下面结合实例对本发明做进一步说明。

本实施例提供了一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，包括以下步骤：

步骤S1：设定一个栅格尺寸的值t和矢量点A(x_A,y_A)最大的点位中误差σ_{A_max}，可根据式(1)计算点A转为栅格点后，点A落入该栅格点的概率p：

式(1)中，(x₀,y₀)为(x_A,y_A)的均值。式(1)的推导过程如下：

假设从实际地图上得到点A(x_A,y_A)。根据测量平差理论有：

其中σ_A为点A(x_A,y_A)的点位中误差。当ρ_A未知时，我们可以假设及ρ_A＝0。将代入得到下式：

将ρ_A＝0及代入背景技术中概率p的计算公式得：

若将一副地图中σ_A的最大值记作σ_{A_max}，σ_{A_max}即为该地图所对应的最大点位中误差(可以根据《GB/T 14912-20051:5001:10001:2000外业数字测图技术规程》--以下简称国标《GB/T 14912-2005》，查得)。由于没有其它任何关于σ_A分布的信息,根据最大熵原理,我们只能假设地图上所有点的点位中误差σ_A(这里σ_A可以视为一个随机变量)均匀分布于区间[0,σ_{A_max}]上(以获得σ_A的最大分布熵)。根据均匀分布的性质可以计算出σ_A的数学期望值为E[σ_A]＝(σ_{A_max}+0)/2＝σ_{A_max}/2，将E[σ_A]＝σ_{A_max}/2代替式中的σ_A，即可得到：

推导完毕。

步骤S2：计算m个矢量点转换为栅格点后，至少n个矢量点落入对应的栅格点的事件的概率p_m,n；

步骤S3：给t赋不同值，并重复步骤S1与步骤S2，计算当t为不同值时对应的p_m,n，进而得到t与p_m,n的对应关系表；

步骤S4：根据用户指定的概率p₁，在步骤S3得到的t与p_m,n的对应关系表中选出大于p₁的p_m,n所对应的t的范围，并将该范围作为t的最佳范围，即最佳栅格尺寸范围。

在本实施例中，步骤S1中式(1)中的定积分部分采用蒙特卡洛法计算。具体计算方法如下：

设存在定积分其中Ω为该定积分的积分区域。那么Ω的大小可由下式计算:而计算定积分I时,可以采用蒙特卡洛法在Ω上均匀地随机采样,以获得N个的随机样本：根据概率论中的大数定理，有所以定积分I可以近似地使用下式来计算(式中N数值越大所得到的I数值越准确)。

在本实施例中，步骤S1中所述σ_{A_max}根据《GB/T 14912-20051:5001:10001:2000外业数字测图技术规程》查得。

在本实施例中，所述步骤S2中，p_m,n的计算采用下式：

以下结合具体实例来说明本发明的计算方法。

步骤S1(包括(1)与(2)两个小步骤)：

(1)假设从一幅1:2000的标准地理底图(用于城市规划)中选取了矢量点A。现需计算A转换为栅格点后，A落入后者这一事件的概率。从表1中查得该地图中允许的最大点位中误差为1m(实际距离)(见表中第三行),那么点A的最大点位中误差也应为σ_{A_max}＝1.0m(那么E[σ_A]＝(σ_{W_max}+0)/2＝(1.0+0)/2＝0.5m)。

表1地物点平面位置精度^*单位为m

表1摘自国标《GB/T 14912-2005》中3.7.1小节内“表3地物点平面位置精度”(部分)。根据国标《GB/T 14912-2005》，当测图用于城市规划或一般用途时可以选取表1中括号内的点位中误差数值。对于其它比例尺地图，也有相关最大点位中误差的国家标准，限于篇幅此处不一一赘述。

(2)假设式(1)中t＝0.5m，即栅格尺寸为0.5m。那么使用蒙特卡洛法计算式(1)中定积分的步骤如下：

步骤1：式(1)中定积分的积分区域的大小为:

步骤2：生成25500000个随机4维向量,每组向量表示为(x_A,x₀,y_A,y₀)。

其中，可以使用数学软件Matlab的随机数发生器命令random产生或excel软件的rand命令产生随机向量(x_A,x₀,y_A,y₀)。也可以使用其它软件产生随机向量(x_A,x₀,y_A,y₀)。随机向量(x_A,x₀,y_A,y₀)的组数25500000与所使用计算机的内存容量以及所使用的软件有关(本实施例中，使用Matlab软件产生随机向量,所使用的计算机内存容量为4GB；经过尝试发现25500000为该计算机产生随机向量组数的最大上限；计算机内存容量越大,该上限值越高，该定积分的计算结果越准确)。

步骤3：将步骤1计算出的积分区域大小V＝0.0625和步骤2中产生的25500000组随机向量代入式(1)中定积分部分得到：

步骤4：将I＝0.0340和t＝0.5m代入式(1)得到p＝0.0340/0.5²≈13.60％。

步骤S2：假设从步骤S1中使用的1:2000的标准地理底图(用于城市规划)中选取了100个矢量点(m＝100)。需要将它们转换为栅格点进行空间分析。计算至少50个点(n＝50)在转换为栅格点后,前者落入对应后者的事件的概率。据步骤S1的计算，p＝13.60％。将m＝100,n＝50和p＝13.60％代入式(2)有：

步骤S3：令t＝1.0，1.25，1.5，...,2.25m，重复步骤S1和步骤S2，得到相应的p_100,50(见表2)。

表2t与p₁₀₀₅₀对应关系表

步骤S4：当用户指定p₁＝95％(那么需要保证p_100,50≥p₁＝95％)时。根据步骤S3中实例部分的t与p_n对应关系表(见表2),建议选取的最佳栅格尺寸范围为大于等于1.75m(实际距离)(当栅格尺寸达到1.75m时，p_100,50＝97.96％＞95％＝p₁)。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (4)

1.一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，包括以下步骤：

步骤S1：设定一个栅格尺寸的值t和矢量点A(x_A,y_A)最大的点位中误差σ_{A_max}，根据下式计算矢量点A转为栅格点后，矢量点A落入该栅格点的概率p：

上式中，(x₀,y₀)为(x_A,y_A)的均值；

步骤S2：计算m个矢量点转换为栅格点后，至少n个矢量点落入对应的栅格点的事件的概率p_m,n；

步骤S3：给t赋不同值，并重复步骤S1与步骤S2，计算当t为不同值时对应的p_m,n，进而得到t与p_m,n的对应关系表；

步骤S4：根据用户指定的概率p₁，在步骤S3得到的t与p_m,n的对应关系表中选出大于p₁的p_m,n所对应的t的范围，并将该范围作为t的最佳范围，即最佳栅格尺寸范围。

2.根据权利要求1所述的一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，其特征在于：步骤S1中公式的定积分采用蒙特卡洛法计算。

3.根据权利要求1所述的一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，其特征在于：步骤S1中所述σ_{A_max}根据《GB/T 14912-2005 1:500 1:1000 1:2000外业数字测图技术规程》查得。

4.根据权利要求1所述的一种矢量点转换为栅格点时最佳栅格尺寸范围的计算方法，其特征在于：所述步骤S2中，p_m,n的计算采用下式：

其中，为组合数符号；