CN106023197B - A kind of method of upright rock stratum automatic identification and extraction - Google Patents
A kind of method of upright rock stratum automatic identification and extraction Download PDFInfo
- Publication number
- CN106023197B CN106023197B CN201610333329.9A CN201610333329A CN106023197B CN 106023197 B CN106023197 B CN 106023197B CN 201610333329 A CN201610333329 A CN 201610333329A CN 106023197 B CN106023197 B CN 106023197B
- Authority
- CN
- China
- Prior art keywords
- line
- point
- straight
- flp
- rock stratum
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20048—Transform domain processing
- G06T2207/20061—Hough transform
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30184—Infrastructure
Landscapes
- Image Analysis (AREA)
- Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)
Abstract
The invention discloses one kind based on geology dignity figure layer and contour figure layer, the method that automatic identification extracts upright rock stratum.This method comprises the following steps：1) the rock stratum face element in geologic map is pre-processed, eliminates the influence on polygon lagoon island (hole) and map sheet boundary, obtain linear bedding margin element；2) Hough transformation is carried out to the point in bedding margin element, bedding margin element is transformed to one group of curve.By the distribution of judgment curves intersection point whether in preset window, to realize the extraction of bedding margin element straight portion；3) filtering that fault line filtering, gentle boundary line filtering and non-parallel pairs of boundary line are carried out to the straight bedding margin extracted, filters out the boundary line of the pairs of upright rock stratum of composition, and draws upright rock stratum.
Description
Technical field
The invention belongs to geographical information technology application fields, and in particular to one kind being based on geology dignity figure layer and contour map
Layer, with Hough transformation algorithm, the method for realizing upright rock stratum automatic identification and extraction.
Background technology
According to construction occurrence type, rock stratum can be divided into horizontal stratum, tilted stratum and upright rock stratum.Wherein, upright rock stratum
Refer to the rock stratum that inclination angle is more than 85 °, is generally present in and constructs strong area.Its geological boundary does linear extension along its trend,
Not by the influence of topography.In conventional GEOLOGICAL APPLICATION, mainly manually finds boundary by reading geologic map and show as straight parallel
The element of line morphology identifies upright rock stratum.And existing computer automation interpreting blueprints technology is mainly used in architectural drawing field
(Yang Huafei, Yang Ruoyu, Lu Tong wait manual reading of drawings Analysis on Mechanism and its application [J] the computers section in developing computer reading graph
Learn, 2008,35 (2)), it is less for the research based on map knowledge automation interpreting blueprints.The interpreting blueprints of geologic map is still by special
Industry personnel combine knowledge rule, the mode for extracting geologic elements on demand to realize.But manual identified method, efficiency is low, accuracy because
People and it is different, quality is difficult to ensure.
Invention content
Technical problem to be solved by the present invention lies in overcome defect of the existing technology, it is proposed that a kind of upright rock
The method of layer automatic identification and extraction, using upright rock stratum, exposure is this feature of straight border in geologic map, based on suddenly
The Straight Line Identification algorithm of husband's transformation extracts bedding margin straight in geologic map, and is filtered by tomography filtering, gentle boundary line
It is filtered with non-parallel pairs of boundary line to screen the boundary line for forming upright rock stratum, finally realize the automatic identification of upright rock stratum and carried
It takes, realizes computer automatic drafting of the upright rock stratum in all kinds of geologic maps.
Differentiate that the Main Basiss of upright rock stratum are the presence of straight geological boundary.Hough transformation in image procossing as identifying
One of basic skills of geometry, initial design are exactly to be extracted for realizing the identification of graph line.For this purpose, this patent intends base
In Hough transformation method, by identification by stages to straight portion in geological boundary and extraction, the automatic knowledge of upright rock stratum is realized
Not.
The main thought of this patent is：1) the geology dignity element processing in geologic map is wanted for linear bedding margin
Element, and eliminate the influence on polygon lagoon island (hole) and map sheet boundary；2) carrying out Hough transformation to the point in bedding margin element can
Bedding margin element is transformed to one group of curve, by analyzing this group of curve, can determine whether current line feature is to meet the requirements
Straight line；3) fault line filtering, the filtering of gentle boundary line and non-parallel pairs of boundary line are carried out to the straight bedding margin extracted
Filtering, filter out the boundary line for forming upright rock stratum, and draw upright rock stratum.
The method of the present invention upright rock stratum automatic identification and extraction, is as follows：
Step 1: formation boundaries are read and pretreatment
Step 1.1 loads the rock stratum face element figure layer StratumLayer of shp formats, obtains all face elements combinations
Stratums=stai | and i=1,2,3 ..., n }, n is the quantity of rock stratum face element, staiIt is each rock stratum face element.
Step 1.2 reads each rock stratum face element sta in set StratumsiBoundary, save as line feature set OriL
={ oli| i=1,2,3 ..., n }, each face element staiBoundary be stored as a line feature oli；Each line feature oli
A corresponding point set OLPi={ opij(xij, yij) | j=1,2,3 ..., mi, point opijFor point set OLPiIn j-th point, j is
Point opijIn point set OLPiIn subscript, coordinate representation be (xij, yij)。miFor line feature oliCorresponding point set OLPiIn point
Quantity；oli" id " attribute is stored to record the rock stratum face element sta of its ownershipiNumber.
There may be part line feature ol in step 1.3 line feature set OriLiFor multi-part element, to multi-part element
oliIt is split, obtains new line feature set SinL={ slk| k=1,2,3 ..., l }, l is rock stratum face element side after splitting
The quantity summation on boundary, line feature slkAttribute " id ", which is inherited, splits preceding oliThe attribute " id " of line feature.
Line feature sl is certainly existed in step 1.4 line feature set SinLk, part is and non-genuine positioned at map sheet boundary
Formation boundaries, need to delete the part positioned at formation boundaries.
Step 1.5：Traverse the line feature set SinL={ sl after formation boundaries are deletedk| k=1,2,3 ..., l }, if depositing
It is identical in attribute " id ", and the line feature sl of Spatial AdjacencyiWith slj(i, j=1,2,3 ..., l, and j ≠ i), to line feature sli
With sljIt merges, obtains new line feature set FinL={ flu| u=1,2,3 ..., p }, p is the rock after the completion of pretreatment
The quantity summation of stratum boundary line.
Step 2: the straight borderline distilling in rock stratum
Straight line set AllSL is established for storing the straight boundary line that all line features extract, to line feature set FinL
In each line feature fluCarry out straight borderline distilling.For single line feature flu, create an empty straight line set
StaightL stores fluIn the straight line that extracts, line feature fluStraight borderline distilling method it is as follows：
Step 2.1 gives a length standard LlimitCome limit extraction straight line shortest length.
The each line feature fl of step 2.2uExpression is point set FLPu={ flpuv| v=1,2,3 ..., qu, quFor line flu
On number of vertices；For each point flpuv(xuv, yuv) carry out Hough transformation.The principle of Hough transformation is cartesian coordinate system
On line be mapped as a bit by Hough transformation, point by Hough transformation is mapped as a curve.Set angle tolerance Δ θ is used as and carries
Limitation is moved towards in the straight boundary line taken, i.e. Δ θ is bigger, and the straight boundary line of extraction may be more bent, conversely, the straight boundary line of extraction
Then closer to straight line.Pure straight line is less in geologic map, and therefore, it is necessary to adjust Δ θ to extract most suitable straight boundary
Line.Hough transformation is carried out to a point to take using angular tolerance Δ θ as interval Corresponding ρ can be obtained according to formula (1)0, ρ1, ρ2, ρ3,…,ρj, point
(θ0,ρ0),(θ1,ρ1),(θ2,ρ2),(θ3,ρ3),…,(θd,ρd) curve the HoughL { (θ after a Hough transformation can be linked to be0,
ρ0),(θ1,ρ1),(θ2,ρ2),…,(θd,ρd)}；
ρ=x cos θ+y sin θs (1)
Then for point flpuv, using Δ θ as interval, by its coordinate (xuv, yuv) bring formula (1) into, point flp can be obtaineduvHough
Curve HoughL after transformationuv={ (0, xuv), (Δ θ, xuvcosΔθ+yuvSin Δ θ), (Δ θ * 2, xuvcos(Δθ*2)+
yuvSin (Δ θ * 2)), (π, xuv)}。
Step 2.3 creates a straight line point set StaightLineP, the point in alignment for storage group；It records first
It is 1 to start subscript s, by flpusSet StaightLineP is added.
Step 2.4 is by FLPuIn point be added to StaightLineP successively, every time into StaightLineP addition newly
When point, i.e., StaightLineP is proceeded as follows：
1) if the points that StaightLineP includes are less than 3, continue to add new point into StaightLineP；Conversely,
Seek all the points flp in StaightLinePutCorresponding HoughLut(t=1,2 ..., c, c are in current StaightLineP
The quantity of point)；
2) ρ is takenmaxFor all HoughLutIn maximum ρ values, take ρminFor minimum ρ values by tolerance Δ θ and Δ ρ=
(ρmax-ρmin)/d obtains margin window (Δ θ, Δ ρ)；
3) HoughL is calculated successivelyutWith HoughLu(t+1)Intersection point pt(θt(t+1), ρt(t+1)) (t=1,2,3 ..., c-1),
Obtain point set InterP={ p1(θ12, ρ12),p2(θ23, ρ23),p3(θ34, ρ34),…,pc-1(θ(c-1)c, ρ(c-1)c), each intersection point
Pt represents the straight line in cartesian coordinate system, coordinate (θt(t+1), ρt(t+1)) be the straight line Hough transformation after coordinate.Sentence
Whether being distributed for all intersection points of breaking is more than the margin window (Δ θ, Δ ρ) set, and only there are one intersection points then without the judgement：
If the distribution of all intersection points is without departing from tolerance, and the point flp newly addedutIt is not line fluUpper last point then continues to add
Add some points, repeat the 1 of step 2.4) -3) step；
If the distribution of all intersection points is without departing from the tolerance of setting, and the point flp newly addedutAs line fluOn last point,
Then by StaightLineP={ flpus, flpu(s+1)..., flputIt is saved as straight boundary line li, its place is calculated according to formula (2)
Hough coordinate (the l θ of straight linex, l ρx), it stores to straight boundary line liAttribute " θ " and " ρ "；By straight boundary line liSet is added
StaightL completes line feature fluThe extraction of cathetus element；
Wherein, l θxWith l ρxFor straight boundary line liThe point coordinates of Hough transformation mapping, θt(t+1), ρt(t+1)For intersection point ptSeat
Mark, c are the quantity at the current midpoints StaightLineP.
If the distribution of all intersection points exceeds the tolerance of setting, and will point flputAfter being removed in StaightLineP, line SL<
flpus,flpu(t-1)>Length greatly be equal to Llimit, then by StaightLineP={ flpus, flpu(s+1)..., flpu(t-1)Deposit
At straight boundary line li, Hough coordinate (the l θ of straight line where calculating it according to formula (2)x, l ρx), it stores to straight boundary line liCategory
Property " θ " and " ρ ".By straight boundary line liSet StaightL is added.Then StaightLineP is emptied, record starts subscript again
S is t-1, by point flpu(t-1)StaightLineP is added, repeats the 1 of step 2.4) -3) step.
If the distribution of all intersection points exceeds the tolerance of setting, and will point flputAfter being removed in StaightLineP, line SL<
flpus,flpu(t-1)>Length be less than Llimit, then will point flpusIt is removed from StaightLineP, record starts subscript s again
For s+1, by point flpu(t-1)StaightLineP is added, repeats the 1 of step 2.4) -3) step.
Step 2.5 is for the straight boundary line l in StaightLiWith lj(j ≠ i), if meeting：Vertex overlaps and attribute " θ " is poor
Value is less than tolerance Δ θ, then by two straight boundary line liWith ljIt merges.
All straight boundary line l in step 2.6StaightLiAttribute " id " inherit line feature fluAttribute of the same name.And
All straight line set AllSL are added in element in StaightL.
Step 3: upright rock stratum identification is drawn
Step 3.1 extracts straight boundary line to all line features and obtains straight line set AllSL={ li| i=1,2,3 ...,
W }, w is the quantity summation to the straight boundary line of all rock stratum.
Step 3.2 carries out dip fault filtering to straight line set AllSL；For the straight line element point of one group of phase mutual overlapping
Rock stratum " stratum symbol " sequence for not recording this group of straight line element both sides, if the sequence of both sides is consistent, then it is assumed that this organizes straight boundary
Line forms a fault line, and operation is filtered to this straight boundary line of group.
Step 3.3 loads contour figure layer, obtains contour set Contour={ clj| j=1,2,3 ..., r }, base
Gentle boundary line filtering is carried out to AllSL in contour；Angle threshold value α is set, for straight boundary line li, journal its with it is contour
The intersection point of line, and calculate point of intersection contour and boundary line liAngle { lci1, lci2..., lcik}；Seek average angle lciIf
It is less than the angle threshold value α, AllSL=SSL of setting1∪SSL2∪SSL3..., SSLg, g is still to retain straight boundary line after filtering
Rock stratum quantity；There is the rock stratum in pairs of parallel straight boundary line to be regarded as upright rock stratum, non-parallel pairs of straight boundary line is carried out
Filtering；It is the threshold value for screening parallel straight boundary line that the present invention, which tries to please and limits Δ θ, judges rock stratum SSLtIn whether meet, exist in pairs
Straight boundary line lxWith ly(y ≠ x) attribute " θ " difference is less than tolerance Δ θ, and the straight boundary line to being unsatisfactory for above-mentioned condition carried out
Filter operation；
Step 3.5 is for meeting the rock stratum of upright rock stratum condition, straight boundary line lxWith lyThe part to fence up is as upright
, straight boundary line l is connected automaticallyxWith lyHead and the tail point, draw face element sltai, as upright rock stratum；Finally obtain upright rock
Layer set LStatum={ lstai| i=1,2,3 ..., f }, f is the quantity for the upright rock stratum that identification obtains.
The method of the present invention is calculated using geometry the process of manual reading of drawings being changed into meter in conjunction with the knowledge that geologic map is interpreted blueprints
Calculate machine-readable figure.In conjunction with address interpreting blueprints knowledge, the automatic identification of upright rock stratum is realized more with respect to manual identified using geometric ways
Accelerate prompt, accuracy higher.Geologic map is read for ordinary person and provides help, can also meet Geologic modeling, tectonic knot etc.
The demand of the production and construction such as academic research and comparison and choice.
Description of the drawings
The flow chart of Fig. 1 the method for the present invention.
Fig. 2 embodiments experimental data (terrain and geologic map in somewhere).
Fig. 3 embodiment formation boundaries line drawing effects.
Fig. 4 example points Hough transformation line charts.
Fig. 5 example lines detection design sketch.
Fig. 6 embodiment lines detection result figures.
Fig. 7 example tomographies filter schematic diagram.
Fig. 8 embodiment tomography filter effect figures.
Fig. 9 example contours filter schematic diagram.
The gentle boundary line filter effect figure of Figure 10 embodiments.
The non-parallel boundary line filter effect figure in Figure 11 embodiments rock stratum.
Extract result in the upright rock stratum of Figure 12 embodiments.
Specific implementation mode
Below in conjunction with the accompanying drawings and by describing an automatic identification and marking the example of upright rock stratum, to further illustrate this
The effect of invention.It is that sample data (mainly divides including contour figure layer and geology that the present embodiment, which selects the terrain and geologic map in somewhere,
Area's figure layer), as shown in Figure 2.
(1) formation boundaries read and pre-process
Step 1.1 loads the geology dignity figure layer data " geologic division " of shp formats, reads rock stratum face element and reads in set
Stratums={ stai| i=1,2,3 ..., 25 }, include 25 rock stratum face element altogether in the present embodiment.
Step 1.2 load terrain and geologic map figure layer basic parameter be：Wide 2316.33m, high 1705.86m, map sheet boundary
Point set BoundaryP=(461.44, -1999.47), (507.27, -1999.54), (716.37, -1999.88), (1,
646.99,-2001.37),(1749.90,-2001.53),(1,891.30,-2001.76),(2378.58,-2,002.53),
(2459.91,-2002.66),(2777.77,-2003.17),(2769.89,-1095.25),(2769.89,-1095.25),
(2763.88,-402.55),(2762.97,297.38),(1941.99,-298.70),(1686.37,-299.11),
(1302.08,-299.73),(1058.81,-300.12),(817.65,-300.50),(700.11,-300.69),
(457.74,-301.08),(459.83,1259.61),(460.39,-1514.49)}.Embodiment straight line full-length LlimitIt takes
Value 230m；The window 15*15 of Hough transformation, tolerance Δ θ values are π/15, i.e. 0.209rad, while this value is provided parallel to line
The threshold value of screening.The angle threshold value of gentle boundary line filtering is set as 20 °, i.e. 0.349rad.
Step 1.3 reads staiBoundary, save as line feature oli, obtain set OriL={ oli| i=1,2,3 ...,
25}.Each face element staiIt is stored as an oli。oliInherit stai" id " attribute.
The case where island (hole) is not present in step 1.4 the present embodiment, all boundary lines are all single loop wire.Obtain SinL=
{sli| i=1,2,3 ..., 25 }, sliInherit oli" id " attribute；The point set positioned at map sheet boundary is recycled to divide line feature,
And according to condition screened and connected, finally obtain set FinL={ flu| u=1,2,3 ..., 28 }, with ol3For, boundary
Line is obtained boundary line fl by boundary segmentation3And fl4, and fl3And fl4" Id " attribute be 3.
Step 1.5 full figure takes boundary to rock stratum, and it is as shown in Figure 3 to obtain pretreated geological boundary figure layer.
(2) the straight borderline distilling in rock stratum
Lines detection extracts respectively for each line feature, with line feature fl3For, point set FLP3=
(507.27, -1999.54), (509.12, -1992.93), (518.21, -1985.494) ..., (1730.64, -
1981.40), (1749.90, -2001.53) }, including 38 points.
Step 2.1 is to point set FLP3Hough transformation is carried out, is the Hough transformation curve that interval calculation obtains each point with π/15.
With point flp31Curve HoughL is obtained for (507.27, -1999.54), after Hough transformation31<(0.507.27), (0.21,
80.46), (0.42, -349.87), (0.63, -764.91), (0.84, -1146.52) ..., (2.93, -911.91), (3.14,
507.27)>, it is as shown in Figure 4 to draw curve.
Step 2.2 creates a straight line point set StaightLineP, by FLP3In first point flp31Rectilinear point is added
Collection includes at this time only an element in StaightLineP.
Step 2.3 is by FLP3In point be added to StaightLineP successively, every time into StaightLineP addition newly
When point, i.e., StaightLineP is operated.With line feature fl3For, by point flp32Set StaightLineP is added, it is right
It is proceeded as follows：
1) only include at this time two points in StaightLineP, then add new point flp into set33。
2) point flp is obtained31, point flp32With flp33Curve HoughL after Hough transformation31、HoughL32With HoughL33。
Obtain ρmax=518.21, ρminΔ ρ=171.78 are calculated in=- 2058.43, determine margin window be (0.21,
171.78)。
3) HoughL is calculated31With HoughL32Intersection point obtain point p1(2.87, -1021.44) calculate HoughL2With
HoughL3Intersection point obtain point p2(2.27, -1857.81), the difference (0.60,836.37) of two intersection points exceed margin window, then
It will point flp33It is removed from StaightLineP, calculates line SL<flp31,flp32>Length, value 6.86m, be less than Llimit,
Then by flp31It is removed from StaightLineP, by point flp33Rejoin StaightLineP.
4) only include at this time two points in StaightLineP, then add new point flp into set34。
5) point flp is obtained32, point flp33With point flp34Set HoughL after Hough transformation32、HoughL33With
HoughL34.Obtain ρmax=542.38, ρminΔ ρ=173.39 are calculated in=- 2058.43, determine that margin window is
(0.21,171.78).
6) HoughL is calculated32With HoughL33Intersection point obtain point p1(2.27, -1857.81) calculate HoughL33With
HoughL34Intersection point obtain point p2(2.20, -1901.60), the difference (0.07,43.79) of two intersection points are no more than margin window,
Continue to add new point flp35, repeat to judge StaightLineP.
7) addition point fpl36、fpl37、flp38,…,flp311, in set StaightLineP after Hough transformation curve friendship
Point is all in margin window, until adding point flp312(875.36, -1751.91), at this point for Hough curve { HoughL32,
HoughL33,…,HoughL312, margin window is (0.21,293.88), intersection point collection InterP={ p1(2.27 ,-
1857.81), p2(2.20, -1901.60), p3(2.27, -1887.48) ..., p9(2.11, -1944.95), p10(2.04 ,-
1950.72) } intersection point distribution has exceeded margin window, then by flp312It is removed from StaightLineP, calculates line SL<flp32,
flp311>Length, value 385.56m, be more than Llimit, then by StaightLineP={ flp32, flp33..., flp311}
Save as straight line element l5, it is added into set StaightL, by calculating, attribute " θ " is assigned a value of 2.09, and attribute " ρ " is assigned
Value is -1918.73.
8) StaightLineP is emptied, by point flp311Set is added, continues to judge a to the last point flp338。
After step 2.4 completes straight borderline distilling to straight line, adjacent and approximate line in set StaightL is carried out
Merge.With line feature fl3For, it can extract out 2 sections of straight boundary line { l5, l6, wherein l5With l6Attribute " θ " and " ρ " gap are small
It in tolerance, and abuts, then to l5With l6It merges, straight borderline distilling effect is carried out to it and merging effect is as shown in Figure 5.
Step 2.5 carries out the element in set StaightL the assignment of attribute " id ", with line feature fl3For, it is right
It should { l in set StaightL5" id " attribute assignment be 3.
(3) upright rock stratum identification is drawn
After step 3.1 completes straight borderline distilling, the straight boundary line of all rock stratum is obtained as shown in fig. 6, the straight boundary of full figure
Line set AllSL={ li| i=1,2,3 ..., 48 }.
Step 3.2 carries out dip fault filtering, the rock symbol of the straight line element both sides overlapped to one group to AllSL
Sequence is judged.As shown in fig. 7, line l3、l8With line l43、l10It is overlapped, and l3、l8Its stratum symbol sequence of the rock stratum of ownership
It is classified as { " C_2 ", " C_3 " }, l43、l10Its stratum symbol sequence of the rock stratum of ownership is { " C_2 ", " C_3 " }, and two sequences are identical,
Then judge line l3、l8、l43And l10One group of fault line is formed, by line l3、l8、l43And l10It is removed from set AllSL.It meets sometimes
To two unequal situations of sequence length, then judged as standard comprising another one using one.It is obtained after tomography filters
Set AllSL={ li| i=1,2,3 ..., 39 }, figure layer is drawn as shown in Figure 8.
Step 3.3 loads contour figure layer, obtains contour set Contour={ clj| j=1,2,3 ..., 78 }.Base
Gentle boundary line filtering is carried out to AllSL in contour.With straight line element l3、l4、l5、l6、l7And l8For (as shown in Figure 9),
Attribute is as shown in the table with estimate of situation：
Gentle straight line filtering is carried out one by one to the straight line element in AllSL, finally obtains set AllSL={ li| i=1,
2,3 ..., 25 }, figure layer is drawn as shown in Figure 10.
Step 3.4 is by the equal straight line l of attribute " Id "xIdentity set SSL is addedt, finally obtain set AllSL=SSL1
∪SSL2∪SSL3..., SSL15, it is 2,4,6,8,10,11,12,13,14,15,16,19,20,24 to correspond to attribute " id " respectively
Rock stratum face element with 25 is through passing fault filtering and the filtered straight boundary line in gentle boundary line.Non-parallel pairs of boundary is carried out to it
Line is filtered, and judges SSL successivelytIn whether there is pairs of straight line element lxWith ly, attribute " θ " difference be less than tolerance Δ θ,
Straight boundary to being unsatisfactory for condition is filtered, and the boundary line for finally obtaining the condition of satisfaction see the table below (as shown in figure 11) information：
Step 3.5 is for meeting the element of upright rock stratum condition, by its pairs of straight line element lxWith lyHead and the tail connection is drawn
New face element, the as upright part in rock stratum.5 upright rock stratum parts are separately connected, it is as shown in figure 12 to obtain upright rock stratum.
Claims (1)
1. a kind of method of upright rock stratum automatic identification and extraction, is as follows：
Step 1: formation boundaries are read and pretreatment
Step 1.1 loads the rock stratum face element figure layer StratumLayer of shp formats, obtains all face elements combination Stratums
={ stai| i=1,2,3 ..., n }, n is the quantity of rock stratum face element, staiIt is each rock stratum face element；
Step 1.2 reads each rock stratum face element sta in set StratumsiBoundary, save as line feature set OriL={ oli
| i=1,2,3 ..., n }, each face element staiBoundary be stored as a line feature oli；Each line feature oliCorresponding one
A point set OLPi={ opij(xij, yij) | j=1,2,3 ..., mi, point opijFor point set OLPiIn j-th point, j be point opij
In point set OLPiIn subscript, coordinate representation be (xij, yij), miFor line feature oliCorresponding point set OLPiIn point number
Amount；oliId attributes are stored to record the rock stratum face element sta of its ownershipiNumber；
Step 1.3 is for multi-part line feature ol in part present in line feature set OriLiIt is split, obtains new line and want
Element set SinL={ slk| k=1,2,3 ..., l }, l is the quantity summation on rock stratum face element boundary after splitting, line feature slk
Attribute id, which is inherited, splits preceding oliThe attribute id of line feature；
Step 1.4 for present in line feature set SinL partly positioned at map sheet boundary and the line of not true formation boundaries
Element slk, the part of formation boundaries is located to it and is deleted；
Step 1.5：Traverse the line feature set SinL={ sl after formation boundaries are deletedk| k=1,2,3 ..., l }, belong to if existing
Property id is identical, and the line feature sl of Spatial AdjacencyiWith slj, i, j=1,2,3 ..., l, and j ≠ i, to line feature sliWith sljInto
Row merges, and obtains new line feature set FinL={ flu| u=1,2,3 ..., p }, p is the bedding margin after the completion of pretreatment
Quantity summation；
Step 2: the straight borderline distilling in rock stratum
Straight line set AllSL is established for storing the straight boundary line that all line features extract；To in line feature set FinL
Each line feature fluCarry out straight borderline distilling；For single line feature flu, create an empty straight line set StaightL and deposit
Put fluIn the straight line that extracts, line feature fluStraight borderline distilling scheme it is as follows：
Step 2.1 gives a length standard LlimitCome limit extraction straight line shortest length；
The each line feature fl of step 2.2uExpression is point set FLPu={ flpuv| v=1,2,3 ..., qu, quFor line fluOn
Number of vertices；The line that cartesian coordinate is fastened is mapped as a bit by Hough transformation, and point is mapped as a curve by Hough transformation,
To line fluOn each of point flpuv(xuv, yuv) carry out Hough transformation；For point flpuv, using Δ θ as interval, take θ0=0, θ1=
Δ θ, θ2=Δ θ * 2, θ3=Δ θ * 3 ..., θd=π (d=[π/Δ θ]), ρ is calculated separately according to formula (1)0, ρ1, ρ2, ρ3...,
ρd, obtain point flpuvCurve point set HoughL after Hough transformationuv={ (0, xuv), (Δ θ, xuvcosΔθ+yuvSin Δ θ),
(Δ θ * 2, xuvcos(Δθ*2)+yuvSin (Δ θ * 2)), (π ,-xuv)}；
ρ=xcos θ+ysin θ (1)
Wherein：Δ θ is the angular tolerance of Hough transformation, and d is the interval quantity that angle π is divided under this tolerance；
Step 2.3 creates a straight line point set StaightLineP, the point in alignment for storage group；Record starts first
Subscript s is 1, by flpusSet StaightLineP is added；
Step 2.4 is by FLPuIn point be added to StaightLineP successively, every time into StaightLineP add newly put when,
StaightLineP is proceeded as follows：
1) if the points that StaightLineP includes are less than 3, continue to add new point into StaightLineP；Conversely, seeking
All the points flp in StaightLinePutCorresponding HoughLut, t=1,2 ..., c, c is the current midpoints StaightLineP
Quantity；
2) ρ is takenmaxFor all HoughLutIn maximum ρ values, take ρminFor minimum ρ values, by tolerance Δ θ and Δ ρ=(ρmax-
ρmin)/d obtains margin window (Δ θ, Δ ρ)；
3) HoughL is calculated successivelyutWith HoughLu(t+1)Intersection point pt(θt(t+1), ρt(t+1)) (t=1,2,3 ..., c-1), it obtains
Point set InterP={ p1(θ12, ρ12), p2(θ23, ρ23), p3(θ34, ρ34) ..., pc-1(θ(c-1)c, ρ(c-1)c), each intersection point ptGeneration
Straight line in table cartesian coordinate system, coordinate (θt(t+1), ρt(t+1)) be the straight line Hough transformation after coordinate；Judge institute
Have whether the distribution of intersection point is more than the margin window (Δ θ, Δ ρ) set, only there are one intersection points then without the judgement：
If the distribution of all intersection points is without departing from tolerance, and the point flp newly addedutIt is not line fluUpper last point then continues to add
Point repeats the 1 of step 2.4) -3) step；
If the distribution of all intersection points is without departing from the tolerance of setting, and the point flp newly addedutAs line fluOn last point, then will
StaightLineP={ flpus, flpu(s+1)..., flputIt is saved as straight boundary line li, straight line where calculating it according to formula (2)
Hough coordinate (l θx, l ρx), it stores to straight boundary line liAttribute θ and ρ；By straight boundary line liSet StaightL is added, it is complete
At line feature fluThe extraction of cathetus element；
Wherein, l θxWith l ρxFor straight boundary line liThe point coordinates of Hough transformation mapping, θt(t+1), ρt(t+1)For intersection point ptCoordinate, c is
The quantity at the current midpoints StaightLineP；
If the distribution of all intersection points exceeds the tolerance of setting, and will point flputAfter being removed in StaightLineP, line SL<
flpus, flpu(t-1)>Length be more than or equal to Llimit, then by StaightLineP={ flpus, flpu(s+1)...,
flpu(t-1) it is saved as straight boundary line li, Hough coordinate (the l θ of straight line where calculating it according to formula (2)x, l ρx), it stores to straight
Boundary line liAttribute " θ " and " ρ "；By straight boundary line liSet StaightL is added；Then StaightLineP is emptied, is remembered again
It is t-1 that record, which starts subscript s, by point flpu(t-1)StaightLineP is added, repeats the 1 of step 2.4) -3) step；
If the distribution of all intersection points exceeds the tolerance of setting, and will point flputAfter being removed in StaightLineP, line SL<
flpus, flpu(t-1)>Length be less than Llimit, then will point flpusIt is removed from StaightLineP, record starts subscript s again
For s+1, by point flpu(t-1)StaightLineP is added, repeats the 1 of step 2.4) -3) step；
Step 2.5 is for the straight boundary line l in StaightLiWith lj(j ≠ i), if meeting：Vertex overlaps and attribute θ differences are less than
Tolerance Δ θ, then by two straight boundary line liWith ljIt merges；
All straight boundary line l in step 2.6StaightLiAttribute id inherit line feature fluAttribute of the same name, and will
All straight line set AllSL are added in element in StaightL；
Step 3: upright rock stratum identification is drawn
Step 3.1 extracts straight boundary line to all line features and obtains straight line set AllSL={ li| i=1,2,3 ..., w }, w is
To the quantity summation in the straight boundary line of all rock stratum；
Step 3.2 carries out dip fault filtering to straight line set AllSL；The straight line element of one group of phase mutual overlapping is remembered respectively
Rock stratum " stratum symbol " sequence for recording this group of straight line element both sides, if the sequence of both sides is consistent, then it is assumed that this organizes straight boundary line group
At a fault line, operation is filtered to this straight boundary line of group；
Step 3.3 loads contour figure layer, obtains contour set Contour={ clj| j=123 ..., r }, it is based on contour
Gentle boundary line filtering is carried out to AllSL；Angle threshold value α is set, for straight boundary line li, the friendship of journal itself and contour
Point, and calculate point of intersection contour and boundary line liAngle { lci1, lci2..., lcik}；Seek average angle lciIf it is less than
The angle threshold value α of setting, then by straight boundary line liIt is considered as gentle geological boundary and is filtered operation；
Step 3.4 is by straight boundary line l equal attribute idxIdentity set SSL is addedt, finally obtain set AllSL=SSL1∪
SSL2∪SSL3..., SSLg, g is the rock stratum quantity for still retaining straight boundary line after filtering；There is the rock in pairs of parallel straight boundary line
Layer is regarded as upright rock stratum, is filtered to non-parallel pairs of straight boundary line；The present invention tries to please limit Δ θ to screen parallel straight boundary
The threshold value of line judges rock stratum SSLtIn whether meet, there are pairs of straight boundary line lxWith ly, y ≠ x, attribute θ differences, which are less than, to be held
Δ θ is limited, the straight boundary line to being unsatisfactory for above-mentioned condition is filtered operation；
Step 3.5 is for meeting the rock stratum of upright rock stratum condition, straight boundary line lxWith lyThe part to fence up be it is upright,
Automatically straight boundary line l is connectedxWith lyHead and the tail point, draw face element sltai, as upright rock stratum；Finally obtain upright rock stratum collection
Close LStatum={ lstai| i=1,2,3 ..., f }, f is the quantity for the upright rock stratum that identification obtains.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610333329.9A CN106023197B (en) | 2016-05-18 | 2016-05-18 | A kind of method of upright rock stratum automatic identification and extraction |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610333329.9A CN106023197B (en) | 2016-05-18 | 2016-05-18 | A kind of method of upright rock stratum automatic identification and extraction |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106023197A CN106023197A (en) | 2016-10-12 |
| CN106023197B true CN106023197B (en) | 2018-11-09 |
Family
ID=57097654
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610333329.9A Active CN106023197B (en) | 2016-05-18 | 2016-05-18 | A kind of method of upright rock stratum automatic identification and extraction |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106023197B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106651936B (en) * | 2016-10-19 | 2019-06-18 | 南京师范大学 | A kind of adaptive determination method of the earth's surface exposure attitude of rocks |
| CN106777117B (en) * | 2016-12-15 | 2020-04-03 | 南京师范大学 | Automatic identification method for horizontal rock stratum structure landform |
| CN106709439B (en) * | 2016-12-16 | 2020-04-03 | 南京师范大学 | Automatic identification method for monoclinic rock stratum structure landform |
| CN106600661B (en) * | 2016-12-20 | 2019-06-07 | 黄河勘测规划设计研究院有限公司 | The method for accurately generating segmental arc geologic section |
| CN108875755B (en) * | 2018-05-28 | 2021-11-12 | 南京师范大学 | Automatic horizontal rock stratum extraction method based on curve parallel characteristics |
| CN112417077B (en) * | 2020-11-25 | 2022-06-10 | 首都师范大学 | Method and device for automatically simplifying geologic body boundary and electronic equipment |
| CN113539051B (en) * | 2021-06-24 | 2022-11-25 | 南京师范大学 | Geological map-based method and device for acquiring stratum boundary point-by-point rock stratum attitude |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2009141572A (en) * | 2007-12-05 | 2009-06-25 | Celsys:Kk | Method and program for creating mask image for frame image of cartoon |
| CN101635051A (en) * | 2008-07-25 | 2010-01-27 | 鸿富锦精密工业（深圳）有限公司 | Boundary element extracting method and computer system thereof |
| CN103942809A (en) * | 2014-05-12 | 2014-07-23 | 福州大学 | Method for detecting joint fissures in rock images |
-
2016
- 2016-05-18 CN CN201610333329.9A patent/CN106023197B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2009141572A (en) * | 2007-12-05 | 2009-06-25 | Celsys:Kk | Method and program for creating mask image for frame image of cartoon |
| CN101635051A (en) * | 2008-07-25 | 2010-01-27 | 鸿富锦精密工业（深圳）有限公司 | Boundary element extracting method and computer system thereof |
| CN103942809A (en) * | 2014-05-12 | 2014-07-23 | 福州大学 | Method for detecting joint fissures in rock images |
Non-Patent Citations (2)
| Title |
|---|
| OPTIMISING THE APPLICATION OF THE HOUGH;N. C. FITTON 等;《Computers & Geosciences》;19981231;第22卷(第10期);全文 * |
| 人工读图机理分析及其在计算机读图中的应用;杨华飞 等;《计算机科学》;20081231;第35卷(第2期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106023197A (en) | 2016-10-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106023197B (en) | A kind of method of upright rock stratum automatic identification and extraction | |
| JP6000899B2 (en) | How to detect text automatically | |
| Pintus et al. | A survey of geometric analysis in cultural heritage | |
| Lee et al. | Linking past to present: Discovering style in two centuries of architecture | |
| Recky et al. | Windows detection using k-means in cie-lab color space | |
| CN106323836B (en) | A kind of borehole wall calculation of permeability and device | |
| AU2010321921B2 (en) | System and method for reservoir analysis | |
| CN104077447B (en) | Urban three-dimensional space vector modeling method based on paper plane data | |
| Patel et al. | The seismic analyzer: Interpreting and illustrating 2d seismic data | |
| CN107423455B (en) | Automatic modeling method and system for highway | |
| CN108876843A (en) | Method and system for component geometry structure extraction | |
| Satari et al. | A multi‐resolution hybrid approach for building model reconstruction from lidar data | |
| Al-Ruzouq et al. | Infrastructure growth assessment of urban areas based on multi-temporal satellite images and linear features | |
| TW201810093A (en) | User background information collection method and device | |
| CN106777117B (en) | Automatic identification method for horizontal rock stratum structure landform | |
| Wang et al. | A topology-preserving simplification method for 3d building models | |
| US11873709B2 (en) | Log based diagenetic rock typing and sweet spot identification for tight gas sandstone reservoirs | |
| CN110634222B (en) | Bank bill information identification method | |
| Đurić et al. | Two-dimensional shape analysis of complex geometry based on photogrammetric models of iconostases | |
| CN107066997A (en) | A kind of electrical equipment price quoting method based on image recognition | |
| CN115100394A (en) | City block function identification method based on interest point Voronoi graph | |
| CN114820960A (en) | Method, device, equipment and medium for constructing map | |
| CN107357923A (en) | The cubical method of the tax is generated based on FreeMarker | |
| Kostrikov et al. | Automated Extraction of Heavyweight and Lightweight Models of Urban Features from LiDAR Point Clouds by Specialized Web-Software | |
| CN107329181B (en) | A method of seeking muddy ore formation water resistivity and litho-electric parameters |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |