CN103969683B - Method for picking position faces in batched mode based on constraint in three-dimensional seismic interpretation - Google Patents

Abstract

The invention discloses a method for picking position faces in a batched mode based on the constraint in three-dimensional seismic interpretation. An upper position face and a lower position face which are picked manually or in other modes serve as constraint faces, known conditions provided by the constraint position faces are obtained, the interaction relation between middle position faces to be picked is combined, and all position faces between the constraint position faces are picked synchronously. The new method for automatically picking a plurality of position faces with constraint information added is provided, the information of the upper constraint position face and the information of the lower constraint position face are combined in the method, and auxiliary information of the position faces to be picked is provided. A seed seismic channel spreading method is adopted on sample points which are obtained through filtering, and the accuracy and the time efficiency can be improved.

Description

A kind of method that in 3-D seismic interpretation, the batch based on constraint picks up layer plane

Technical field

The invention belongs to computerized algorithm, Computer Image Processing and Review of geologic model building techniques field, it is a kind of based on the overall situation Optimum, spreads search carried out while automatic Picking multiple layers as auxiliary information with reference to seed seismic channel using restraint layer plane The method of plane.

Background technology

Seismic exploration technique is the important technical of oil-gas exploration.The explanation efficiency of seismic interpretation and the essence of explanation results Degree directly influences exploration benefit, therefore earthquake interpretation algorithms is studied, and improves the computational efficiency of seismic interpretation algorithm, changes The Explanation Accuracy of kind seismic interpretation algorithm, can effectively shorten seismic survey period, save exploration cost.

Form reflection due to meeting geological interface (layer position) in artificial earthquake wave direction underground propagation, shake hence with seismic wave The information change such as width, frequency can reflect subsurface geologic structures.In general seismic profile, amplitude is maximum or smallest point is considered It is the most believable part of interpretation horizon, from the back wave at same interface, on adjacent common depth point, the time is very close, pole Property same phase is consistent.The line in earthquake record Shang Ge road vibration phase identical extreme value (being commonly called as crest or trough) is referred to as homophase Axle, same reflection waveform, the lineups continuing phase place are parallel to each other.It can be said that knowing the homophase in the reflecting layer on seismic profile Axle is known that its corresponding geologic horizon.

On time section, the Wave data near the upper and lower both sides of reflecting layer lineups just reflects the geology letter of geologic horizon Breath.If the geological information wondering geologic horizon, it is necessary to the data comprising this information is extracted, is analyzed further And process.Obtain this data, just extract the waveform feature parameter near the upper and lower both sides of homophase track shaft, the method adopting now Exactly find the track wave crest point (or trough point) of this lineups.When following the trail of lineups, mainly follow following criterion：

(1), with phasic property criterion (continuity criterion)：The time that the back wave at same interface occurs on adjacent record is near As, lineups should be smooth, and extends longer；

(2), energy criterion：Because back wave has certain energy, he is typically larger than the energy of jamming pattern, and reflection is shaken Width is bigger, more can judge the presence reflected.It is thus determined that the extreme value of record is most important；

(3), waveform similarity criterion：The wave character being reflected on adjacent record at same interface is similar.For example, The apparent cycle (or video frequency) of back wave, number of phases, the strong or weak relation of each phase place are similar.

Currently, increasingly increased by the geological data amount obtaining of seismic acquisition, traditional artificial horizon picking side Method not only inefficiency, and easily affected by human factor and noise.In order to overcome the weakness of Traditional Man pickup, perhaps The automatic pick method of many seismic horizons arises at the historic moment, such as neural net method, inclination angle sector method, mode identification method, company Matching algorithm, cross-correlation method etc..Since the seventies in last century, lot of domestic and international scholar utilizes waveform similarity principle research Seismic horizon automatic pick method.Artificial neural network is a kind of new technology developed in recent years, it have self study, from The function such as adapt to and associate, can be used to solve the problems, such as lineups automatic Picking.Seismic profile has multiform expression, permissible Seismic profile is regarded as a two-dimentional digital figure, equally being capable of accurate pickup earthquake layer using image processing method Position, such as edge detection method.

The layer plane automatic pick method that presently, there are is all based on the pickup of individual layer of plane mostly, and these methods are being picked up Take very high in the efficiency of individual layer of plane and precision, but shortcoming is cannot to obtain multiple layer of plane simultaneously.The present invention is to base Carry out improving it is proposed that oneself of multiple layer of plane is picked up in a kind of synchronization in the algorithm that seed point spreads method PICKING UP THE EVENTS AUTOMATICALLY Dynamic picking algorithm.

Content of the invention

The problem to be solved in the present invention is：The tomography of conventional three-dimensional seismic interpretation extracts and relies on manually, and artificial demand is big, when Between consume big, therefore efficiency is low.Existing extraction method is limited to the extraction of individual layer of plane it is impossible to obtain multiple layers simultaneously Plane.

The technical scheme is that：A kind of method that in 3-D seismic interpretation, the batch based on constraint picks up layer plane, Identify that upper and lower two smooth layer plane, as containment surfaces, using the known conditions of containment surfaces, have phase between layer plane The mutually relation of constraint, all layers of plane between synchronous pickup containment surfaces, comprise the following steps：

Step 1：Read three-dimensional data volume data, obtain extreme point with interpolation method, and the extreme point obtaining is filtered Process, obtain sample point；

Step 2：According to the layering of geological type, identify that upper and lower two smooth layer plane, as containment surfaces, obtain constraint The information of each sample point on face, as the constraint information of pickup layer plane, described constraint information includes each sample in containment surfaces The obliquity information of this point and amplitude information；

Step 3：With two containment surfaces of step 2 as boundary, seed seismic channel is set, is searched based on " seed point spreads method " Rope, progressively extends pickup sample point according to pickup standard, obtains layer plane to be picked up, described pickup standard refer to sample point it Between similarity set threshold range in：

Primary earthquake road, referred to as seed seismic channel are set first, neighbor search is carried out respectively to its adjacent seismic channel, to plant Based on sample point on sub- seismic channel, judge whether to pick up the sample point of adjacent seismic channel, the binding characteristic condition bag of search Include pickup type, amplitude similarity, depth similarity and inclination angle constraint, for two current sample points, it is right that pickup type refers to Two sample points are maximum or the judgement of minimizing uniformity, are unanimously then picked up；Amplitude similarity refers to two samples Whether the amplitude similarity of point in the threshold range setting, then picks up；Depth similarity refers to the depth value of two sample points Whether difference in the threshold range setting, then picks up；Inclination angle constraint refers to adjacent upper in sample point is for 3D data volume During lower two sample points, whether the inclination angle similarity between two sample points in the threshold range setting, then picks up；

Qualified sample point, quilt are picked up on the adjacent seismic channel of seed seismic channel according to above-mentioned binding characteristic condition The seismic channel that the sample point of pickup is located becomes new seed seismic channel, to new seed earthquake to proceeding neighbor search, If all seismic channels in the neighborhood of a seed seismic channel are all searched, this seed seismic channel is labeled as Search, the search after being no longer participate in, when all of seismic channel is all searched once, iteration ends；Pickup in sample point During, sample point divides automatically according to binding characteristic condition, constitutes different layer planes；

Step 4：The step 3 layer plane that obtain of pickup is smoothed and plugs a hole process so as to get layer plane Completely, complete the batch pickup of all layers of plane between containment surfaces.

Step 1 is specially：For one piece of given 3D data volume, quadratic polynomial interpolation is adopted to each seismic channel Method, obtains an extreme point according to three serial number points on seismic channel every time, i.e. crest or trough is adopted to the extreme point obtaining Use digital averaging filtering method, when occur two or more amplitude be continuously on the occasion of or be continuously negative value extreme point when, this is several Individual continuous extreme point synthesizes a synthetic point, and the coordinate of synthetic point is the coordinate mean value of the extreme point being synthesized, amplitude It is worth the mean value of amplitude for the extreme point being synthesized, after all extreme points are completed with filtering, obtain the sample point of each seismic channel.

The search Ji Yu " seed spread method " for the step 3 includes spreading operation and pick-up operation, spreads operation first from current The seed seismic channel processing is traveled through to the adjacent seismic channel not being searched；Pick-up operation be traversal current certain During one adjacent seismic channel, according to the pickup standard setting, sample point pick-up operation is carried out to described adjacent seismic channel, step 3 has Body is：

1) spread operation

Spread the operation traversal BFS method based on figure for the flow process, construct a figure corresponding with 3D data volume： Each of 3D data volume seismic channel is configured to one of figure point, is any two adjacent seismic channel Trace_iAnd Trace_j Construction a line E_ij, from seed seismic channel corresponding summit, S proceeds by breadth first traversal, is adjacent until having traveled through All summits, that is, traveled through all seismic channels being adjacent；

Wherein, three Mark Arrays BF, BNF and T are arranged to seismic channel, T records the search condition of seismic channel, if current position Shake road is marked as 1, representing that this seismic channel has been searched, if being marked as 0, representing that this seismic channel is not yet searched；BF Record current seed seismic channel, referred to as boundary trace, if seismic channel carries out neighbor search as seed seismic channel is designated as 1, no It is then 0；The BNF new seismic channel that searches from boundary trace of record, BNF be marked as this seismic channel of 1 expression will as under The boundary trace once searched for；Spread operation to store the value of BNF in BF first, represent and the seismic channel becoming new boundary trace is set Be set to current border road, carry out spreading search by current border road, i.e. BFS, and by all values of BNF empty for 0, in follow-up pick-up operation, recording again new will become the seismic channel of boundary trace；

2) pick-up operation

Current border road is referred to as source road, says from source and the seismic channel traverse is sent to its adjacent direction is referred to as target Road, pick-up operation finds, in target track, the sample point meeting pickup standard, has processed a target track just by this target every time On road, picked sample point and corresponding source road sample point add same layer plane, and it is special that described correspondence refers to meet the constraint Levy condition, pick-up operation completes the automatic division of sample point according to binding characteristic condition, obtains each layer of plane；

When pick-up operation is carried out, from Yuan Dao, if the T of the target track traversing is labeled as 1, continue with next Individual adjacent seismic channel；If the T of target track is labeled as 0, then target track is picked up operate, picks up meet the constraint feature bar The sample point of part, this point is added in corresponding layer plane；Setting threshold value t, when having been added to a layer plane in a target track The percentage of the number of sample point exceed threshold value t, then arranging the BNF mark of this target track and T flag is 1, represents this target Road is searched to be processed and became the boundary trace next time searched for；After road corresponding all target track search in source terminate, Complete a pick-up operation, now, in BNF array, recorded new boundary trace.

In pick-up operation, design energy function is used for describing pickup standard, specially：

Choose a seismic channel as source road first between containment surfaces, from Yuan Dao, judge that target track is about No searched, if source road is Trace₁, the not searched target track of it some is Trace₂, n₁And n₂Represent respectively this two The sample point quantity containing on individual seismic channel, is that the sample point in source road and target track puts on sequence number, respectively S_1l,S_2g, wherein L=1,2 ..., n₁, g=1,2 ..., n₂If m represents pickup mode, pickup mode refers to the sample point picking up to set, m={ m₁, m₂,...,m_k, m_pRepresent a sample point pair, m_p=(S_1l,S_2g),l∈[1,n₁],g∈[1,n₂], p ∈ [1, k],

Setting pickup constraints：(1) layer plane can not intersect；(2) phase polarity of layer plane is consistent, design energy Function f describes pickup mode, and energy function f illustrates the binding character between the correlation between sample point and layer plane：

F (m)=exp (- U (m)) (1)

U (m) is a linear combination of the restrictive function between relevance function and layer plane between sample point：

U (m)=ω_cC(m)+ω_rR(m) (2)

Function C (m) represent sample point in source road and target track between total amplitude coefficient correlation, function R (m) represents source Total inclination angle coefficient correlation of all adjacent sample points, ω on road_cWith ω_rThe weights of C (m) and R (m) respectively, described weights pass through Experimental debugging obtains；

Described amplitude coefficient correlation is：

For a sample point S (x, y, z), its place seismic channel is selected centered on S (x, y, z) each K up and down Sample point, to be described as the wave character of S (x, y, z) using this 2K+1 sample point { S (x, y, z') | | z'-z |≤K }, that is, this 2K+1 sample point only has depth value z different, and x and y coordinates are identical；Normalized crosscorrelation technology is selected to describe sample point local The similitude of waveform, function E (x, y, z) represents the amplitude of sample point S (x, y, z), for two sample point s₁(x₁,y₁, z₁), s₂(x₂,y₂,z₂), the normalized-cross-correlation function between them is：

For a pickup mode m={ m₁,m₂,...,m_k, wherein there is k to sample point pair, by the phase of each sample point pair Close coefficient to add up, obtain total sample point amplitude coefficient correlation in pickup mode m in current source road and target track, that is,

Described inclination angle coefficient correlation is：

With the angle information between two inclination angles as judging both similarity degrees, if dip is (S_1l,S_2g) and dip (S_1l′, S_2g′) it is respectively the obliquity information of adjacent two sample points on same seismic channel, wherein l=1,2 ..., n₁, g=1, 2,…,n₂, l '=1,2 ..., n₁, g '=1,2 ..., n₂, will be similar as inclination angle for the cosine value of the angle between this two inclination angles Degree, on the road of source, the inclination angle similarity factor of two adjacent sample points is defined as following form：

Mode m={ m thus can be picked up₁,m₂,...,m_kTotal inclination angle coefficient correlation：

Wherein R (S_1,u,S_1,1) between containment surfaces are intersected with source road in expression sample point and first sample point in source road Inclination angle similarity factor, R (S_1,k,S_1,d) represent the sample that on the road of source, last sample point and lower containment surfaces are intersected with source road The similarity factor at the inclination angle between point；

After the solution to amplitude coefficient correlation and inclination angle coefficient correlation, obtain different pickup modes between seismic channel Corresponding energy function, finds an optimum pickup modeMake energy functionFor minimum, that is,To meet Following form：

Wherein, M represents the set of all possible pickup mode, for optimum pickup modeIncluded in all samples This point is right, each sample point centering, and the sample point of target is assigned in the layer plane belonging to source road sample point, thus simultaneously Sample point in multiple layers of plane of pickup, realizes batch and picks up；

In conjunction with formula (1) to formula (7), target track is picked up operate, after having executed pick-up operation, an iteration Process terminates, and judges whether the road also not searched according to T flag, if there are the road not searched, then continues to change next time For process, otherwise algorithm terminates.

Further, in pick-up operation, a layer plane label is carried out to sample point, is easy to the differentiation of layer plane, for sample The layer plane label set-up mode of this point is as follows：

First, it is that the currently affiliated layer plane of all of sample point demarcation is numbered 0, then from primary earthquake road, Process first seismic channel adjacent with it, after obtaining the optimum pickup mode between them, split from first sample point Begin Distribution Layer plane label one by one, and layer plane label from the beginning of 1, that is, gives the sample point in optimum pickup mode corresponding to distributing Layer plane label, now just contains two sample points in each layer of plane, and for the not sample in optimum pickup mode Point then still retaining layer plane be numbered 0；Spread in operation and pick-up operation in follow-up, in each optimum pickup mode Each sample point pair, if the sample point belonging on the road of source has had a layer plane label, then corresponding with it belong to target Sample point on road just puts on the layer plane label the same with it, if the layer plane label of two sample points is all 0, according to work as The front layer plane number having existed, increases a layer plane label newly, and described two layers of plane label is all that 0 sample point just becomes For two sample points of new layer plane, the layer plane label updating of two sample points is new layer plane label simultaneously.

In step 4, smoothing processing takes mean value method, and the depth to each sample point on the layer plane of pickup, by institute State the mean value that depth value is the 8- neighbour of sample point and sample point itself totally 9 points；If a certain seismic channel not by Pick up, but its adjacent seismic channel has more than 6 to be all searched, then carry out process of plugging a hole, this is not picked up The value of the seismic channel got is disposed adjacent to the mean value of seismic channel sample point, through smoothing processing and process of plugging a hole, obtains Last layer plane.

In step 2, the layer bit-plane positions of pickup according to required for explanation personnel, by artificial pickup or other existing sides Formula corresponding pickup containment surfaces.

Because in seismic profile, amplitude is maximum or smallest point is considered as the most believable part of interpretation horizon, the present invention exists The just sample point at the crest just for seismic wave and trough when geologic data is processed, that is, extreme point.To extreme value Point reduces Local Extremum to the interference working afterwards after being filtered.Increase the information of upper and lower containment surfaces, can improve The precision of pickup layer plane is so that the layer plane of pickup is suitable for whole geological environment.Adopt during processing each seismic channel The pickup mode of optimum, with pickup mode is converted into searching Optimal solution problem, is obtained from global sense.By spreading and picking up Take two operations, all seismic channels in 3-d seismic data set are picked up.Finally the layer plane tentatively picking up is entered Row smoothing processing and plugging a hole, in conjunction with the higher layer plane of precision achieved above.

In addition, the seismic image of layer position region has certain similitude, using this characteristic, carrying out the knowledge of layer position Just hunting zone can be limited in the range of very little when other, so both can improve search efficiency and can also avoid layer " string layer " problem during search of position.

The automatic Picking layer plane method of the present invention is a little compared to existing method：Constrained based on upper and lower two In the case of face, each sample point in containment surfaces contains the geologic feature itself having, such as obliquity information and move towards letter Breath.These by with prior art for example OpendTect instrument pickup individual restraint layer plane there is higher confidence level, because Information on this restraint layer plane has relative directive significance.From formula (6), containment surfaces are in the process of pickup sample point In provide similitude aspect constraint so that picked sample point must be to the similar system of the sample point in these containment surfaces Number is in threshold range；Secondly, in pick process, every time two adjacent seismic channels are processed, pick up all , as candidate, to set, each sample point is to respectively for the sample point on one of pickup mode as two seismic channels for mode Any is respectively taken to constitute by two seismic channels.In processing procedure, an energy function is set, energy is solved to each pickup mode Functional value, the pickup mode obtaining minimum energy value can be considered as optimal solution.Carry out spreading method pickup according to said process Sample point on all seismic channels, realizes the automatic division of sample point, the batch pickup of complete layer plane.

Brief description

Fig. 1 is the flow chart of the present invention.

Fig. 2 is Dip countion form on inline direction for the sample point in the present invention.

Fig. 3 is the schematic flow sheet of pick process in the present invention.

Specific embodiment

Due to the diversity of geological type, there are some geological layering can be with instrument as OpendTect come very accurate True identifies and extracts.Therefore, it can the two smooth layer planes up and down that on this basis, will identify that as constraint Face, the layer plane having multiple complexity between this two containment surfaces is it is also possible to contain tomography.Containment surfaces can provide many letters Breath, the such as trend of current layer position, obliquity information of each sample point etc., adjacent layer be can be seen that on amplitude section figure In the trend of plane, and seismic channel identical with each, there are similitude, the therefore inclination angle of containment surfaces in the inclination angle of intersecting sample point Information can be used as the tutorial message of layer plane to be identified.These contents can be formed pickup intermediate layer position pact to be met Bundle condition, thus so that the layer plane of pickup is more accurate.Fig. 1 gives the main flow of the present invention.

In conjunction with the drawings and specific embodiments, the present invention is described in further detail：

Step 1：Read three-dimensional data volume data, obtain extreme point as pending sample point with interpolation method, and to To extreme point be filtered process.

For one piece of given 3D data volume, according to the phase type of lineups, including crest (MAX) and trough (MIN) two kinds, will taking Find Max operation as a example be introduced below, that is, determine crest, the reality of another operation (trough) Now it is similar to.The interpolation method adopting in the present invention is quadratic polynomial interpolation method, thus can be obtained according to three sample points every time One extreme point.

The extreme point that obtained according to above-mentioned interpolation method it is possible that crest be negative value and trough be on the occasion of result, this When need result is filtered.The filtering method adopting in the present invention is digital averaging filtering method, that is, several continuous in appearance Be on the occasion of or be continuously negative value extreme point when, these sample points are synthesized a sample point, the position of this sample point is The mean place of these points, the amplitude of sample point is the mean amplitude of tide of these points.After after filtering, pseudo- extreme value can be eliminated The interference of point.

Step 2：Read upper and lower containment surfaces, obtain the information of each sample point in containment surfaces, as the pact of pickup layer plane Bundle information.

Containment surfaces are the known layer plane with high confidence level by OpendTect instrument pickup, this two constraints On face, all of sample point provides the information currently itself being comprised, can be according to each sample point and its adjacent sample Put and to obtain some features of itself, such as obliquity information, layer plane trend etc..

The obliquity information taken in the present invention is：If S (x, y, z) represents the sample point on seismic channel, dip_inline, dip_crosolineIt is respectively the inclination value on this sample point inline and crossline direction： Wherein inline direction and crossline direction are two line directions in 3-d seismic data set, its Middle inline represents main profile direction, and crossline represents cross-track direction；△ t represents current sample point and adjacent earthquake Depth difference between sample point on same layer plane is belonged on road, △ x represents current sample point and belongs to same on adjacent seismic channel Alternate position spike on inline line direction between sample point on one layer of plane, △ y represents current sample point and adjacent seismic channel On belong to alternate position spike on the crossline line direction between sample point on same layer plane.Referred to herein as inclination value dip_inlineAnd dip_crslineWith dip (S described later_1l,S_2g) implication is consistent, calculation is consistent.

Step 3：Setting seed seismic channel, progressively extends the suitable sample of pickup to spread method to its eight adjacent seismic channels Point, obtains layer plane to be picked up.

The algorithm picking up layer plane in the present invention is a process being constantly iterated, as shown in figure 3, in figure round dot table Show seismic-data traces, the round dot of in figure mark " S " represents and is presently processing seismic channel, Ji Yuan road, and inl represents inline Line direction, crl represents crossline line direction.Scan for its adjacent road from current seismic road, adjacent seismic channel can From 4- neighbour or 8- neighbour.The present invention adopts 8- neighbor search.

According to pickup type, amplitude similarity, depth similarity, wave-form similarity, the feature such as inclination angle constraint, it is being adjacent to Shake picks up suitable sample point on road.For two current sample points, pick up type and refer to two sample points for maximum also It is that minimizing uniformity judges, be unanimously then picked up；Whether the amplitude similarity that amplitude similarity refers to two sample points exists In the threshold range setting, then pick up；Whether depth similarity refers to the depth value difference of two sample points in the threshold value setting In the range of, then pick up；When inclination angle constraint refers to adjacent two sample points up and down in sample point is for 3D data volume, two samples Whether the inclination angle similarity between this point in the threshold range setting, then picks up.

Search procedure includes spreading operation and pick-up operation, and spreading is to neighboring track from the seismic channel being presently processing The process being traveled through, pick-up operation is the process of the sample point that search meets pickup standard in neighboring track, picked sample The seismic channel that this point is located will become new seed seismic channel.The round dot of " F " is marked to represent new seed seismic channel in Fig. 3, he Represent current search border.If all roads in the neighborhood of a seed seismic channel are all searched, then this Individual seed seismic channel will be labeled as searching for, and is no longer participate in next step search, when all of seismic channel is all searched once, Algorithm terminates.Here the search Ji Yu " seed spread method " includes spreading operation and pick-up operation, spreads operation first from current The seed seismic channel processing is traveled through to the adjacent seismic channel not being searched；Pick-up operation be traversal current certain During one adjacent seismic channel, according to the pickup standard setting, sample point pick-up operation is carried out to described adjacent seismic channel, specifically retouches State as follows：

(1) spread operation

Spreading operation is the process that the road not being searched in the road searched for is to its field is traveled through, its time Go through BFS (Breadth First Search, the BFS) method similar to figure for the flow process.Construction one and 3-D seismics The corresponding figure of data volume：Each of data volume seismic channel is configured to one of figure point, is the adjacent seismic channel of any two Trace_iAnd Trace_jConstruction a line E_ij.From seed seismic channel corresponding summit, S proceeds by breadth first traversal, until time Go through all summits being adjacent, that is, travel through all seismic channels being adjacent.

Spreading operation is based on BFS, for this reason, three Mark Arrays of setting, to mark whether respectively as side Boundary, whether it is new border and whether is searched, by array BF, BNF, T are representing.Wherein T records searching of all seismic channels Strand state, if current seismic road is marked as 1, representing that this seismic channel has been searched, if being marked as 0, representing this earthquake Road is not yet searched；BF records current seed seismic channel, referred to as boundary trace, if seismic channel is carried out closely as seed seismic channel Adjacent search is designated as 1, otherwise for 0；The new seismic channel that BNF record searches from boundary trace, BNF is marked as 1 expression should Seismic channel is using as the boundary trace searched for next time.Spread operation first the value of BNF to be stored in BF, that is, representing to become new The road of border seismic channel is formally set to border seismic channel, thus can carry out spreading search from this seismic channel.Then by BNF's All values empty as 0, thus can be in pick-up operation next, and recording new will become the seismic channel of boundary trace.

(2) pick-up operation

The work of pick-up operation is to find, in target track, the sample point meeting feature constraint condition, has processed one every time Picked sample point on this road will be added in corresponding layer plane for target track.It is described below for convenience, will The seed seismic channel being presently processing is referred to as source road, says from source and is sent to the seismic channel newly traverse its adjacent direction and claims For target track.After spreading operation, can get current desired seismic channel to be processed, Ji Yuan road.From Yuan Dao, such as The T of the target track that fruit traverses is labeled as 1, then continue with next adjacent seismic channel；If the T of target track is labeled as 0, that Target track is picked up operate, pickup meets the sample point of feature constraint condition, and this point is added to corresponding layer In plane set.One threshold value t of setting, when having been added to a layer percentage for the number of the sample point of plane set in a seismic channel Ratio exceedes threshold value t, then arranging the BNF mark of this seismic channel and T flag is 1, represents that this target track is searched and processed simultaneously And become the boundary trace next time searched for；After road corresponding all target track search in source terminate, complete a pick-up operation, this When, have recorded new boundary trace in BNF array.

Core concept in pick-up operation is used for describing pickup standard for design energy function, and specific design scheme is：

Choose a seismic channel first as initial boundary road, referred to as source road, from Yuan Dao, judge about adjacently Whether shake road has not searched seismic channel, if currently processed source road is Trace₁, the not searched target of it some Seismic channel is Trace₂, n₁And n₂The sample point quantity containing is represented respectively on this two seismic channels.For the source being presently processing Sample point on road and target track puts on sequence number, respectively S_1l,S_2g, wherein l=1,2 ..., n₁, g=1,2 ..., n₂.If m table Show pickup mode, the sample point being expressed as current pickup is to set, i.e. m={ m₁,m₂,...,m_k, represent in current pickup mode Comprise k to sample point pair, wherein m_p=(S_1l,S_2g),l∈[1,n₁],g∈[1,n₂], p ∈ [1, k], represent a sample point Right, the now S on the road of source_1lWith the S in target track_2gIt is allocated and belong to same layer plane.Setting pickup constraints：(1) layer position Face can not intersect；(2) phase polarity of layer plane is consistent.According to these constraintss, can be significantly reduced those improper Pickup mode.

Design energy function f describes pickup mode, energy function f include correlation between sample point and layer plane it Between binding character, that is,：

F (m)=exp (- U (m)) (1)

Wherein U (m) is a linear group of the restrictive function between relevance function and layer plane between sample point Close：

U (m)=ω_cC(m)+ω_rR(m) (2)

Wherein function C (m) and R (m) show respectively sample point in source road and target track between total amplitude phase relation Total inclination angle coefficient correlation of all adjacent sample points, ω on Shuo Heyuan road_cWith ω_rIt is respectively respective weights, all weights lead to Cross experimental debugging to obtain.Because in above-mentioned formula, the initial value unit of each data and scope are all different, so counting for convenience Calculate and the more preferable effect that embodies has to these data are normalized；

Described amplitude coefficient correlation is：

For a sample point S (x, y, z), its place seismic channel is selected centered on S (x, y, z) each K up and down Sample point, to be described as the wave character of S (x, y, z) using this 2K+1 sample point { S (x, y, z') | | z'-z |≤K }, that is, this 2K+1 sample point only has depth value Z different, and x and y coordinates are identical.Normalized crosscorrelation technology is selected to describe sample point local The similitude of waveform, E (x, y, z) represents the amplitude of sample point, for two sample point s₁(x₁,y₁,z₁), s₂(x₂,y₂,z₂), Normalized-cross-correlation function between them is：

For a pickup mode m={ m₁,m₂,...,m_k, wherein there is k to sample point pair, by the phase of each sample point pair Close coefficient to add up total sample point amplitude coefficient correlation in a pickup mode obtaining in current source road and target track, that is,

Described inclination angle coefficient correlation is：

With the angle information between two inclination angles as judging both similarity degrees.If dip is (S_1l,S_2g) and dip (S_1l′, S_2g′) wherein l=1,2 ..., n₁, g=1,2 ..., n₂, l '=1,2 ..., n₁, g '=1,2 ..., n₂, respectively same The obliquity information of two adjacent sample points on seismic channel, then using the cosine value of the angle between this two inclination angles as inclination angle phase Like spending, on the road of source, the inclination angle similarity factor of two adjacent sample points is defined as following form：

A pickup mode m={ m thus can be obtained₁,m₂,...,m_kTotal inclination angle similarity factor：

S in formula_1,1Represent Trace between containment surfaces₁On first sample point, the k in formula just refers to pickup mode Middle k is to k sample point on sample point Zhong Yuan road.U refers to up, point on the road of source for the containment surfaces in expression, and d represents under down about Point on Shu Mianyuan road.R(S_1,u,S_1,1) represent the sample point that on the road of source, first sample point and upper containment surfaces are intersected with source road Between inclination angle similarity factor, R (S_1,k,S_1,d) represent sample in last sample point and lower containment surfaces and source road on the road of source The similarity factor at the inclination angle between this point.Seismic channel is perpendicular to a virtual noun on ground, is by a series of x, y-coordinate Identical, the different sample point of z coordinate is constituted, contain the sample point between containment surfaces and outside sample point.And layer plane (comprising containment surfaces) is the similar curved surface parallel to horizontal plane, so seismic channel and each layer of plane have intersecting sample point, Here S_1,uAnd S_1,dIt is the sample point that upper and lower containment surfaces are intersected with source road.

After the solution to amplitude coefficient correlation and inclination angle coefficient correlation, obtain each energy corresponding to pickup mode Flow function, finds an optimum pickup modeMake energy functionFor minimum, that is,To meet following form：

Wherein, M represents the set of all possible pickup mode.

In conjunction with formula (1) to formula (7), target track is picked up operate, after having executed pick-up operation, an iteration Process terminates, and judges whether the also promising road searching according to T flag, if there are the road not searched, then continues to change next time For process, otherwise algorithm terminates.

Further, in pick-up operation, a layer plane label is carried out to sample point, be easy to the differentiation of succeeding layer plane, right As follows in the layer plane label set-up mode of sample point：

First, it is that the currently affiliated layer plane of all of sample point demarcation is numbered 0, then from primary earthquake road, Process first seismic channel adjacent with it, after obtaining the optimum pickup mode between them, split from first sample point Begin Distribution Layer plane label one by one, and layer plane label from the beginning of 1, that is, gives the sample point in optimum pickup mode corresponding to distributing Layer plane label, now just contains two sample points in each layer of plane, and for the not sample in optimum pickup mode Point then still retaining layer plane be numbered 0；Spread in operation and pick-up operation in follow-up, in each optimum pickup mode Each sample point pair, if the sample point belonging on the road of source has had a layer plane label, then corresponding with it belong to target Sample point on road just puts on the layer plane label the same with it, if the layer plane label of two sample points is all 0, according to work as The front layer plane number having existed, increases a layer plane label newly, and described two layers of plane label is all that 0 sample point just becomes For two sample points of new layer plane, the layer plane label updating of two sample points is new layer plane label simultaneously.

Step 4：The preliminary layer plane of pickup is plugged a hole so that the layer plane being obtained is complete.

Last in the present invention, need the layer plane of pickup is plugged a hole and smoothing processing, smoothing processing is taken Mean value method, is his 8- neighbour and he itself totally 9 to the depth of each sample point on the layer plane of pickup by its value The mean value of individual point.If a certain seismic channel do not have picked, but his adjacent road major part is all picked, then just Plugged a hole, the value of this seismic channel is disposed adjacent to the mean value of seismic channel sample point.Through smoothing processing and the place that plugs a hole Reason, obtains last layer plane.

The invention provides a kind of method of new multiple layer of plane of automatic Picking of increase constraint information, the method combination The information of given containment surfaces, there is provided the auxiliary information to the layer plane that will pick up.The sample obtaining after after filtering Carry out spreading method based on seed seismic channel on point and pick up a layer plane, accuracy and time efficiency can be improved.In pick process In, every time two adjacent seismic channels are carried out with the screening of optimum pickup mode, in a pickup mode, contains multipair sample point Right, they belong to different layer planes, when obtaining an optimum pickup mode, that is, these sample points are distributed simultaneously In the layer plane being located to paired sample point, that is, complete the purpose simultaneously picked up.

Claims (6)

1. a kind of method that in 3-D seismic interpretation, the batch based on constraint picks up layer plane, is characterized in that identifying that upper and lower two are put down Sliding layer plane is as containment surfaces, using the known conditions of containment surfaces, for the relation that there is mutually constraint between layer plane, same All layers of plane between step pickup containment surfaces, comprise the following steps：

Step 1：Read three-dimensional data volume data, obtain extreme point with interpolation method, and the extreme point obtaining is filtered process, Obtain sample point；

Step 2：According to the layering of geological type, identify that upper and lower two smooth layer plane, as containment surfaces, obtain in containment surfaces The information of each sample point, as the constraint information of pickup layer plane, described constraint information includes each sample point in containment surfaces Obliquity information and amplitude information；

Step 3：With two containment surfaces of step 2 as boundary, seed seismic channel is set, is scanned for based on " seed point spreads method ", Pickup sample point is progressively extended according to pickup standard, obtains layer plane to be picked up, described pickup standard refers between sample point Similarity set threshold range in：

Primary earthquake road, referred to as seed seismic channel are set first, respectively neighbor search are carried out to its adjacent seismic channel, with seed-breeding field Based on sample point on shake road, judge whether to pick up the sample point of adjacent seismic channel, the binding characteristic condition of search includes picking up Take type, amplitude similarity, depth similarity and inclination angle constraint, for two current sample points, pick up type and refer to two Sample point is maximum or minimizing uniformity judges, is unanimously then picked up；Amplitude similarity refers to two sample points Whether amplitude similarity in the threshold range setting, then picks up；Depth similarity refers to the depth value difference of two sample points Whether in the threshold range setting, then pick up；Inclination angle constraint refers to adjacent up and down two in sample point is for 3D data volume During individual sample point, whether the inclination angle similarity between two sample points in the threshold range setting, then picks up；

Qualified sample point is picked up on the adjacent seismic channel of seed seismic channel according to above-mentioned binding characteristic condition, picked Sample point be located seismic channel become new seed seismic channel, neighbor search is proceeded to new seed seismic channel, if All seismic channels in the neighborhood of one seed seismic channel are all searched, then be labeled as searching by this seed seismic channel Rope, the search after being no longer participate in, when all of seismic channel is all searched once, iteration ends；Pickup in sample point Cheng Zhong, sample point divides automatically according to binding characteristic condition, constitutes different layer planes；

Described include spreading operation and pick-up operation based on the search of " seed point spreads method ", spread operation first from current The seed seismic channel processing is traveled through to the adjacent seismic channel not being searched；Pick-up operation be traversal current some During adjacent seismic channel, according to the pickup standard setting, sample point pick-up operation is carried out to described adjacent seismic channel, specially：

1) spread operation

Spread the operation traversal BFS method based on figure for the flow process, construct a figure corresponding with 3D data volume：Three-dimensional Each of data volume seismic channel is configured to one of figure point, is any two adjacent seismic channel Trace_iAnd Trace_jConstruction A line E_ij, from seed seismic channel corresponding summit, S proceeds by breadth first traversal, until having traveled through be adjacent all Summit, that is, traveled through all seismic channels being adjacent；

Wherein, three Mark Arrays BF, BNF and T are arranged to seismic channel, T records the search condition of seismic channel, if current seismic road Being marked as 1, represent that this seismic channel has been searched, if being marked as 0, representing that this seismic channel is not yet searched；BF records Current seed seismic channel, referred to as boundary trace, if seismic channel carries out neighbor search as seed seismic channel is designated as 1, otherwise for 0；The new seismic channel that BNF record searches from boundary trace, BNF is marked as this seismic channel of 1 expression will be used as next time The boundary trace of search；Spread operation to store the value of BNF in BF first, represent and the seismic channel becoming new boundary trace is set to Current border road, carries out spreading search by current border road, i.e. BFS, and all values of BNF are emptied as 0, In follow-up pick-up operation, recording again new will become the seismic channel of boundary trace；

2) pick-up operation

Current border road is referred to as source road, says from source and the seismic channel traverse is sent to its adjacent direction is referred to as target track, Pick-up operation finds, in target track, the sample point meeting pickup standard, has processed a target track just by this target track every time Picked sample point and corresponding source road sample point add same layer plane, and described correspondence refers to meet the constraint feature bar Part, pick-up operation completes the automatic division of sample point according to binding characteristic condition, obtains each layer of plane；

When pick-up operation is carried out, from Yuan Dao, if the T of the target track traversing is labeled as 1, continue with next neighbour Connect seismic channel；If the T of target track is labeled as 0, then target track is picked up operate, pickup meet the constraint characteristic condition Sample point, this point is added in corresponding layer plane；Setting threshold value t, when the sample having been added to layer plane in a target track The percentage of the number of this point exceedes threshold value t, then arranging the BNF mark of this target track and T flag is 1, has represented this target track Processed and became the boundary trace next time searched for through searched；After road corresponding all target track search in source terminate, complete Pick-up operation, now, recorded new boundary trace in BNF array；

Step 4：The step 3 layer plane that obtain of pickup is smoothed and plugs a hole process so as to get layer plane complete, Complete the batch pickup of all layers of plane between containment surfaces.

2. the method that in a kind of 3-D seismic interpretation according to claim 1, the batch based on constraint picks up layer plane, its Feature is that step 1 is specially：For one piece of given 3D data volume, quadratic polynomial interpolation method is adopted to each seismic channel, Obtain an extreme point according to three serial number points on seismic channel every time, i.e. crest or trough, the extreme point obtaining is adopted Digital averaging filtering method, when occur two or more amplitude be continuously on the occasion of or be continuously negative value extreme point when, by these Continuous extreme point synthesizes a synthetic point, and the coordinate of synthetic point is the coordinate mean value of the extreme point being synthesized, amplitude The mean value of amplitude of the extreme point for being synthesized, after all extreme points are completed with filtering, obtains the sample point of each seismic channel.

3. the method that in a kind of 3-D seismic interpretation according to claim 1, the batch based on constraint picks up layer plane, its Feature is in pick-up operation, and design energy function is used for describing pickup standard, specially：

Choose seismic channel as source road first between containment surfaces, from Yuan Dao, judge about target track whether by Search, if source road is Trace₁, the not searched target track of it some is Trace₂, n₁And n₂Represent this two ground respectively The sample point quantity containing on shake road, is that the sample point in source road and target track puts on sequence number, respectively S_1l,S_2g, wherein l= 1,2,…,n₁, g=1,2 ..., n₂If m represents pickup mode, pickup mode refers to the sample point picking up to set, m={ m₁, m₂,...,m_k, m_pRepresent a sample point pair, m_p=(S_1l,S_2g),l∈[1,n₁],g∈[1,n₂], p ∈ [1, k],

Setting pickup constraints：(1) layer plane can not intersect；(2) phase polarity of layer plane is consistent, design energy function F describes pickup mode, and energy function f illustrates the binding character between the correlation between sample point and layer plane：

F (m)=exp (- U (m)) (1)

U (m) is a linear combination of the restrictive function between relevance function and layer plane between sample point：

U (m)=ω_cC(m)+ω_rR(m) (2)

Function C (m) represent sample point in source road and target track between total amplitude coefficient correlation, function R (m) represents on the road of source Total inclination angle coefficient correlation of all adjacent sample points, ω_cWith ω_rThe weights of C (m) and R (m) respectively, described weights pass through experiment Debugging obtains；

Described amplitude coefficient correlation is：

For a sample point S (x, y, z), each K sample up and down is selected centered on S (x, y, z) on its place seismic channel Point, to be described as the wave character of S (x, y, z) using this 2K+1 sample point { S (x, y, z') | | z'-z |≤K }, i.e. this 2K+1 Individual sample point only has depth value z different, and x and y coordinates are identical；Normalized crosscorrelation technology is selected to describe sample point partial waveform Similitude, function E (x, y, z) represent sample point S (x, y, z) amplitude, for two sample point s₁(x₁,y₁,z₁), s₂ (x₂,y₂,z₂), the normalized-cross-correlation function between them is：

$C(s_1,s_2)=\frac{\underset{\mathrm{\Δ}z=-K}{\overset{+K}{\Σ}}\left(E\right(x_1,y_1,z_1+\mathrm{\Δ}z)*E(x_2,y_2,z_2+\mathrm{\Δ}z\left)\right)}{\sqrt{\underset{\mathrm{\Δ}z=-K}{\overset{+K}{\Σ}}{E}^2(x_1,y_1,z_1+\mathrm{\Δ}z)}*\sqrt{\underset{\mathrm{\Δ}z=-K}{\overset{+K}{\Σ}}{E}^2(x_2,y_2,z_2+\mathrm{\Δ}z)}}---\left(3\right)$

For a pickup mode m={ m₁,m₂,...,m_k, wherein there is k to sample point pair, by the phase relation of each sample point pair Number is cumulative, obtains total sample point amplitude coefficient correlation in pickup mode m in current source road and target track, that is,

$C\left(m\right)=\underset{p=1}{\overset{k}{\Σ}}C\left(m_p\right)---\left(4\right)$

Described inclination angle coefficient correlation is：

With the angle information between two inclination angles as judging both similarity degrees, if dip is (S_1l,S_2g) and dip (S_1l',S_2g') It is respectively the obliquity information of adjacent two sample points on same seismic channel, wherein l=1,2 ..., n₁, g=1,2 ..., n₂, l ' =1,2 ..., n₁, g '=1,2 ..., n₂, using the cosine value of the angle between this two inclination angles as inclination angle similarity, on the road of source The inclination angle similarity factor of two adjacent sample points is defined as following form：

$R(S_{1,l},S_{1,l^{\′}})=\frac{dip(S_{1,l},S_{2,g})*dip(S_{1,l^{\′}},S_{2,g^{\′}})+1}{\sqrt{1+{\mathrm{dip}}^2(S_{1,l},S_{2,g})}*\sqrt{1+{\mathrm{dip}}^2(S_{1,l^{\′}},S_{2,g^{\′}})}}---\left(5\right)$

Mode m={ m thus can be picked up₁,m₂,...,m_kTotal inclination angle coefficient correlation：

$R\left(m\right)=\underset{p=1}{\overset{k-1}{\mathrm{\Σ}}}R(S_{1,p},S_{1,p+1})+R(S_{1,u},S_{1,1})+R(S_{1,k},S_{1,d})---\left(6\right)$

Wherein R (S_1,u,S_1,1) inclining between containment surfaces are intersected with source road in expression sample point and first sample point in source road The similarity factor at angle, R (S_1,k,S_1,d) represent sample point that on the road of source, last sample point and lower containment surfaces are intersected with source road it Between inclination angle similarity factor；

After the solution to amplitude coefficient correlation and inclination angle coefficient correlation, obtain different pickup mode institutes between seismic channel right The energy function answered, finds an optimum pickup modeMake energy functionFor minimum, that is,To meet following shape Formula：

$\tilde{m}=\arg \underset{m\∈M}{min}f\left(m\right)---\left(7\right)$

Wherein, M represents the set of all possible pickup mode, for optimum pickup modeIncluded in all sample points Right, each sample point centering, the sample point of target is assigned in the layer plane belonging to source road sample point, thus picks up simultaneously Sample point in multiple layers of plane, realizes batch and picks up；

In conjunction with formula (1) to formula (7), target track is picked up operate, after having executed pick-up operation, an iteration process Terminate, judge whether the road also not searched according to T flag, if there are the road not searched, then continue next iteration mistake Journey, otherwise algorithm terminate.

4. the method that in a kind of 3-D seismic interpretation according to claim 3, the batch based on constraint picks up layer plane, its Feature is in pick-up operation, carries out a layer plane label to sample point, is easy to the differentiation of layer plane, for the layer plane of sample point Label set-up mode is as follows：

First, it is that the currently affiliated layer plane of all of sample point demarcation is numbered 0, then from primary earthquake road, process First seismic channel adjacent with it, obtain they between optimum pickup mode after, from first sample point to beginning by Individual Distribution Layer plane label, layer plane label from the beginning of 1, that is, gives the sample point in optimum pickup mode to the corresponding layer position of distribution Face label, now just contains two sample points in each layer of plane, and for the not sample point in optimum pickup mode then Still retaining layer plane is numbered 0；Spread in operation and pick-up operation in follow-up, for every in each optimum pickup mode One sample point pair, if the sample point belonging on the road of source has had a layer plane label, then corresponding with it belong in target track Sample point just put on the layer plane label the same with it, if the layer plane label of two sample points is all 0, according to current Layer plane number through existing, increases a layer plane label newly, and described two layers of plane label is all that 0 sample point just becomes new Two sample points of layer plane, the layer plane label updating of two sample points is new layer plane label simultaneously.

5. the method that in a kind of 3-D seismic interpretation according to claim 1, the batch based on constraint picks up layer plane, its Feature is in step 4, and smoothing processing takes mean value method, the depth to each sample point on the layer plane of pickup, will be described Depth value is the mean value of the 8- neighbour of sample point and sample point itself totally 9 points；If a certain seismic channel is not picked up Get, but its adjacent seismic channel has more than 6 to be all searched, then carry out process of plugging a hole, this is not picked To the value of seismic channel be disposed adjacent to the mean value of seismic channel sample point, through smoothing processing and process of plugging a hole, obtain Layer plane afterwards.

6. the method that in a kind of 3-D seismic interpretation according to claim 1, the batch based on constraint picks up layer plane, its Feature is in step 2, and the layer bit-plane positions of pickup according to required for explanation personnel, by artificial pickup or other existing modes Corresponding pickup containment surfaces.