巧用等值線追蹤法求取地震采集設計區域邊界

鐘家均，李華科，唐雪梅

（中石化西南油氣分公司 勘探開發研究院，成都 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。