巧用等值线追踪法求取地震采集设计区域边界

钟家均，李华科，唐雪梅

（中石化西南油气分公司 勘探开发研究院，成都 610041）

巧用等值线追踪法求取地震采集设计区域边界

钟家均，李华科，唐雪梅

（中石化西南油气分公司 勘探开发研究院，成都 610041）

在地震采集设计时，观测系统优选论证是最基础也是最重要的一步。通常会根据不同的炮检点布设方案求取炮检点或者不同覆盖次数面元格点分布的限定多边形边界，然后求取其面积或者其内部点分布密度等参数，更多的是利用这些边界以各种底图为背景进行投影成图。如果是复杂的边界形状且没有一种反求边界的方法，则只能逐个读取点的坐标，然后手工输入坐标序列，计算、成图等操作，这个过程费时又费力。这里讨论在给定不同点集的情况下，利用矩形网格数据格点等值追踪法，求取并自动有序排列边界多边形控制点坐标，直接转换输出地震采集设计软件可用的障碍物数据格式，方便设计过程中的使用。研究证明，这种区域边界的求取方法能够满足大部分理论设计的需要，减少手工勾勒及人工读取边界坐标值的繁琐过程，实现坐标求取和应用一步到位，大大提高了边界坐标求取的准确度及采集方案设计工作的效率。

采集设计；坐标计算；障碍物；等值点追踪；区域边界

0 引言

在进行采集方案设计时，经常遇到利用排列信息反求炮检点及各种覆盖次数面元网格边界多边形的问题。一般拿到的只是针对目标区域或者构造设定的大致部署范围，需要设计者进行相关区域的资料收集整理之后，提出针对目标的多种观测系统方案（包括根据目的层深度、构造方向等，设计不同的排列长度、面元网格、规定观测方向等），进行工作量及方案优选。在这种调整和优选的过程中，会很频繁地更改观测系统模板形式或者增减部署设计工作量。而每做出一种设计方案改变，都必须快速的获得此设计方案的炮检点分布面积、分布密度，以及设计方案的满覆盖次数边界、资料边界甚至计算指定覆盖次数的边界坐标序列。

另外将这些限定边界多边形准确的投影到地形图、地质图、水文图、地理卫星图等一切可用的采集设计底图中，就必须知道这些多边形各个控制点的准确坐标。在专业设计软件没有求取限定多边形功能的情况下，也许会选择逐个的读取，然后逐个输入多边形控制点坐标，形成拐点坐标序列的文件，再进行后续的面积、长度、点分布密度等参数的计算，这种过程非常繁琐，而且很容易导致计算结果不准确。这种操作对于简单的边界形状还可以接受，如果遇到几十上百个控制点的不规则形状的限定边界形状，这种逐个点操作的方法无疑显得笨拙和相当费时。

在设计方案中，还会出现设计的炮检点因为多个不规则障碍物禁止激活操作重叠后出现炮检点分布的不规则图形，如何快速获取这种区域内的炮检点数目、该区域的面积等问题？

如果有一种利用点集反求限定边界坐标并得到其拐点序列的方法，那么无论进行部署方案如何改变，设计面积如何调整，观测系统排列和方向如何变化，都能快速准确地获得炮检点的限定边界坐标序列并形成边界控制点坐标序列文件。从而直接完成各种限定边界面积、点分布密度以及限定多边形在各种设计底图上投影成图的操作过程，能够为不同采集设计方案的优化和比较，不同观测方案的成本核算等方面大幅度提高效率。

作者借助绿山软件的地震采集设计功能，利用其设计数据，灵活选择设置不同类型的炮点、检波点、不同覆盖次数的面元网格点等，然后利用网格结点的等值线追踪法的思路计算指定类型点集的限定边界坐标并按序输出边界坐标序列，并直接转换指定格式的文本数据文件，可直接应用于面积、周长、区域内点的分布密度等参数的计算操作，从而减少手工操作的强度和提高边界连接的准确程度，大大地提高了采集设计工作效率。

1 设计过程常见的边界类型

目前常用的地震采集设计软件有GMG、Omni、KLseis、SeisWay等。这些软件支持编辑的障碍物类型一般为点状、线状和区域状，并提供了支持这些障碍物对象的外部文本格式接口。这些障碍物文本格式也较容易理解，通常涉及了障碍物的类型、对炮检点的禁止或者激活状态、控制点坐标数据个数、以及后续的坐标序列等［10－11］。如果已经获得了某个障碍物的边界各个控制点坐标序列，则很容易通过手工编辑或者编制小程序处理制作出该专业软件需要的障碍物处理格式。在不同的地震采集设计软件中，障碍物对象存档的数据格式也不同，但它们属于多边形边界中比较容易处理的简单多边形类型。图1为某工区绿山软件设计图中障碍物编辑及覆盖次数分析图。

地震采集设计中编辑处理难度较大的区域边界类型，主要为连续折线或者闭合的多边形区域。这两种区域边界类型，本身涉及的控制点坐标数据通常很多，录入编辑工作量也很大。其中连续折线类型可以当作闭合多边形的一种特例，因为只需要将闭合的连续折线（即封闭的多边形）在需要的结点处断开就可以获得该图形。

图1 障碍物编辑及覆盖次数分析图Fig．1 Obstacle editing and coverage analysis chart during seismic acquisition design

2 设计中获取区域边界的可行性

求取包围对象的边界坐标，关键是必须知道包围对象的各个部分位置坐标。采集设计过程中的一切对象（点、线、面）都有其明确的相对关系和确定的坐标位置。在地震采集观测系统设计方案中，炮检点的布设一般都是在一定的面积范围内按照一定的规则重复滚动布设。这样形成的炮检点排列方式必定在这个范围内形成规则的排列图案（图2（a）～（d））。无论是多个线性排列、正交排列或者斜交的网状排列，都可以通过一定的数据代换规则，从而转换成数据网格这种矩阵的形式。

图2 几种常见观测系统模板模拟放炮之后炮检点分布（a－d）及规则网格化处理方法（e－h）Fig．2 Distribution of SR after simulation shot of several observation system template（a－d）and regular grid data processing method（e－h）

图2的（e）－（h）表示如何将（a）－（d）的点进行移动换算，最终构造成标准的网格数据方法。在四种地震采集观测系统设计方案中，（a）方案是目前最常用的方案之一，布设之后的炮检点分布是规则的网格状分布，很容易就可以转换成按照单位增量的线号和桩号记录的网格数据，以满足后续边界追踪要求；（b）方案可以设置一个增量数组，将交错的炮线移动对齐，然后按照线号和桩号按单位1增量的形成标准网格数据（f）。同理，（c）、（d）图案也可以按照图2中与（f）类似的（g）、（h）方案进行处理。

矩形网格数据一般是以行列方式的，读取数据的时候一般是按列或者按行索引。在输出数据序列时，只有当行（列）输完时，才开始输出下一行（列）的点［4］。因此就没有办法知道已知各点之间的拓扑关系。但利用网格结点进行边界构造有大量的成熟算法和思路，可以选择实用的算法，对网格中的数据点进行区域边界构造。利用邻点追踪或者插值算法，对这些结点关系进行查找，将处于某个轨迹上或者区域边界上的点按一定的顺序连接起来［3］，并按照一定的排列顺序（顺时针或者逆时针）输出坐标数据，就可以实现不同点类型的区域圈定，进而获得输出其区域边界坐标，最终获得各种区域边界的数据文件，以供后续处理之用。

3 网格数据中区域边界求取方法

图3是一种正交观测系统的设计方案。针对需要获取图中水绿色炮点集的区域边界，求取覆盖的面积，以及其中包含的点集信息等。作者提出了网格数据中采用等值点追踪法［7－8］，该方法能够快速获取点集的圈定范围，并能够满足正交观测系统设计的绝大部分需要。

图3 具有不同类型不规则（非直角）区域炮点边界方案Fig．3 With various types of shot in irregular regions（non orthogonal）boundary scheme

3．1 网格点等值线追踪法原理及适用条件

如何转换设计数据为矩形网格数据，是等值线追踪方法求取区域边界的重点。按照图2中提供的思路，我们可以将图3中的炮点集按照线号和点号的变化规律，将所有的线和点放置到行列号增量为“1”的矩形网格中（图4）。需要求取区域边界的点集，则将其网格的结点值设置为与区域边界外不同的结点值。笔者采用的是区域外部结点置为“0”，区域内部结点全部置为“1”。区域的边界可以不是直角弯折，而可能是任意角度变化，但需要保证每个区域的点集是连续的，否则等值法追踪就会失效。

图4 布设的炮点按照标准网格数据转换后的图示Fig．4 Some text data of standard grid data graph converted from shot point array

网格结点的等值线追踪基本原理［7－8］是在网格数据中，如果用z（x，y）表示第i行第j列的高程值，Ek表示第k个需要跟踪的高程值。对于指定的高程值Ek，任意相邻两个数据点A（i1，i2），B（i2，j2），如果有（z（i1，j1）－Ek）（z（i2，j2）－Ek）＜0，则在数据点A（i1，i2），B（i2，j2）之间有一条高程值等于Ek的等值线［7－8］。等值线跟踪是在已知格网点的基础上，内插出等值线点，然后跟踪等值点，形成封闭或非封闭等值线。非封闭等值线的端点必然落在边界线上，而边界线必然是封闭的，通过跟踪边界线，再把非封闭等值线端点插入边界线，通过边界线建立等值线间的拓扑关系，追踪获得需要的区域边界。因此求不同类型点集的边界问题，可以转化为求取不同高程值的等值线问题，并且可以一次性求出不同的多个边界。

3．2 网格点等值线追踪法主要步骤

首先将设计数据进行标准网格化，对不同区域的点设置为不同的结点值，设置等值线增量间隔，并保证不同类型的标记值均落在等值线上，这种方法的网格结点数据不需要进行区域内的冗余点处理，根据不同的等值线间隔值可将多个区域连续追踪出来，并可形成专业软件需要的障碍物数据数据格式。网格点的等值线追踪算法其追踪分为以下步骤：

1）按线号（行）、点号（列）单独将炮点（检波点或者面元信息）数据转换为标准的网格数据，求取区域边界内的点置为“1”，边界外的点置为“0”（图4），并保存好网格值与实际坐标值的对应关系。标准网格中的x、y值代表原始坐标数据所在的行列值。

2）利用指针数组（列表）的方式，按照矩形等值点追踪的方法，将所有等值线值为1的线与矩形网格边的交点（用（x、y）表示）全部追踪出来（包括与边界相交的开口等值线或者内部闭合的等值线），并逐个追加存储在等值线结点数组（链表）中。考虑到追踪时结点值与等值线值相同时可能造成等值线走向混乱，因此涉及到与追踪的等值线值相同的结点，需要适当加上一个偏移量以避开这种情况。

3）对已经搜索到的等值线（边界）链表进行重复结点处理，即两个控制点之间结点按相同的X、Y增量均匀变化，则可将这些中间值标记为删除，输出时可以忽略这些结点而不会影响区域边界的形状。

4）相同的等值线高程值可能存在很多条（图5（a）），根据等值线算法原理可知，非封闭等值线端点必然落在边界上，这条不封闭等值线段必与其他不封闭等值线段构造成闭合路径，且这种等值线最多存在4条。需要新开辟一个区域边界数组（链表），将搜索到的与边界相交的等值线，将其首尾交点按照逆时针（左、下、右、上）的方向进行排列，然后调整结点顺序，追加到区域边界数组（链表）中，如果结点序列相反，则进行逆向输出。

5）考虑到存在开口等值线存在的可能性，在进行标准化网格数据转换之前，将边界进行首末行增加和首末列增加的办法对行列数据进行扩展，从而使不同点集的区域就远离边界。利用等值点法追踪区域边界时，就可以直接追踪到封闭的等值线（图5（b）），这样可以省去步骤4）。

6）逐个取出区域边界数组（链表）中的每对x（即行）、y（即列）值，换算成原始的坐标值X、Y，按专业软件规定的障碍物格式进行外部文件输出，加载即可绘图（图6）。

图5 标准网格中边界数据结点图形化显示Fig．5 Standard grid data node graphical display

图6 等值点追踪方法获得的不同类型炮点区域的叠加显示Fig．6 The different region obtained by contour point tracing method displayed with shot points

4 求取区域边界应用实例

根据上述探讨的思路，作者根据自己的地震采集设计工作内容编写了一个辅助处理工具，用于求取不同点集的区域边界，经实际应用，效果很好。图7是一个涉及各种边界比较复杂的某三维地震采集设计方案，里面涉及了拐点很多的炮检点边界、覆盖次数边界求取，以及多个禁炮区域边界。

炮点、检波点、资料边界（覆盖次数≥1）、满覆盖边界（覆盖次数≥满覆盖次数）的区域边界求取，可以分别将全部的炮点、检波点、资料边界（覆盖次数≥1的面元）、满覆盖边界（覆盖次数≥满覆盖次数的面元）矩形网格结点数据矩阵进行转换，然后进行区域内的点集进行处理，可以得到图8所示的各种边界连续离散点的网格点数据。然后利用等值线追踪法可以求出炮点、检波点、一次覆盖、满覆盖边界和三个禁炮点集的区域边界。需要注意的是，求取整体边界时，需要忽略区域内部的小空洞，即忽略三个禁炮点集。而求取内部的小空洞区域需要进行单独的转换和处理（图8）。继而完成区域的周长、面积及内部点集统计等操作。最终计算完成的区域边界如图9所示。

图7 某三维的复杂炮检点边界及多个禁炮区域Fig．7 A 3Dseismic acquisition design with complex boundaries and some regions prohibbited to shoot

图9 最终获得的某三维多种边界图案Fig．9 Finally obtained boundaries

图8 某三维边界按网格标准化处理之后的结点图Fig．8 The regular grid standardized data node graphical display of a 3Ddesign boundaries

5 结束语

在矩形网格中，利用等值线追踪法获取不同类型点集的区域边界具有较大的实用价值。该方法能够准确、快速建立边界离散点之间的拓扑关系，可应用于精确统计特定区域的周长、面积、以及区域内包含的点信息等。这种方法与专业软件融合，直接利用其设计项目文件或者转换数据计算生成区域边界障碍物数据，很大程度地减少手工操作强度或者转换数据步骤，可以在类似的设计工作中进行推广应用，并达到有助于提高其设计工作效率之目的。

［1］周培德．计算几何算法设计与分析（第四版）［M］．北京．清华大学出版社，2011．

ZHOU P D．Computational geometry：algorithm design and analysis（fourth edition）［M］．Beijing：Tsinghua University Press，2011．（In Chinese）

［2］韩同春．面向对象技术在勘察软件开发中的应用［J］．物探化探计算技术，2003，24（3）：277－278．

HAN T CH．Object oriented technology applying in software development［J］．Computing Techniques for Geophysical and Geochemical Exploration，2003，25（3）：277－278．（In Chinese）

［3］马永卓，周之平，黎明．基于离散点的任意多边形构造算法研究［J］．南昌航空大学学报：自然科学版，2012，26（4）：50－55．

MA Y ZH，ZHOU ZH P，LI M．Study on arbitrary polygon algorithm based on discrete points［J］．Journal of Nanchang University of Aeronautics：Natural Science Edition，2012，26（4）：50－55．（In Chinese）

［4］刘吉余，许洪东，王长生，等．任意面积储量计算方法研究［J］．物探化探计算技术，2003，22（1）：75－76．

LIU J Y，XU H D，WANG CH SH，et al．Research on the Calculation Method of Arbitrary Area Reserves ［J］．Computing Techniques for Geophysical and Geochemical Exploration，2003，22（1）：75－76．（InChinese）

［5］郭水平，陈锦昌．基于白色与黑色像素区域相间明显的快速边界获取的方法［J］．工程图学学报，2005（04）：104－108．

GUO SH P，CHEN J CH．A Method of quickly acquiring of edge base on white and black pixel region［J］．Journal of Engineering Graphics．2005（04）：104－108．（In Chinese）

［6］郭志宏．一种实用的等值线型数据网格化方法［J］．物探与化探，2001，25（3）：203－208．

GUO ZH H．A practical contour type data gridding method［J］．GEOPHYSICAL AND GEOCHEMICAL EXPLORATION，2001，25（3）：203 －208．（In Chinese）

［7］宋丽娟，龚晓峰，钟猛．基于网格法的等值线绘制方法［J］．现代电子技术，2005（14）：65－67．

SONG L J，GONG X F，ZHONG M．Method of drawing isoline based on grid method［J］．Modern electronic technology，2005（14）：65－67．（In Chinese）

［8］吴卫华，袁宁，陈爱斌，等．基于格网的等值线生成算法的研究［J］．实验与研究，2003（4）：28－30．

WU W H，YUAN N，CHEN AI B，et al．Study on contour generation algorithm based on grid［J］．Experiment and research，2003（4）：28－30．（In Chinese）

［9］周培德，王树武，李斌．连接不相交线段成简单多边形（链）的算法［J］．工程图学学报，2002，14（6）：522－525．

ZHOU P D，WANG SH W，LI B．Connecting non intersecting line segments into a simple polygon （chain）algorithm［J］．Journal of Engineering Graphics，2002，14（6）：522－525．（In Chinese）

Using grid contour tracking method to obtain seismic acquisition design area boundary

ZHONG Jia-jun，LI Hua-ke，TANG Xue-mei
（Exploration ＆Production Research Institute，SWPB，Chengdu 610041，China）

Observing system preferred demonstration is the most basic and important step in seismic acquisition design process．Designers usually obtained the polygon of boundary depending on the distribution of shots and receivers points or bin grid nodes with different fold number，and then calculate its area or internal point density and so on．It is more useful to map these boundary and projection them onto a variety of background image．You can only read point coordinates one by one and manually enter the coordinates sequence，which is time－consuming and laborious，for subsequent calculations，mapping and other operations when boundary shape is complex and there is no way to obtain the boundary coordinates．This article discusses how to calculate and automatically sort the point coordinates of boundary polygon in different case of a given set of points by equivalent tracing of rectangular grid data node，and then how to directly output available obstacle data format for the seismic acquisition design software to facilitate the design process．It is practical proved that method for calculating the boundary of point set can meet the needs of most of the theoretical design．The labour of hand outlining boundary and one by one reading boundary coordinate values is significantly reduced．Calculating and applications of boundary coordinates can be finished at the same time．It is greatly improved the accuracy of calculating the boundary coordinates and efficiency of acquisition program design．

acquisition design；coordinate calculating；obstacle；contour point tracing；region boundary

TP 274

A

10．3969／j．issn．1001-1749．2015．03．19

1001-1749（2015）03-0397-06

2014-07-04 改回日期：2014-10-11

钟家均（1973－），男，高级工程师，从事地震勘探部署、采集方法研究及地震采集项目技术设计等相关工作，E-mail：jinzhongyilang＠163．com。