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

矿业权管理信息系统中勘查区块数计算的实现算法一、引  言    根据国务院制定的《矿产资源区块登记管理办法》，我国对矿产资源勘查区域采用区块方式进行管理矿产资源勘查区块的划分均以 1 : 5 万图幅为基础，按经差 1’、纬差 1’划分成基本单位区块；以基本单位区块为基础，按经差 30 "、纬差 30"划分成 A、B、C、D四个四分之一区块；以四分之一区块为基础，按经差 15 "、纬差 15" 划分成 1、2、3、4四个小区块”（摘自《矿产资源勘查区块划分及编号办法》）    与之相应，《矿业权管理信息系统》的探矿权子系统和油气勘查开采子系统在数据库结构中，采用经纬度坐标保存地质调查区域、勘查区域和油气采矿区域等的空间位置，并设“基本区块数、四分之区块数、小区块数和折合基本区块”等字段记录其区块信息，以上区块数的计算由系统根据区域坐标自动完成    区块数的计算属于特定条件下的特殊应用，现在比较流行的 GIS（地理信息系统）软件，没有提供能直接获得结果的现成函数《矿业权管理信息系统》在基于 GIS 平台计算区块数没有达到预期效果的情况下，从速度、效率、成本等方面考虑，决定采用自行研制的算法解决区块数计算问题。

    二、算法思想    本算法从分析勘查区块划分的特点入手，通过建立区域坐标矩阵，将矢量坐标区域转化为以 15 " 为基本单元的栅格矩阵，并采用自下而上扫描的方式，计算出勘查区域内的栅格个数，从而得到坐标区域所覆盖的区块数     1．勘查区块管理的特点    依据《矿产资源勘查区块登记管理办法》通过对勘查区域坐标的分析，发现其具有一定的特点：    ① 区域范围由经纬度坐标表示，无须投影变换；    ② 区域坐标排列有一定规律：多数区域由直线段组成，坐标对交叉相同，即（ X1 , Y1 ）（ X2, Y1 ）（ X2 , Y2 ）（ X3 , Y2 ） … （ X1 , Y1 ）；     ③ 区域坐标值多数为15"的整数；     ④ 区域范围均表现为由一个主区域，并 n （ 0 , l … ）个挖空区域组成；     ⑤ 区域坐标个数一般不多（ < 1000 个），且范围有限制（小于 2500 个基本区块——见《矿产资源勘查区块登记管理办法》）     2．建立区域矩阵    首先找出区域范围中的极值坐标（东经起、东经止、北纬起、北纬止，即最大 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 . 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 ），建立的 12×6 二维矩阵如图 1 所示。

    注：以上坐标 XXX YYZZ为简写，表示 XXX°YY"ZZ"    3．纵向扫描    以 Y 轴最小值（北纬起）为起点，以15"为步长，用一条宽 15"、无限长的矩形块对建立的栅格矩阵按坐标轴 Y 方向进行自下而上的扫描，得到交点坐标     ① 扫描线密度的确定    由于区域坐标不一定是15"的整数倍，要保证扫描精度，必须首先确定扫描线的密度    对坐标区域进行预扫描，得到与扫描区边界最近的距离（见图 2 ），确定扫描线个数和扫描步长，以保证扫描线与扫描区域有交点    ② 与区域坐标求交点    遍历区域坐标的每一条边，与扫描线求交点，则问题转化为线段与线段求交点    利用初等代数的知识即可求解：已知线段（ xl , yl ）（ x2 , y2 ），和扫描线段 y = kv （ k ），求交点坐标 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 所示）     ④ 修改矩阵值    根据所得到一组交点对，将二维栅格矩阵所对应的行，两个交点之间的数值设置为 l ，表示被坐标区域所覆盖     ⑤ 挖空区域的处理    对于挖空区域，与主区域扫描过程相同但将二维栅格矩阵所对应的行，两个交点之间的数值设置为 0 ，表示被坐标区域所挖空     4．计算区块数    循环以上过程，可得到一个扫描完成的二维栅格数值矩阵（用数值 0 、 1 填充），可根据此数值矩阵计算区块数     ① 基本区块数    遍历数值矩阵由于基本区块定义为经差 l'、纬差 l ' ，因此首先找到起点（ 1 ' 整数倍的最小值），步长为 l ' （ 4 个基本单元），考察每一点相临的 4×4 个单元如果均为 l ，则基本区块数+1 ，且将该 4×4 个单元数值置为 0 ，表示已扫描过；如果出现 0 ，则继续遍历。

    ② 四分之一区块数    重新遍历数值矩阵由于四分之一区块定义为经差 30"、纬差 30"，因此首先找到起点（ 30"整数倍的最小值），步长为 30" （ 2 个基本单元），考察每一点相临的 2×2 个单元如果均为 l ，则四分之一区块数+l ，且将该 2×2 个单元数值置为 0 ，表示已扫描过；如果出现 0 ，则继续遍历     ③ 小区块数    最后遍历数值矩阵，步长为15"（ 1 个基本单元），统计剩余数值为 1 的单元，赋值给小区块数，完成区块数计算     5．特殊情况说明    对于出现斜线组成的区域范围（见图 4 ），采取以下方法处理：     ① 增加扫描线数：均按“最短距离 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_Poly As Polygon_Struct    *主区域     Dim Sub_Poly（ ）As Polygon_Struct   *挖空区域     Dim Polygon_XY （ ）As Byte     *二维矩阵     2．对区域坐标进行转换     ① 计算极值坐标对区域坐标进行遍历，得到极值坐标（ Min_Px , Min-Py） （ Max _ PX , Max——Py）。

    ② 坐标转换    为增加运行速度，对区域坐标进行转换：利用function Tran_XY_Mexy(xy as POINTAPI) as POINTAPI函数转换为以秒为单位的坐标对并将结果分别减 （ Min_Px， Min_Py ）转换为相对坐标，保证最小值为（ 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 % ）过程对区域进行扫描，并填充二维矩阵    其中 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）＞ = min （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     Po。