矿业权管理信息系统中勘查区块数计算的实现算法

1、矿业权管理信息系统中勘查区块数计算的实现算法一、引 言 根据国务院制定的矿产资源区块登记管理办法，我国对矿产资源勘查区域采用区块方式进行管理。“矿产资源勘查区块的划分均以 1 : 5 万图幅为基础，按经差 1、纬差 1划分成基本单位区块；以基本单位区块为基础，按经差 30 、纬差 30划分成 A、B、C、D四个四分之一区块；以四分之一区块为基础，按经差 15 、纬差 15 划分成 1、2、3、4四个小区块”（摘自矿产资源勘查区块划分及编号办法）。 与之相应，矿业权管理信息系统的探矿权子系统和油气勘查开采子系统在数据库结构中，采用经纬度坐标保存地质调查区域、勘查区域和油气采矿区域等的空间位置，并

2、设“基本区块数、四分之区块数、小区块数和折合基本区块”等字段记录其区块信息，以上区块数的计算由系统根据区域坐标自动完成。 区块数的计算属于特定条件下的特殊应用，现在比较流行的 GIS（地理信息系统）软件，没有提供能直接获得结果的现成函数。矿业权管理信息系统在基于 GIS 平台计算区块数没有达到预期效果的情况下，从速度、效率、成本等方面考虑，决定采用自行研制的算法解决区块数计算问题。 二、算法思想 本算法从分析勘查区块划分的特点入手，通过建立区域坐标矩阵，将矢量坐标区域转化为以 15 为基本单元的栅格矩阵，并采用自下而上扫描的方式，计算出勘查区域内的栅格个数，从而得到坐标区域所覆盖的区块数。 1

3、勘查区块管理的特点 依据矿产资源勘查区块登记管理办法通过对勘查区域坐标的分析，发现其具有一定的特点： 区域范围由经纬度坐标表示，无须投影变换； 区域坐标排列有一定规律：多数区域由直线段组成，坐标对交叉相同，即（ X1 , Y1 ）（ X2, Y1 ）（ X2 , Y2 ）（ X3 , Y2 ） （ X1 , Y1 ）； 区域坐标值多数为15的整数； 区域范围均表现为由一个主区域，并 n （ 0 , l ）个挖空区域组成； 区域坐标个数一般不多（ 1000 个），且范围有限制（小于 2500 个基本区块见矿产资源勘查区块登记管理办法）。 2建立区域矩阵 首先找出区域范围中的极值坐标（东经起、东经

4、止、北纬起、北纬止，即最大 XY 和最小 XY 值） , 以此建立一个以 l5 为基本单元的二维矩阵。 例如：主区域坐标为（ 120 . 3645 , 36 . 5130 ）（ 120 . 3900 , 36 . 5130 ）（ 120 . 3900 , 36 . 5200 ）（ 120 . 3930 , 36 . 5200 ）（ 120 . 3930 , 36 . 5300 ）（ 120 . 3700 , 36 . 5300 ）（ 120 . 3700 , 36 . 5200 ）（ 120 . 3630 , 36 . 5200 ）（ 120 . 3630 , 36 . 5145 ）（ 120

5、 . 3645 , 36 . 5145 ）（ 120 . 3645 , 36 . 5130 ），挖空区域坐标为（ 120 . 3800 , 036 . 5200 ）（ 120 . 3815 , 036 . 5200 ）（ 120 . 3815 , 036 . 5145 ）（ 120 . 3800 , 036 . 5145 ）（ 120 . 3800 , 036 . 5200 ），其极值坐标为（ 120 . 3630 , 36 . 5130 ）和（ 120 . 3930 , 36 . 5300 ），建立的 126 二维矩阵如图 1 所示。 注：以上坐标 XXX YYZZ为简写，表示 XXXYYZ

6、Z 3纵向扫描 以 Y 轴最小值（北纬起）为起点，以15为步长，用一条宽 15、无限长的矩形块对建立的栅格矩阵按坐标轴 Y 方向进行自下而上的扫描，得到交点坐标。 扫描线密度的确定 由于区域坐标不一定是15的整数倍，要保证扫描精度，必须首先确定扫描线的密度。 对坐标区域进行预扫描，得到与扫描区边界最近的距离（见图 2 ），确定扫描线个数和扫描步长，以保证扫描线与扫描区域有交点。 与区域坐标求交点 遍历区域坐标的每一条边，与扫描线求交点，则问题转化为线段与线段求交点。 利用初等代数的知识即可求解：已知线段（ xl , yl ）（ x2 , y2 ），和扫描线段 y = kv （ k ），求交点坐

7、标 Sx 其中 yl k v （ k ） y2 （对于那些 y 值没有落在扫描线区域内的线段，肯定无交点，无须进行计算，排除了 y2 = yl 的情况）。 对交点进行排序 通过以上计算，得到 n 个交点坐标（ Sx （ 1 ）, Sx （ 2 ） Sx （ n ）, Sx （ i ） Sx （ j ），且 n 为偶数（ 0 , 2 ）。因为区域坐标最小单位为 l ，而扫描线 y = kv （ k ）为小数值（ 7 . 5 , 4 . 5 , 3 . 5 ），如果与区域有交点，必成对出现。 由于区域坐标点的先后顺序不可知，因此求出的交点坐标应从小到大进行排序，从而得到一组交点对（如图 3 所示）

8、。 修改矩阵值 根据所得到一组交点对，将二维栅格矩阵所对应的行，两个交点之间的数值设置为 l ，表示被坐标区域所覆盖。 挖空区域的处理 对于挖空区域，与主区域扫描过程相同。但将二维栅格矩阵所对应的行，两个交点之间的数值设置为 0 ，表示被坐标区域所挖空。 4计算区块数 循环以上过程，可得到一个扫描完成的二维栅格数值矩阵（用数值 0 、 1 填充），可根据此数值矩阵计算区块数。 基本区块数 遍历数值矩阵。由于基本区块定义为经差 l、纬差 l ，因此首先找到起点（ 1 整数倍的最小值），步长为 l （ 4 个基本单元），考察每一点相临的 44 个单元。如果均为 l ，则基本区块数+1 ，且将该 4

9、4 个单元数值置为 0 ，表示已扫描过；如果出现 0 ，则继续遍历。 四分之一区块数 重新遍历数值矩阵。由于四分之一区块定义为经差 30、纬差 30，因此首先找到起点（ 30整数倍的最小值），步长为 30 （ 2 个基本单元），考察每一点相临的 22 个单元。如果均为 l ，则四分之一区块数+l ，且将该 22 个单元数值置为 0 ，表示已扫描过；如果出现 0 ，则继续遍历。 小区块数 最后遍历数值矩阵，步长为15（ 1 个基本单元），统计剩余数值为 1 的单元，赋值给小区块数，完成区块数计算。 5特殊情况说明 对于出现斜线组成的区域范围（见图 4 ），采取以下方法处理： 增加扫描线数：均按“

10、最短距离 2 ，扫描线 k = 7 条， y （ 1 ）= l . 5 步长 2 ”处理，以保证扫描精度； 区块计数的调整：只要斜线经过的基本单元，即视为有效单元（见图 4 ，主区域设置为 l ，挖空区域设置为 0 ） 三、算法过程描述 该算法在矿业权管理信息系统中用 VB6.0 实现。 1结构定义 Public Type POINTAPI *用于存放坐标点 x As Long y As Long End Type Public Type Polygon_truct *用于存放坐标区域 P_Num As Integer xy（ ）As POINTAPI End Type Dim Main_Po

11、ly As Polygon_Struct *主区域 Dim Sub_Poly（ ）As Polygon_Struct *挖空区域 Dim Polygon_XY （ ）As Byte *二维矩阵 2对区域坐标进行转换 计算极值坐标对区域坐标进行遍历，得到极值坐标（ Min_Px , Min-Py） （ Max _ PX , MaxPy）。 坐标转换 为增加运行速度，对区域坐标进行转换：利用function Tran_XY_Mexy(xy as POINTAPI) as POINTAPI函数转换为以秒为单位的坐标对。并将结果分别减 （ Min_Px， Min_Py ）转换为相对坐标，保证最小值为（

12、 0 , 0 ），同时调整最大值（ Max_Px， Max_Py）为相对坐标。 确定主区域和挖空区域根据坐标区域范围，分别将主区域和挖空区域放入相应的变量中（ Main_Poly 和 Sub_Poly （ ）。 3建立二维矩阵 确定二维矩阵的大小： X_num = CInt （ Max _ Px + 7 ）/ 15 ）Y _ Num = CInt （ （ Max Py 7 ）/ 15 ），并对二维矩阵 Polygon_XY （ ）赋初值 0 。 4区域扫描 使用 Sub Adjust _ Poly_XY （ y % , S _ xy As Polygon_Struct , Flag % ）过程

13、对区域进行扫描，并填充二维矩阵。 其中 y = i * 15 （ i = l Y _ Num ，步长为 15 ）； S _ xy 分别是 Main _ Poly 和 Sub _ Poly （ ）； Flag 是区域标识（ l 表示主区域， 0 表示挖空区域）。 确定扫描线数和步长 遍历 S_xy 的坐标对，得到与基本单元的最短距离 Step_V（Step_V=S_xy(i).Y mod 15），根据最短距离，确定扫描初值、扫描线数和步长（见表 l ）。 与区域坐标求交点 对于区域坐标中的线段（ xl , yl ）（ x2 , y2 ），满足以下条件进行交点计算： y + kv （k） = mi

14、n （y l , y2 ）And y + kv （k） = max （ y1 , y2 ） If xl = x2 Then Sx （ i ）= CInt （ （ xl + 7 ）/ 15 ） Else Sx （ i ）= Int （x2 一 xl ）/ （ y2 一 yl ）* （ y + kv （ k ）- yl ）+ xl ）/ 15 ） End If 对交点进行排序，填充坐标二维矩阵 对交点队列 Sx （ ）进行排序，并对二维矩阵进行赋值： For i = 0 To S _ Num - 1 Step 2 For j = Sx （ i ）To Sx （ i + l ）- l Polygon

15、XY （j， y / 15 ）= Flag Next j Next i 5计算区块数 对于基本区块，查找起点（ 1整数倍的最小值）坐标（ X _ Sta , Y _ Sta ）公式： X_Sta = （ 4 - Int （ （ Min_Px Mod 60 ）/ 15 ）Mod 4 Y _Sta = （ 4 - Int （ （ Min _ Py Mod 60 ）/ 15 ）Mod4 对于四分之一区块，查找起点（ 30整数倍的最小值）坐标（ X _ Sta , Y _ Sta ）公式： X_Sta = （ 2 - Int （ （ Min-Px Mod 30 ）/ 15 ）Mod 2 Y_sta =

16、 （ 2 - Int （ （ Min_Py Mod 30 ）/ 15 ）Mod 2 四、其他应用 本算法在矿业权管理信息系统其他方面也得到一定应用。 1.计算区块编号 描述：矿产资源勘查区块编号由其所在1：5 万图幅号、各区块行列号和四分之一区块号、小区块号共 16 位码组成（基本单位区块编码的无两位数字码为 0 ）； 基本思路：将矿产资源勘查区域用本算法建立二维矩阵，对维矩阵进行循环扫描。如果满足基本区块的要求（相临的 44 个单元均为 1 ），则根据其所在经纬度计算该基本区块编号；同理，再扫描矩阵两次，分别计算该四分之一区块编号和小区块编号； 说明：计算区块编号时由于勘查区域均由规范的基本

17、区块、四分之区块、小区块组成，计算比较简单。本算法在计算区块编号过程中的核心作用是确定勘查区域内的基本区块、四分之一区块、小区块的分布及各区块经纬度。 2基于经纬度区域的重叠检查 描述：重叠检查是指检查一个区域（指定勘查区域少是否与其他区域（其他勘查区域）存在区域重叠情况； 基本思路：将被检查的区域（指定勘查区域）用本算法建立二维矩阵；将检查区域（其他勘查区域）建立同样大小和范围的二维矩阵；将两个矩阵进行比较；得出是否重叠的结论； 说明：重叠检查比区块计算复杂如果检查区域为直角坐标，要先将坐标坐标为经纬度坐标；由于检查区域情况复杂（采矿区域、国家规划区域、地质调查区、划定矿区范围等），异形范围较多，检查速度受到一定影响；对于切割基本单元的情况（斜线穿过基本单元）不能直接比较，要冉进步细化（细化到 l1 ），算法复杂程度增加本文不做详细论述。 五、实验对比 矿业权系统早期版本的区块计算，是利用图衫平台提供的函数，将区域坐标全部铺成小区块，再通过对小区块的分析，形成四分之一区块数和基本区块数。该方法对于范围较小的勘查区域，计算速度和准确性尚可，但对大面积的区域计算速度和准确性无法保证，虽经多次优化，仍无法满足要求。 经过反复测试，使用本算法实现的区块计算现已应用到矿业权管理信息系统中，收到了较好的效果。计算速度较快，且占用内存空间不大（小于 250250 个字节）经