浅析钻井布局问题的数学模型(doc 18页).doc

1、钻井布局问题的数学模型摘要勘探部门在某地区找矿时，首先进行初步勘探，取几个位置钻井，取得地质资料；然后进行系统勘探，进行纵横等距的撒网式钻井。显然如果能尽可能多的在系统勘探时利用初步勘探的钻井资料，就能有效的节约费用。 在不考虑网格方向的情况下，本文首先给出了两个结论，即网格的位置由节点唯一确定，与原始矿井节点以及单位方格内的矿井映射点有相同的性质。这样就将问题等效为在单位方格内确定网格的一个节点。要解决这个问题，首先我们提出运用一般的搜索法对网格节点在单位方格内进行遍历（模型一）。通过对遍历算法进行有效的优化，大量减少了搜索的次数，进而初步计算得到了原井位最多有4个可被利用，并给出了方格节点

2、的坐标为：考虑到搜索算法的复杂度，我们给出了模型二，即在单位方格内通过确定每个矿井映射节点被利用时节点的区域，来找出方格内被这些区域覆盖次数最高的部分，显然如果将节点放在这部分内，将会有最多的点被利用，从而也就确定了节点的位置范围。运用MATLAB进行计算与判别，得到最多有4个可被利用，并求出了网格节点坐标具体的范围：当网格方向可以改变时，我们建立了模型三。考虑到判别条件是欧氏距离，可以将原题简化为一个圆形进行覆盖，圆的半径为，再用类比利用模型二进行判断，那么就能相应的找到最优规划。模型三首先进行了误差分析，根据假设的误差使用夹逼法则，然后，为了减小搜索范围，我们证明了时，最多有6个矿井可被利

3、用。对于第三问的判定算法，我们仍然根据模型三，建立假设模型四。构造出两个极端情况，此时所有矿井均可被利用。具体算法的见问题三分析步骤。 最后我们对模型四的一个假设进行了检验。虽然这个假设严格的说并不成立，但通过我们用蒙特卡罗方法进行多次模拟，发现假设成立的概率极高。综上，我们可以先用模型四进行计算，再对结果进行检验，对极少数不成立的，可以综合特殊情况进行考虑。关键词 勘探 矿井 遍历算法 蒙特卡罗 数值分析 误差分析 假设检验一 问题重述勘探部门在某地区找矿。初步勘探时期已零散地在若干位置上钻井，取得了地质资料。进入系统勘探时期后，要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面

4、钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井的地质资料，不必打这口新井。因此，应该尽量利用旧井，少打新井，以节约钻探费用。比如，钻一口新井的费用为500万元，利用旧井资料的费用为10万元，则利用一口旧井就节约费用490万元。 设平面上有n个点Pi，其坐标为(ai, bi),i=1,2,n，表示已有的n个井位。新布置的井位是一个正方形网格N的所有结点（所谓“正方形网格”是指每个格子都是正方形的网格；结点是指纵线和横线的交叉点）。假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。整个网格是可以在平面上任意移动的。若一个已知点Pi与某个网格结

5、点Xi的距离不超过给定误差（0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。为进行辅助决策，勘探部门要求我们研究如下问题： 1)假定网格的横向和纵向是固定的（比如东西向和南北向），并规定两点间的距离为其横向距离（横坐标之差绝对值）及纵向距离（纵坐标之差绝对值）的最大值。在平面上平行移动网格N，使可利用的旧井数尽可能大。试提供数值计算方法，并对下面的数值例子用计算机进行计算。 2)在欧氏距离的误差意义下，考虑网格的横向和纵向不固定（可以旋转）的情形，给出算法及计算结果。 3)如果有n口旧井，给出判定这些井均可利用的条件和算法（你可以任意选定一种距离）。 n=12个点的坐标

6、如下表所示： i123456789101112ai0.50 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50 bi2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80 二 符号说明( 原始勘测中任意一点 将旧井移至单位正方形里面时候的新井f可利用的最大钻井数n 正方形的候选位置的数目 旧井离新表格节点的距离Q（s,t） 网格的一个节点 判断旧井节点是否在表格范围内的标志位 a,b等原始节点需要转化至单位正方形内部时候需要移动的整数大小。 三 模型假设1. 初步勘探时，所

7、取的点比较分散，不存在两个点的资料被系统勘探时同一个点利用的情况。2. 系统勘探所取的点数比初步勘探时取的点多。3. 初步勘探所取的点在系统勘探的勘探范围之内。4. 假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。整个网格是可以在平面上任意移动的。5. 若一个已知点Pi与某个网格结点Xi的距离不超过给定误差（=0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。四 问题的分析及说明考虑到勘探部门在某地区找矿时是分为两步进行的，即：第一步，初步勘探，取得地质资料。第二步，系统勘探，进行“撒网式”全面钻探。那么，为了节约钻探费用，我们自然希望充分利用第一步的数据

8、，来减少勘探次数。对实际问题进行分析，可以认为初步勘探时钻井的位置可对应与二维坐标系中的点，成为初始点。系统勘探时的全面钻探可认为是二维坐标系中横坐标与纵坐标间距相等的点为节点所构成的网格。这样，问题的实质就转化为如何定位网格，使尽可能多的初始点位于网格节点的误差范围内。为了更好的表述问题，下面给出两条结论：结论一：假定网格的横向和纵向是固定的，那么只要确定其中一个节点，就可确定整个网格。证明：很显然，当网格节点中的一个点的坐标确定时，其他点的坐标可由如下公式确定：结论二：对于坐标系里的每一个点，可以定义其映射点，映射点与原来的点对网格有相同的位置关系。证明：意为对向下取整。由于网格横向和纵向

9、的单位都是1，所以当点位于网格中一个节点的误差范围内时，则其映射点必然位于的误差范围内。反之如果点不在任何一个节点的误差范围内，那么映射点也不在任何一个节点的误差范围内。五 模型的建立模型一 遍历法考虑已经给定点（矿井）以及网格的横向与纵向固定，那么问题就变为如何确定出网格的最优平移位置。由结论一，定出网格的位置只要确定出一个节点，同时，由于网格以1为单位，那么在单位方格内，必然有且只有一个节点，这样我们在四个点所构成的方块（单位方格）内对节点进行搜索，就能将网格的全部可能的平移情况进行遍历。考虑到题目中所给数据精确到0.01，那么可取0.01为两个坐标的步长进行搜索。对于搜索到的每一个坐标，

10、作为网格的节点，来确定整个网格，再计算到底有多少初始点能进入网格节点的误差范围，最后进行比较，选出容纳初始点最多的网格作为解决方案。搜索法从理论上是可行的，然而这种方法的计算量往往比较大。如上例，要搜索的节点个数为（即循环次数）：10010010000（次）而每一次循环都要对12个点依次进行判断，这必然十分复杂，当实际的数据量比较大，数据精度要求比较高时，这种直接搜索的方法计算起来就十分耗时，甚至是不可行的。为了减少循环次数，必须想办法对算法进行优化。这里我们考虑对于每一个初始点，由结论二，首先在单位方格内找到该点的映射点，这样就将所有点都映射到了单位方格内。然后进行搜索，在映射点附近搜索网格

11、节点，使映射点位于网格节点的误差范围内，再对节点构成的网格进行分析、比较，选出其中的最优。经过这样的过程，循环的次数为：1111121452（次）这样就大大减少了循环次数，优化了搜索算法，使之耗时短，使方案变得可行。 运用优化后的搜索算法，对所给数据进行计算，得到最多可以有4口井不必打而可以利用初步勘探的资料。可以利用的初始点为： (1.41,3.50) (3.37,3.51) (3.40,5.50) (8.38,4.50)网格坐标为：模型二 框图分析法考虑到遍历节点法毕竟是一种复杂度极高的算法，当数据量大的时候就会有严重缺陷，因此这里给出一种简单易行而又十分精确的算法：方框图分析法。 首先，

12、根据定理二，确定出初始点在单位方格内的坐标（即映射点），变换前后点的坐标如下：表1：坐标映射表n123456789101112初始点0.50 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50 2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80 变换点0.50 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50 0.00 0.50 0.50 0.51 0.50 0.00 0.24 0.10 0.01 0.50 0

13、.41 0.80 使用matlab软件将变换前的点画在一坐标系下，如图1：点图：1同时将映射之后的点画于一个坐标系下面得到如图2：映射点图：2将所有的旧井节点平移至单位正方形可以得到Q=（x,y）|0=x=1,0=y 在上述变换后，问题1大致等价于用一个变长为2的正方形去覆盖尽可能多的,正方形的中心就是网格的一个节点所在位置。如图-3所示：观察可得，正方形有个候选位置。若正方形左边经过，右边经过，则【,】,【】代入数据得：【,】,【】（1）根据上表1可知：=【0.5,0.41,0.00,0.37,0.40,0.72,0.72,0.43,0.57,0.38,0.98,0.50】将以上数据代入（1

14、）式得：经过将上述数据使用matlab数学软件分析可以得到，满足条件的旧井数目为4，方格节点的纵横坐标为（0.42,0.46），亦即（0.42+,0.46+）,，模型三 动态分析法网格平行移动，对网格上面的任一点，当网格移动整数个单位时，网格节点的距离不变。则现在把所有旧井节点平移至正方形：Q=(x,y)|0=x=1,0=y=1;设网格的一个节点为(s,t)，0=s,t设:ui=1否则ui=0；则以上可以归纳为：（2）等价于以某一角度为步长旋转网格，一个半径为的圆覆盖更多的旋转模式：旋转的网格图4：图4可利用当且仅当（3）最多可以利用旧井位为max,其中，ui=即 化简该算式得：设定flag为

15、节点数目的标志位，经过C语言编程以及数据的录入，得到结果flag=6，即旧井的节点数目为6个可以使用，采用欧氏距离的条件下。模型四 分类讨论法问题3的分析：假定坐标轴方向已定，考虑纵横距离。n 口井均可利用的充要条件是对 i = 1, n , 不等式组 (1) 同时有解，即如下两个不等式组有解：两组不等式，只是横坐标与纵坐标的问题，考虑证明一个即可，选择不等式组（4） , 结论.当且仅当如下二条件之一成立时，不等式组 (4) 有解：证. 当 (i) 成立时，取 s 满足,当 (ii) 成立时，取 s 满足反之，设不等式组 (4) 有解。即 (i) 成立。类似可证：不等式组 (5) 有解当且仅当

16、成立。所有旧井可利用的充要条件（网格的横纵向固定时）对纵横距离：所有旧井中的任意两口可以同时利用；对欧氏距离：所有旧井中的任意三口可以同时利用。六、模型检验 由上可以看出，无论是搜索的遍历算法还是分情况讨论，都有疏忽的地方。下面，我们就用蒙特卡罗模拟方法进行计算。 蒙特卡罗模拟方法是指运用随机抽样的方法求出事件的频率，近似的作为事件的概率的方法。具体的，对于本例，首先在一个区间内随机的取三个点作为圆心，分别判断三个圆两两相交的情况和三个圆有公共区域的情况，对这一过程做多次，统计两种情况发生的次数分别为与。其次，将作为两种情况同时出现的概率。 运用MATLAB进行模拟得到一组数据如下：4573

17、4385 0.96测试得数据足够多的情况下，我们可认为三个圆两两相交时，有96的可能三个圆有公共区域。这个概率十分大，所以我们计算时完全可以用模型四的方法进行计算，只是在得到结果后，我们需要对结果进行检验。检验的结果正确的可能性很大，少数结果不正确，只需对各角度区间按照利用个数从大到小的顺序运用模型三进行搜索即可，这样问题就全面的解决了。七、模型评价针对问题一中网格只可以平行移动的情况，我们建立了模型一搜索法和模型二方框图分析法。针对问题二中网格可以旋转的复杂情况，我们建立了模型三动态搜索和模型四。在模型一中，我们根据结论二，采用了将原钻井坐标映射为点，从而把原来网格内的点映射到一个以为顶点的

18、正方形内，然后使网格原点再在映射点附近的可行域内进行离散化搜索。这样就避免了搜索了整个网格，大大简化了计算。然而由于属于离散模型受搜索精度的影响并不可能搜索到每个映射点的可行域。在模型二中，同样用了模型一的思路将原钻井坐标映射到了一个正方形内。根据定理二，这样原来点的误差范围就转化成为了以映射点为几何中心，0.1为边长的正方形。然后通过作图的方法将映射点的可行域可视化，在图3中，观察每个映射点的可行域的范围即可确定究竟那几个钻井可被共同利用。然后就可确定最大可利用的钻井数。并求出网格原点坐标可取的具体范围。相对于模型一，它采用图形的方式给出了问题的精确解并克服了模型受精度影响的缺点。然而，随着

19、问题规模的增加，不可能通过看图的方法得出精确解，模型二并没有提供数值计算的方法。所以有一定的局限性。在模型检验中，我们利用蒙特卡罗模拟的方法对结果进行分析。但检验发现正确率为96。从而说明模型四对于求解问题仍比较可行。对于少数不成立的情况我们可以在其附近进行小范围搜索即可。参考文献1 楼顺天，MATLAB 5.x 程序设计语言，西安：西安电子科技大学，2000年4月。2 郭嗣琮,陈刚，信息科学中的软计算方法，沈阳：东北大学出版社，2001年11月。3 Patric R.J.，A fast algorithm for the maximum clique problem，Discrete App

20、lied Mathematics，120(2002)， page197-2074 卢开澄，图论及其应用，北京：清华出版社，1995年8月5姜启源，数学模型，高等教育出版社。附录：程序以及数据1. 第一问：使用Matlab软件编写的程序代码： a=0.5 1.41 3.00 3.37 3.40 4.72 5.43 7.57 8.38 8.98 9.50; b=2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80; a=0.5 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.5

21、0; plot(a,b,r) plot(a,b,*) hold on grid mkdir chapter_1999 cd chapter_1999 a=0.5 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50; b=0.00 0.50 0.50 0.51 0.50 0.00 0.24 0.10 0.01 0.50 0.41 0.80; plot(a,b,*) grid on plot(a,b,*) grid on a=0.5 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50; b=0.00 0.50 0.50 0.51 0.50 0.00 0