鐘家均，李華科，唐雪梅

（中石化西南油氣分公司 勘探開發(fā)研究院，成都 610041）

巧用等值線追蹤法求取地震采集設(shè)計區(qū)域邊界

鐘家均，李華科，唐雪梅

（中石化西南油氣分公司 勘探開發(fā)研究院，成都 610041）

在地震采集設(shè)計時，觀測系統(tǒng)優(yōu)選論證是最基礎(chǔ)也是最重要的一步。通常會根據(jù)不同的炮檢點布設(shè)方案求取炮檢點或者不同覆蓋次數(shù)面元格點分布的限定多邊形邊界，然后求取其面積或者其內(nèi)部點分布密度等參數(shù)，更多的是利用這些邊界以各種底圖為背景進行投影成圖。如果是復(fù)雜的邊界形狀且沒有一種反求邊界的方法，則只能逐個讀取點的坐標，然后手工輸入坐標序列，計算、成圖等操作，這個過程費時又費力。這里討論在給定不同點集的情況下，利用矩形網(wǎng)格數(shù)據(jù)格點等值追蹤法，求取并自動有序排列邊界多邊形控制點坐標，直接轉(zhuǎn)換輸出地震采集設(shè)計軟件可用的障礙物數(shù)據(jù)格式，方便設(shè)計過程中的使用。研究證明，這種區(qū)域邊界的求取方法能夠滿足大部分理論設(shè)計的需要，減少手工勾勒及人工讀取邊界坐標值的繁瑣過程，實現(xiàn)坐標求取和應(yīng)用一步到位，大大提高了邊界坐標求取的準確度及采集方案設(shè)計工作的效率。

采集設(shè)計；坐標計算；障礙物；等值點追蹤；區(qū)域邊界

0 引言

在進行采集方案設(shè)計時，經(jīng)常遇到利用排列信息反求炮檢點及各種覆蓋次數(shù)面元網(wǎng)格邊界多邊形的問題。一般拿到的只是針對目標區(qū)域或者構(gòu)造設(shè)定的大致部署范圍，需要設(shè)計者進行相關(guān)區(qū)域的資料收集整理之后，提出針對目標的多種觀測系統(tǒng)方案（包括根據(jù)目的層深度、構(gòu)造方向等，設(shè)計不同的排列長度、面元網(wǎng)格、規(guī)定觀測方向等），進行工作量及方案優(yōu)選。在這種調(diào)整和優(yōu)選的過程中，會很頻繁地更改觀測系統(tǒng)模板形式或者增減部署設(shè)計工作量。而每做出一種設(shè)計方案改變，都必須快速的獲得此設(shè)計方案的炮檢點分布面積、分布密度，以及設(shè)計方案的滿覆蓋次數(shù)邊界、資料邊界甚至計算指定覆蓋次數(shù)的邊界坐標序列。

另外將這些限定邊界多邊形準確的投影到地形圖、地質(zhì)圖、水文圖、地理衛(wèi)星圖等一切可用的采集設(shè)計底圖中，就必須知道這些多邊形各個控制點的準確坐標。在專業(yè)設(shè)計軟件沒有求取限定多邊形功能的情況下，也許會選擇逐個的讀取，然后逐個輸入多邊形控制點坐標，形成拐點坐標序列的文件，再進行后續(xù)的面積、長度、點分布密度等參數(shù)的計算，這種過程非常繁瑣，而且很容易導(dǎo)致計算結(jié)果不準確。這種操作對于簡單的邊界形狀還可以接受，如果遇到幾十上百個控制點的不規(guī)則形狀的限定邊界形狀，這種逐個點操作的方法無疑顯得笨拙和相當費時。

在設(shè)計方案中，還會出現(xiàn)設(shè)計的炮檢點因為多個不規(guī)則障礙物禁止激活操作重疊后出現(xiàn)炮檢點分布的不規(guī)則圖形，如何快速獲取這種區(qū)域內(nèi)的炮檢點數(shù)目、該區(qū)域的面積等問題？

如果有一種利用點集反求限定邊界坐標并得到其拐點序列的方法，那么無論進行部署方案如何改變，設(shè)計面積如何調(diào)整，觀測系統(tǒng)排列和方向如何變化，都能快速準確地獲得炮檢點的限定邊界坐標序列并形成邊界控制點坐標序列文件。從而直接完成各種限定邊界面積、點分布密度以及限定多邊形在各種設(shè)計底圖上投影成圖的操作過程，能夠為不同采集設(shè)計方案的優(yōu)化和比較，不同觀測方案的成本核算等方面大幅度提高效率。

作者借助綠山軟件的地震采集設(shè)計功能，利用其設(shè)計數(shù)據(jù)，靈活選擇設(shè)置不同類型的炮點、檢波點、不同覆蓋次數(shù)的面元網(wǎng)格點等，然后利用網(wǎng)格結(jié)點的等值線追蹤法的思路計算指定類型點集的限定邊界坐標并按序輸出邊界坐標序列，并直接轉(zhuǎn)換指定格式的文本數(shù)據(jù)文件，可直接應(yīng)用于面積、周長、區(qū)域內(nèi)點的分布密度等參數(shù)的計算操作，從而減少手工操作的強度和提高邊界連接的準確程度，大大地提高了采集設(shè)計工作效率。

1 設(shè)計過程常見的邊界類型

目前常用的地震采集設(shè)計軟件有GMG、Omni、KLseis、SeisWay等。這些軟件支持編輯的障礙物類型一般為點狀、線狀和區(qū)域狀，并提供了支持這些障礙物對象的外部文本格式接口。這些障礙物文本格式也較容易理解，通常涉及了障礙物的類型、對炮檢點的禁止或者激活狀態(tài)、控制點坐標數(shù)據(jù)個數(shù)、以及后續(xù)的坐標序列等［10－11］。如果已經(jīng)獲得了某個障礙物的邊界各個控制點坐標序列，則很容易通過手工編輯或者編制小程序處理制作出該專業(yè)軟件需要的障礙物處理格式。在不同的地震采集設(shè)計軟件中，障礙物對象存檔的數(shù)據(jù)格式也不同，但它們屬于多邊形邊界中比較容易處理的簡單多邊形類型。圖1為某工區(qū)綠山軟件設(shè)計圖中障礙物編輯及覆蓋次數(shù)分析圖。

地震采集設(shè)計中編輯處理難度較大的區(qū)域邊界類型，主要為連續(xù)折線或者閉合的多邊形區(qū)域。這兩種區(qū)域邊界類型，本身涉及的控制點坐標數(shù)據(jù)通常很多，錄入編輯工作量也很大。其中連續(xù)折線類型可以當作閉合多邊形的一種特例，因為只需要將閉合的連續(xù)折線（即封閉的多邊形）在需要的結(jié)點處斷開就可以獲得該圖形。

圖1 障礙物編輯及覆蓋次數(shù)分析圖Fig．1 Obstacle editing and coverage analysis chart during seismic acquisition design

2 設(shè)計中獲取區(qū)域邊界的可行性

求取包圍對象的邊界坐標，關(guān)鍵是必須知道包圍對象的各個部分位置坐標。采集設(shè)計過程中的一切對象（點、線、面）都有其明確的相對關(guān)系和確定的坐標位置。在地震采集觀測系統(tǒng)設(shè)計方案中，炮檢點的布設(shè)一般都是在一定的面積范圍內(nèi)按照一定的規(guī)則重復(fù)滾動布設(shè)。這樣形成的炮檢點排列方式必定在這個范圍內(nèi)形成規(guī)則的排列圖案（圖2（a）～（d））。無論是多個線性排列、正交排列或者斜交的網(wǎng)狀排列，都可以通過一定的數(shù)據(jù)代換規(guī)則，從而轉(zhuǎn)換成數(shù)據(jù)網(wǎng)格這種矩陣的形式。

圖2 幾種常見觀測系統(tǒng)模板模擬放炮之后炮檢點分布（a－d）及規(guī)則網(wǎng)格化處理方法（e－h(huán)）Fig．2 Distribution of SR after simulation shot of several observation system template（a－d）and regular grid data processing method（e－h(huán)）

圖2的（e）－（h）表示如何將（a）－（d）的點進行移動換算，最終構(gòu)造成標準的網(wǎng)格數(shù)據(jù)方法。在四種地震采集觀測系統(tǒng)設(shè)計方案中，（a）方案是目前最常用的方案之一，布設(shè)之后的炮檢點分布是規(guī)則的網(wǎng)格狀分布，很容易就可以轉(zhuǎn)換成按照單位增量的線號和樁號記錄的網(wǎng)格數(shù)據(jù)，以滿足后續(xù)邊界追蹤要求；（b）方案可以設(shè)置一個增量數(shù)組，將交錯的炮線移動對齊，然后按照線號和樁號按單位1增量的形成標準網(wǎng)格數(shù)據(jù)（f）。同理，（c）、（d）圖案也可以按照圖2中與（f）類似的（g）、（h）方案進行處理。

矩形網(wǎng)格數(shù)據(jù)一般是以行列方式的，讀取數(shù)據(jù)的時候一般是按列或者按行索引。在輸出數(shù)據(jù)序列時，只有當行（列）輸完時，才開始輸出下一行（列）的點［4］。因此就沒有辦法知道已知各點之間的拓撲關(guān)系。但利用網(wǎng)格結(jié)點進行邊界構(gòu)造有大量的成熟算法和思路，可以選擇實用的算法，對網(wǎng)格中的數(shù)據(jù)點進行區(qū)域邊界構(gòu)造。利用鄰點追蹤或者插值算法，對這些結(jié)點關(guān)系進行查找，將處于某個軌跡上或者區(qū)域邊界上的點按一定的順序連接起來［3］，并按照一定的排列順序（順時針或者逆時針）輸出坐標數(shù)據(jù)，就可以實現(xiàn)不同點類型的區(qū)域圈定，進而獲得輸出其區(qū)域邊界坐標，最終獲得各種區(qū)域邊界的數(shù)據(jù)文件，以供后續(xù)處理之用。

3 網(wǎng)格數(shù)據(jù)中區(qū)域邊界求取方法

圖3是一種正交觀測系統(tǒng)的設(shè)計方案。針對需要獲取圖中水綠色炮點集的區(qū)域邊界，求取覆蓋的面積，以及其中包含的點集信息等。作者提出了網(wǎng)格數(shù)據(jù)中采用等值點追蹤法［7－8］，該方法能夠快速獲取點集的圈定范圍，并能夠滿足正交觀測系統(tǒng)設(shè)計的絕大部分需要。

圖3 具有不同類型不規(guī)則（非直角）區(qū)域炮點邊界方案Fig．3 With various types of shot in irregular regions（non orthogonal）boundary scheme

3．1 網(wǎng)格點等值線追蹤法原理及適用條件

如何轉(zhuǎn)換設(shè)計數(shù)據(jù)為矩形網(wǎng)格數(shù)據(jù)，是等值線追蹤方法求取區(qū)域邊界的重點。按照圖2中提供的思路，我們可以將圖3中的炮點集按照線號和點號的變化規(guī)律，將所有的線和點放置到行列號增量為“1”的矩形網(wǎng)格中（圖4）。需要求取區(qū)域邊界的點集，則將其網(wǎng)格的結(jié)點值設(shè)置為與區(qū)域邊界外不同的結(jié)點值。筆者采用的是區(qū)域外部結(jié)點置為“0”，區(qū)域內(nèi)部結(jié)點全部置為“1”。區(qū)域的邊界可以不是直角彎折，而可能是任意角度變化，但需要保證每個區(qū)域的點集是連續(xù)的，否則等值法追蹤就會失效。

圖4 布設(shè)的炮點按照標準網(wǎng)格數(shù)據(jù)轉(zhuǎn)換后的圖示Fig．4 Some text data of standard grid data graph converted from shot point array

網(wǎng)格結(jié)點的等值線追蹤基本原理［7－8］是在網(wǎng)格數(shù)據(jù)中，如果用z（x，y）表示第i行第j列的高程值，Ek表示第k個需要跟蹤的高程值。對于指定的高程值Ek，任意相鄰兩個數(shù)據(jù)點A（i1，i2），B（i2，j2），如果有（z（i1，j1）－Ek）（z（i2，j2）－Ek）＜0，則在數(shù)據(jù)點A（i1，i2），B（i2，j2）之間有一條高程值等于Ek的等值線［7－8］。等值線跟蹤是在已知格網(wǎng)點的基礎(chǔ)上，內(nèi)插出等值線點，然后跟蹤等值點，形成封閉或非封閉等值線。非封閉等值線的端點必然落在邊界線上，而邊界線必然是封閉的，通過跟蹤邊界線，再把非封閉等值線端點插入邊界線，通過邊界線建立等值線間的拓撲關(guān)系，追蹤獲得需要的區(qū)域邊界。因此求不同類型點集的邊界問題，可以轉(zhuǎn)化為求取不同高程值的等值線問題，并且可以一次性求出不同的多個邊界。

3．2 網(wǎng)格點等值線追蹤法主要步驟

首先將設(shè)計數(shù)據(jù)進行標準網(wǎng)格化，對不同區(qū)域的點設(shè)置為不同的結(jié)點值，設(shè)置等值線增量間隔，并保證不同類型的標記值均落在等值線上，這種方法的網(wǎng)格結(jié)點數(shù)據(jù)不需要進行區(qū)域內(nèi)的冗余點處理，根據(jù)不同的等值線間隔值可將多個區(qū)域連續(xù)追蹤出來，并可形成專業(yè)軟件需要的障礙物數(shù)據(jù)數(shù)據(jù)格式。網(wǎng)格點的等值線追蹤算法其追蹤分為以下步驟：

1）按線號（行）、點號（列）單獨將炮點（檢波點或者面元信息）數(shù)據(jù)轉(zhuǎn)換為標準的網(wǎng)格數(shù)據(jù)，求取區(qū)域邊界內(nèi)的點置為“1”，邊界外的點置為“0”（圖4），并保存好網(wǎng)格值與實際坐標值的對應(yīng)關(guān)系。標準網(wǎng)格中的x、y值代表原始坐標數(shù)據(jù)所在的行列值。

2）利用指針數(shù)組（列表）的方式，按照矩形等值點追蹤的方法，將所有等值線值為1的線與矩形網(wǎng)格邊的交點（用（x、y）表示）全部追蹤出來（包括與邊界相交的開口等值線或者內(nèi)部閉合的等值線），并逐個追加存儲在等值線結(jié)點數(shù)組（鏈表）中?？紤]到追蹤時結(jié)點值與等值線值相同時可能造成等值線走向混亂，因此涉及到與追蹤的等值線值相同的結(jié)點，需要適當加上一個偏移量以避開這種情況。

3）對已經(jīng)搜索到的等值線（邊界）鏈表進行重復(fù)結(jié)點處理，即兩個控制點之間結(jié)點按相同的X、Y增量均勻變化，則可將這些中間值標記為刪除，輸出時可以忽略這些結(jié)點而不會影響區(qū)域邊界的形狀。

4）相同的等值線高程值可能存在很多條（圖5（a）），根據(jù)等值線算法原理可知，非封閉等值線端點必然落在邊界上，這條不封閉等值線段必與其他不封閉等值線段構(gòu)造成閉合路徑，且這種等值線最多存在4條。需要新開辟一個區(qū)域邊界數(shù)組（鏈表），將搜索到的與邊界相交的等值線，將其首尾交點按照逆時針（左、下、右、上）的方向進行排列，然后調(diào)整結(jié)點順序，追加到區(qū)域邊界數(shù)組（鏈表）中，如果結(jié)點序列相反，則進行逆向輸出。

5）考慮到存在開口等值線存在的可能性，在進行標準化網(wǎng)格數(shù)據(jù)轉(zhuǎn)換之前，將邊界進行首末行增加和首末列增加的辦法對行列數(shù)據(jù)進行擴展，從而使不同點集的區(qū)域就遠離邊界。利用等值點法追蹤區(qū)域邊界時，就可以直接追蹤到封閉的等值線（圖5（b）），這樣可以省去步驟4）。

6）逐個取出區(qū)域邊界數(shù)組（鏈表）中的每對x（即行）、y（即列）值，換算成原始的坐標值X、Y，按專業(yè)軟件規(guī)定的障礙物格式進行外部文件輸出，加載即可繪圖（圖6）。

圖5 標準網(wǎng)格中邊界數(shù)據(jù)結(jié)點圖形化顯示Fig．5 Standard grid data node graphical display

圖6 等值點追蹤方法獲得的不同類型炮點區(qū)域的疊加顯示Fig．6 The different region obtained by contour point tracing method displayed with shot points

4 求取區(qū)域邊界應(yīng)用實例

根據(jù)上述探討的思路，作者根據(jù)自己的地震采集設(shè)計工作內(nèi)容編寫了一個輔助處理工具，用于求取不同點集的區(qū)域邊界，經(jīng)實際應(yīng)用，效果很好。圖7是一個涉及各種邊界比較復(fù)雜的某三維地震采集設(shè)計方案，里面涉及了拐點很多的炮檢點邊界、覆蓋次數(shù)邊界求取，以及多個禁炮區(qū)域邊界。

炮點、檢波點、資料邊界（覆蓋次數(shù)≥1）、滿覆蓋邊界（覆蓋次數(shù)≥滿覆蓋次數(shù)）的區(qū)域邊界求取，可以分別將全部的炮點、檢波點、資料邊界（覆蓋次數(shù)≥1的面元）、滿覆蓋邊界（覆蓋次數(shù)≥滿覆蓋次數(shù)的面元）矩形網(wǎng)格結(jié)點數(shù)據(jù)矩陣進行轉(zhuǎn)換，然后進行區(qū)域內(nèi)的點集進行處理，可以得到圖8所示的各種邊界連續(xù)離散點的網(wǎng)格點數(shù)據(jù)。然后利用等值線追蹤法可以求出炮點、檢波點、一次覆蓋、滿覆蓋邊界和三個禁炮點集的區(qū)域邊界。需要注意的是，求取整體邊界時，需要忽略區(qū)域內(nèi)部的小空洞，即忽略三個禁炮點集。而求取內(nèi)部的小空洞區(qū)域需要進行單獨的轉(zhuǎn)換和處理（圖8）。繼而完成區(qū)域的周長、面積及內(nèi)部點集統(tǒng)計等操作。最終計算完成的區(qū)域邊界如圖9所示。

圖7 某三維的復(fù)雜炮檢點邊界及多個禁炮區(qū)域Fig．7 A 3Dseismic acquisition design with complex boundaries and some regions prohibbited to shoot

圖9 最終獲得的某三維多種邊界圖案Fig．9 Finally obtained boundaries

圖8 某三維邊界按網(wǎng)格標準化處理之后的結(jié)點圖Fig．8 The regular grid standardized data node graphical display of a 3Ddesign boundaries

5 結(jié)束語

在矩形網(wǎng)格中，利用等值線追蹤法獲取不同類型點集的區(qū)域邊界具有較大的實用價值。該方法能夠準確、快速建立邊界離散點之間的拓撲關(guān)系，可應(yīng)用于精確統(tǒng)計特定區(qū)域的周長、面積、以及區(qū)域內(nèi)包含的點信息等。這種方法與專業(yè)軟件融合，直接利用其設(shè)計項目文件或者轉(zhuǎn)換數(shù)據(jù)計算生成區(qū)域邊界障礙物數(shù)據(jù)，很大程度地減少手工操作強度或者轉(zhuǎn)換數(shù)據(jù)步驟，可以在類似的設(shè)計工作中進行推廣應(yīng)用，并達到有助于提高其設(shè)計工作效率之目的。

［1］周培德．計算幾何算法設(shè)計與分析（第四版）［M］．北京．清華大學(xué)出版社，2011．

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

［2］韓同春．面向?qū)ο蠹夹g(shù)在勘察軟件開發(fā)中的應(yīng)用［J］．物探化探計算技術(shù)，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］馬永卓，周之平，黎明．基于離散點的任意多邊形構(gòu)造算法研究［J］．南昌航空大學(xué)學(xué)報：自然科學(xué)版，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］．物探化探計算技術(shù)，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］郭水平，陳錦昌．基于白色與黑色像素區(qū)域相間明顯的快速邊界獲取的方法［J］．工程圖學(xué)學(xué)報，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］郭志宏．一種實用的等值線型數(shù)據(jù)網(wǎng)格化方法［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］宋麗娟，龔曉峰，鐘猛．基于網(wǎng)格法的等值線繪制方法［J］．現(xiàn)代電子技術(shù)，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］吳衛(wèi)華，袁寧，陳愛斌，等．基于格網(wǎng)的等值線生成算法的研究［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］．工程圖學(xué)學(xué)報，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－），男，高級工程師，從事地震勘探部署、采集方法研究及地震采集項目技術(shù)設(shè)計等相關(guān)工作，E-mail：jinzhongyilang＠163．com。