CN105180890A - Rock mass structural plane attitude measuring method integrating laser point cloud and digital image - Google Patents

Abstract

A full-automatic digital measuring method for rock mass structural plane attitude by fusing laser point cloud and digital image. The contour line of the track is extracted by adopting a hybrid global threshold value method and a local threshold value method by using a photogrammetric technology, the joint track is extracted by skeleton extraction and track mark, and the information acquisition of the rock joint track attitude is realized by combining other image processing methods; and finally, judging the position relation between the trace line and the structural surface according to the distance between the trace line and the structural surface, merging and grouping the trace line and the structural surface, and representing the occurrence of the rock structural surface by the occurrence of each group. The invention integrates the measuring methods of the laser point cloud and the digital image, fully exerts the properties of rapidness, high efficiency, non-contact and invisible terrain of the two measuring methods, simultaneously improves the measuring accuracy, and has the advantages of high automation degree, strong universality, comprehensive expression and the like.

Description

Merge the ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method of laser point cloud and digitized video

Technical field

The present invention relates to a kind of rock mass structure occurrence Full-automatic digital measuring method merging laser point cloud and digitized video, belong to engineering geological investigation technical field.

Background technology

ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method is summed up and mainly contains four classes:

1., namely scan line method, window statistic law measure structural plane occurrence information one by one by tape measure and compass manual site.The method is in the feature that rock mass structure EDS maps is wide, quantity is many, randomness is strong in Practical Project, poor efficiency, effort, consuming time making it are difficult to meet engineering demand, and run into the impact of slope, uninhabited cliff danger by field condition, measuring method is just difficult to carry out, thus causes the comprehensive of measurement data and representativeness to be also restricted.

2., hole directed core sample method (Rosengren, 1968), photographic process (Eoek & Pentz in hole, 1968), digital panoramic borehole camera method (Wang Chuanying, 2001) namely by after boring and coring to the measurement of structural plane occurrence contained in core or carry out digital picture by deep dynamic consolidation to hole wall inside structure face texture and measure the structural plane occurrence information obtaining and intersect within the scope of drilling depth.The method can obtain deep rock structural surface occurrence information, but pore diameter range measurement size makes the method lower for the measuring accuracy in large-scale structure face, and the broken rock of boring simultaneously and position, hole verticality also can have influence on measuring accuracy.

3., photogrammetry (Ross-Brown & At-kinson, 1972; Hagan, 1978) namely based on digital picture and photogrammetric ultimate principle, the multidisciplinary theory and means such as appliance computer 3 Dimension Image Technique, Image Matching, image interpolation, pattern-recognition merges mutually, the relative coordinate space geometry information on atural object surface can being obtained, can obtain exposing earth's surface ORIENTATION OF DISCONTINUITY IN ROCK MASS information by calculating these geological informations.The method adopts non-contacting measurement means, can not by the restriction of site contour condition, but the method precision controls by photographs quality and binocular camera shooting spacing, still need certain manual intervention simultaneously, the work efficiency of wall scroll rock mass discontinuity increases than accurate measurement collimation method, and a large amount of manual intervention of numerous structural plane remains unacceptable in Practical Project.

4., Three Dimensional Ground laser scanning method (TerrestrialLaserScanning, be called for short TLS), namely by can obtain the high precision cloud data (laser data) describing the relative geometric coordinate information in atural object surface fast to the pointwise laser scanning and ranging of ground object target, can be calculated those by certain manual intervention and data cutout and be exposed to ORIENTATION OF DISCONTINUITY IN ROCK MASS on crag surface.The method eliminates photogrammetric algorithm relatively loaded down with trivial details in photogrammetry, improve measuring speed and precision, but still need manual intervention Manual interception in mass cloud data to describe the cloud data exposing ORIENTATION OF DISCONTINUITY IN ROCK MASS, and the method is for not exposing rock mass discontinuity or being to measure with the ORIENTATION OF DISCONTINUITY IN ROCK MASS that trace form manifests.

In sum, the features such as rock mass discontinuity quantity is many in the Practical Project, distribution is wide and randomness is strong, current technical bottleneck is how to take into account contactless (cancelling place restriction), full-automatic (cancellation manual intervention), measuring speed fast (algorithm is succinct), universality (being applicable to multiple Environmental effect rock mass discontinuity).

Summary of the invention

Technical matters: the object of the invention is to get by the fast accurate space measurement means of Three Dimensional Ground laser scanner technique and the target information robotization sieve of computer visual image recognition technology, data fusion and geometrical registration is carried out by sieving the Digital Image Data after getting containing precise geometrical information laser point cloud data and target information robotization, propose a kind of full-automatic, contactless rock mass structure occurrence digital measuring method based on TLS three-dimensional laser scanner, namely merge the ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method of laser point cloud and digitized video.

Technical scheme: the features such as the quantity mainly for field rock mass discontinuity is many, distribution wide and randomness is strong, traditional ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method work efficiency is low, even runs into slope, uninhabited cliff danger and cannot implement at all.The accurate space measurement technology of the present invention by means of ground laser and the Computer Vision Recognition technology of digitized video, propose a kind of full-automatic, contactless ORIENTATION OF DISCONTINUITY IN ROCK MASS digital measuring method.

Merge the ORIENTATION OF DISCONTINUITY IN ROCK MASS Full-automatic digital measuring method of laser point cloud and digitized video, it measures the extraction that content includes the cloud data extraction of ORIENTATION OF DISCONTINUITY IN ROCK MASS, the view data of ROCK MASS JOINT trace occurrence merges mutually with cloud data, and its form of presentation adopts trace and structural plane matched packet, the common method stated.

In the extraction of cloud data, first carry out denoising with filter method, again triangulation carried out to cloud data thus set up preliminary Three-dimensional Rock model, give full play to the feature that cloud data degree of accuracy is high afterwards, structural plane is tentatively divided into groups based on the DBSCAN clustering method in high density UNICOM region by constructing eigenmatrix to the calculating of triangular facet, triangular facet normal vector and adopting, connectedness subsequently by triangular facet differentiates the secondary grouping carrying out structural plane, finally calculates the occurrence of structural plane; Rock joint trace need extract from view data, first gray processing process is carried out to digitized video, the outline line of rock is extracted again by image blend thresholding method, peel off and remove through pixel the framework information that amputation extracts trace, finally trace be marked and carry out occurrence calculating; After the strike-dip survey of structural plane and trace completes, trace and structural plane are rationally divided into groups, carry out differentiation merger by the position relationship that calculates between trace and structural plane, characterize the occurrence information of rock mass discontinuity with the net result after merger.

Owing to being comparatively difficult by the occurrence of cloud data extraction trace, therefore fusion image data of the present invention and cloud data carry out the extraction of trace occurrence.Rock joint trace need extract from view data, the present invention first carries out gray processing process to digitized video, the outline line of rock is extracted again by image blend thresholding method, peel off and remove through pixel the framework information that amputation extracts trace, finally trace be marked and carry out occurrence calculating.

After the strike-dip survey of structural plane and trace completes, for stating ORIENTATION OF DISCONTINUITY IN ROCK MASS more comprehensively, exactly, trace and structural plane rationally divide into groups by the present invention, carry out differentiation merger by the position relationship that calculates between trace and structural plane, characterize the occurrence information of rock mass discontinuity with the net result after merger.

The ORIENTATION OF DISCONTINUITY IN ROCK MASS Full-automatic digital measuring process merging laser point cloud and digitized video mainly comprises the steps:

The first step: collection site data;

Second step: cloud data and image data are carried out mating, merging;

3rd step: carry out triangulation to cloud data, calculates the parameter of subdivision Vee formation shape method phasor, carries out cluster analysis to normal vector;

4th step: the result after cluster analysis is merged, the data after being combined carry out occurrence calculating;

5th step: to correct image data, gray processing, then extracts the outline line at joint and skeleton and mark skeleton in a cloud;

6th step: calculate any 3 of trace corresponding point cloud and form leg-of-mutton interior angle, occurrence differentiation is carried out to trace and carries out occurrence calculating;

7th step: the structural plane that approximate fits extracts the distance of computation structure face and each trace, then divides into groups structural plane and trace according to the position relationship of structural plane and trace;

8th step: surveying work completes.

Beneficial effect: the features such as the present invention is many mainly for the quantity of field rock mass discontinuity, distribution wide and randomness is strong, traditional ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method work efficiency is low, even runs into slope, uninhabited cliff danger and cannot implement at all.The accurate space measurement technology of the present invention by means of ground laser and the Computer Vision Recognition technology of digitized video, propose a kind of full-automatic, contactless ORIENTATION OF DISCONTINUITY IN ROCK MASS digital measuring method.

In the process of laser point cloud data process, the present invention adopts the normal vector clustering algorithm of density based to extract main structural plane occurrence information with the form of robotization, and is divided into groups by structural plane by the connectedness differentiation of TIN intermediate cam shape.In the process, do not have existing K value cluster, fuzzy clustering etc. to need the shortcoming of artificial supposition cluster result group numbers, data processed result, by the impact of manual intervention, therefore robotization can identify the rock mass discontinuity of various shape.

In the process of digital image data process, the present invention adopts the mixing overall situation and local thresholding method to extract the outline line of trace, and is gone out the occurrence information of trace by skeletal extraction, trace marker extraction.Through practice test, the mixing overall situation and local thresholding method can reach degree of precision in the outline line extraction of ROCK MASS JOINT trace, decision condition that skeletal extraction adopts also is applicable to trace process completely, and trace mark subsequently more can play digital image data abundant information and the high advantage of laser point cloud data precision simultaneously.

In sum, the present invention is prosthetic intervention in strike-dip survey process, can obtain high-precision ORIENTATION OF DISCONTINUITY IN ROCK MASS information in conjunction with a series of digital processing method after image data, be a kind of full-automatic, contactless ORIENTATION OF DISCONTINUITY IN ROCK MASS digital measuring method.

Accompanying drawing explanation

Fig. 1 is the technical scheme process flow diagram of this measuring method.

Embodiment

The ORIENTATION OF DISCONTINUITY IN ROCK MASS Full-automatic digital measuring method of fusion laser point cloud of the present invention and digitized video, mainly comprise on-site data gathering, data fusion and registration, cloud data feature extraction, the occurrence acquisition of information of some cloud part, image feature extract, the occurrence acquisition of information of image part, laser point cloud and digitized video merge, concrete steps are as follows:

The first step: on-site data gathering

1) instrument in-site installation

Light may be caused according to measured zone position and size, structure characteristics of rock mass and complex engineering environment to block and the impact of reflecting destination object, determine height, scanning angle that the website quantity of TLS three-dimensional laser scanner erection measuring point, survey station point particular location and instrument set up, and according to above-mentioned requirements assembling and erection TLS three-dimensional laser scanner, if desired the scan-data of different website is spliced, also need to determine rational control points layout position, to carry out registration to the cloud data of different website;

2) field scan is measured

According to the rock mass structure integrality intending measurement zone, select suitable scanning resolution, correction parameter, the scanning angle of TLS three-dimensional laser scanner.According to instrument decorating position apart from distance and the light of surveying district, determining the parameters such as the exposure of high-resolution digital camera, aperture, aperture time, focal length, obtaining quality digital view data to guaranteeing;

After setting parameter, carry out laser scanning with scanner to needing the region of scanning, thus obtain cloud data, target object range scanner is far away, and the interval of target object surface analyzing spot is larger, suitably adjusts scanning accuracy as required;

While laser scanning, utilize the high-resolution digital camera being placed in scanner top to carry out the image acquisition of structure characteristics of rock mass, thus obtain view data;

Second step: Registration of Measuring Data and fusion

1) camera internal parameter acquiring

The digital camera that carries of three-dimensional laser scanner is generally fixed focal length camera, its focal length can directly obtain from camera lens mark, the inner parameter of digital camera is only relevant with camera lens with camera internal structure, can obtain by the shop instructions of camera and camera lens and Technical Specification Sheet;

2) camera external parameter obtains

When utilizing three-dimensional scanner system image data, camera rotates along with scanner and rotates, therefore need when carrying out camera external parameter and determining to consider the rotation matrix of camera relative to initial coordinate system, beam scanner system can relate to four coordinate systems: camera coordinates system CMCS, scanner coordinate system SOCS, engineering project coordinate system PRCS, world coordinate system GLCS; Spatial digitizer manufacturer provides rotation of coordinate, translation correlation matrix parameter, mainly contains 1. Mounting matrix, the coordinate transformation parameter matrix between camera coordinates system and scanner coordinate system; 2. COP matrix, the rotation matrix of camera moment camera coordinates system and initial camera coordinate system; 3. SOP matrix, each scanning movement point coordinate system is equivalent to the rotation translation matrix of engineering coordinate system;

3) the mutual conversion of coordinate system

1. world coordinate system and camera coordinates system

World coordinates P _w(X _w, Y _w, Z _w), camera coordinates P _u(X _u, Y _u, Z _u), both transformational relations:

$\left[\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]=M\·{\mathrm{COP}}^{-1}\·{\mathrm{SOP}}^{-1}\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]---(2-1)$

Wherein, rotation matrix t is the position coordinates of world coordinate system initial point in camera coordinates system, and rotation matrix R is orthogonal rotation matrix, meets:

$\left\{\begin{array}{c}{r}_{11}^2+r_{12}^2+r_{13}^2=1\\ r_{21}^2+r_{22}^2+r_{23}^2=1\\ r_{31}^2+r_{32}^2+r_{33}^2=1\end{array}\right.---(2-2)$

2. camera coordinates system and image coordinate system

Camera coordinates P _u(X _u, Y _u, Z _u), image coordinate is (x, y), both transformational relations:

$Z_u\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right]---(2-3)$

Wherein, f is camera focus.

3. world coordinate system and image coordinate system

World coordinates P _w(X _w, Y _w, Z _w), image coordinate is (x, y), both transformational relations:

$Z_u\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right]=\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]M\·{\mathrm{COP}}^{-1}\·{\mathrm{SOP}}^{-1}\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]---(2-4)$

3rd step: some cloud feature information extraction

1) cloud data denoising

The present invention adopts general median filtering method to carry out denoising to cloud data, and medium filtering ratio juris is exactly carry out the coordinate figure of alternative point cloud impact point by the intermediate value of all point coordinate in some cloud target neighborhood of a point,

2) triangulation

Triangulation is method the most basic in topology, scattered some cloud subdivision can be countless curved line trangle by it, triangulation will be carried out with Delaunay Triangulation Method at this, the delaunay function in Matlab is used to realize the Delaunay subdivision of three-dimensional point cloud, the data that delaunay function exports have been each leg-of-mutton numbering and apex coordinates after dividing, and the usage of delaunay function in Matlab is as follows:

I, input x a little, y, z coordinate;

II, use delaunay function to carry out triangle division, obtain each leg-of-mutton ID numbering and apex coordinate:

Tri＝delaunay(x,y)；

III, call z coordinate, use trimesh or trisurf instruction to draw diagram of block.

3) leg-of-mutton normal vector calculates

Utilize the coordinate on three summits of subdivision Vee formation shape, calculate leg-of-mutton normal vector coordinate and leg-of-mutton centre of form coordinate, and calculate the included angle in normal vector and coordinate Z axis, X-Z face further _t, θ _t, the centre of form is apart from the distance r of true origin _t,

If certain leg-of-mutton three summit of calculating in Matlab are

A (x _t1, y _t1, z _t1), B (x _t2, y _t2, z _t2), C (x _t3, y _t3, z _t3), then

$\stackrel{\→}{AB}=(x_{T2}-x_{T1},y_{T2}-y_{T1},z_{T2}-z_{T1}),\stackrel{\→}{AC}=(x_{T3}-x_{T1},y_{T3}-y_{T1},z_{T3}-z_{T1}),$

Leg-of-mutton centre of form coordinate is then centre of form O is apart from the distance of true origin $r_T=\sqrt{{(\underset{k=1}{\overset{3}{\Σ}}{x}_{Tk}/3)}^2+{(\underset{k=1}{\overset{3}{\Σ}}{y}_{Tk}/3)}^2+{(\underset{k=1}{\overset{3}{\Σ}}{z}_{Tk}/3)}^2}---(3-1)$

If the normal vector of Δ ABC is normal vector coordinate can be asked by building following system of equations:

$\left\{\begin{array}{c}(x_{T2}-x_{T1})a_T+(y_{T2}-y_{T1})b_T+(z_{T2}-z_{T1})c_T=0\\ (x_{T2}-x_{T3})a_T+(y_{T2}-y_{T3})b_T+(z_{T2}-z_{T3})c_T=0\\ (x_{T3}-x_{T1})a_T+(y_{T3}-y_{T1})b_T+(z_{T3}-z_{T1})c_T=0\end{array}\right.---(3-2)$

Trigonometric function relation is utilized to try to achieve normal vector subsequently φ _t, θ _tvalue:

${\mathrm{\θ}}_T=\arccos \frac{c_T}{\sqrt{{a_T}^2+{b_T}^2+{c_T}^2}}$

4th step: namely occurrence acquisition of information puts cloud part

1) normal vector cluster analysis

1. data center

Data center is the necessary preparation work before cluster analysis, the method for data center in the present invention, is with single leg-of-mutton φ _t, θ _t, r _tvalue deducts all triangle φ _t, θ _t, r _tthe mean value of value;

Cloud data is n triangle by subdivision altogether, so can calculate and obtain the coordinate (a of a corresponding n normal vector _ti, b _ti, 1) and (i=1,2,3...n) and φ _t, θ _t, r _tvalue φ _ti, θ _ti, r _ti(i=1,2,3....n); By the φ of n normal vector _t, θ _t, r _tbe worth with following rank, n × 3 matrix T statement, and retain ID in data processing:

$T=\left(x_{ij}\right)=\left[\begin{array}{ccc}{x}_{11}& x_{12}& x_{13}\\ x_{21}& x_{22}& x_{23}\\ \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}\\ x_{n1}& x_{n2}& x_{n3}\end{array}\right],$ φ _Ti＝x _i1,θ _Ti＝x _i2，r _T＝(x _i3),(i＝1,2,…n),(j＝1,2,3)

The mean value of note jth row is:

$\stackrel{\‾}{x_j}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_{ij}(i=1,2,...n),(j=1,2,3)---(4-1)$

So the centralization that n data object of jth row is implemented is transformed to:

${x^{,}}_{ij}=x_{ij}-{\stackrel{\‾}{x}}_j,(i=1,2...n),(j=1,2,3)---(4-2)$

After the conversion of above formula, often the average of row variable is 0, and namely often the value of row variable has identical basic point, then the result after central transformation is regenerated as new matrix T:

$T^{,}=\left({x^{,}}_{ij}\right)=\left[\begin{array}{ccc}{x^{,}}_{11}& {x^{,}}_{12}& {x^{,}}_{13}\\ {x^{,}}_{21}& {x^{,}}_{22}& {x^{,}}_{23}\\ \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}\\ {x^{,}}_{n1}& {x^{,}}_{n2}& {x^{,}}_{n3}\end{array}\right],(i=1,2,...n),(j=1,2,3).$

2. data normalization

Adopt minimax normalization method to carry out standardization to data at this, by a certain data in matrix divided by the maximal value of this column data absolute value, specifically see following formula:

${x_{ij}}^{,,}=\frac{{x_{ij}}^{,}}{{\max }_i\left|{x_{ij}}^{,}\right|}---(4-3)$

The matrix that final generation is new:

$T^{,,}=\left({x^{,,}}_{ij}\right)=\left[\begin{array}{ccc}{x^{,,}}_{11}& {x^{,,}}_{12}& {x^{,,}}_{13}\\ {x^{,,}}_{21}& {x^{,,}}_{22}& {x^{,,}}_{23}\\ \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}\\ {x^{,,}}_{n1}& {x^{,,}}_{n2}& {x^{,,}}_{n3}\end{array}\right]$

3. feature generates

Above-mentioned matrix T " belong to low-dimensional matrix, easy, vivid for describing, by T " each row vector be projected in three-dimensional cartesian coordinate system and form n point, and with (x " _i1, x " _i2, x " _i3) as this n point coordinate, be designated as (x _i, y _i, z _i); If the point having coordinate to repeat and T " in row vector numerical value identical because its ID is different, retained;

4. the cluster analysis of normal vector

DBSCAN clustering methodology is a kind of clustering algorithm based on high density UNICOM region, and this method is used and be have modified this clustering algorithm, thus reaches the object extracting structural surface information;

First, get ε=3, MinPts=5, T " in 2 q ₁(x _q1, y _q1, z _q2), q ₂(x _q2, y _q2, z _q2), wherein, parameter ε is field radius; With any for the center of circle, ε is the circle of radius, is called the ε field of this point; Parameter MinPts is lease core number of objects; If the quantity of the point comprised in the ε field of certain point is more than or equal to MinPts, then this point is called kernel object; In concrete engineering, the value of parameter ε, MinPts needs adjustment; Then, calculation level q ₁to a q ₂distance d, if d≤3, then will put q ₂be attributed to a q ₁ε field in,

$d=\sqrt{{(x_{q1}-x_{q2})}^2+{(y_{q1}-y_{q2})}^2+{(z_{q1}-z_{q2})}^2}---(4-4)$

Afterwards, input T " in coordinate a little, calculate and find out and can be attributed to a q ₁ε field in institute a little;

Suppose to obtain q ₁ε field in comprise m point be { q ₂, q ₃... q _m+1, then q ₁ε field in comprise count as m;

If m>=5, so by q ₁be classified as kernel object set A ₁interior point; If 0 < m < 5, then by q ₁be classified as edge object set B ₁interior point; If m=0, then by q ₁be classified as noise;

Again to q ₁field in other points { q ₂, q ₃... q _m+1whether be kernel object, will q be put by as above method is carried out calculating and differentiated ₁field in all points being classified as kernel object be attributed to set A ₁, all points being classified as edge object are attributed to set B ₁, delete the point being classified as noise;

Repeat input T " in all point coordinate, calculate the kernel object set { A of each point ₁, A ₂a _n, edge object set { B ₁, B ₂b _n, wherein the kernel object set of noise and edge object set are all designated as 0, then these two big collections are merged respectively as follows;

Get { A ₁, A ₂a _nin two set A ₁, A ₂if, set namely repeat a little, then by A ₂inside not with A ₁the sample intersected is incorporated to A ₁middle formation bunch A ₁', and by A ₁' participate in merging judgement as new combining objects next time; If set namely do not repeat a little, then by A ₁assignment is to a bunch A ₁' in, A ₂assignment is to a bunch A ₂' in, and by A ₁', A ₂' participate in merging judgement as new combining objects, then get { A. next time ₁, A ₂a _nmiddle A ₃respectively with A ₁', A ₂' carry out above-mentioned merging differentiation, if and then by A ₃assignment is to a bunch A ₃' in, and by A ₃' participate in merging judgement as new combining objects next time.Repetition like this, until { A ₁, A ₂a _nin all samples all carried out merge differentiate, finally obtain new set { A ₁', A ₂' ... A _v';

Edge object set { B ₁, B ₂b _nmerging also merge by mentioned above principle, obtain { B ₁', B ₂' ... B _w'.

{ A ₁', A ₂' ... A _v' and { B ₁', B ₂' ... B _w' be then the net result of cluster, in recycling cluster result, the ID numbering of each element just can extract the φ of correspondent method vector _t, θ _t, r value, also can mark corresponding tri patch;

2) cluster result process

1. tri patch merges

This step except invalid tri patch, will make up the deficiency of DBSCAN algorithm, and by the result extraction after cluster, sign, construct rock mass discontinuity by merging triangular facet.

I, utilize ID to mark tri patch corresponding in each bunch by the cluster result of DBSCAN, delete the tri patch repeated, and label is carried out to each bunch.

II, at a bunch A _h' (A _h' be { A ₁', A ₂' ... A _v' in any sample) in appoint and get a tri patch i, if its contiguous existence belongs to A _h' other tri patchs, then i is retained; Otherwise, delete.

III, the new bunch of A generated through II _h' in, appoint and get two tri patch ii, iii, if ii, iii are by belonging to A _h' in any number of tri patchs be connected, then ii, iii are included into structural plane X ₁; Otherwise ii is included into X ₁, iii is included into X ₂.So, by A _h' in tri patch carry out structural plane grouping;

IV, respectively by { A ₁', A ₂' ... A _v' in corresponding tri patch carry out II, III operation, obtain X ₁, X ₂x _s;

V, respectively by { B ₁', B ₂' ... B _w' in corresponding tri patch carry out II, III operation, obtain Y ₁, Y ₂y _r.

The X finally obtained ₁, X ₂x _s, Y ₁, Y ₂y _rform the tri patch of rock mass discontinuity exactly, each structural plane of rock mass after being marked by ID by these tri patchs, can be obtained;

2. excavation face is rejected

In Practical Project, contain a large amount of excavation face information in the occurrence information of rock mass discontinuity, and excavation face does not belong to the structural plane of rock mass.Because excavation face occurrence has repeatability, monistic feature, only need just can be identified in RiscanPro by the artificial method rejected and the tri patch of excavation face is deleted;

3. structural plane occurrence calculates

The occurrence information of structural plane comprises inclination angle and the tendency of structural plane, as long as to X ₁, X ₂x _sthe parameter in intermediate cam face calculates, and just can obtain the occurrence information of structural plane, specific algorithm is as follows:

Suppose X ₁corresponding H triangular facet, calculates X ₁normal vector coordinate (a in intermediate cam face _td, b _td, c _td), normal vector φ _t, θ _tvalue (φ _td, θ _td) (d and X ₁the ID of middle corresponding triangular facet is identical) arithmetic average:

$a_{X1}=\frac{1}{H}{\mathrm{\&Sigma;a}}_{Td}$

$b_{X1}=\frac{1}{H}{\mathrm{\&Sigma;b}}_{Td}---(4-5)$

$c_{X1}=\frac{1}{H}{\mathrm{\&Sigma;c}}_{Td}$

${\mathrm{\&phi;}}_{X1}=\frac{1}{H}{\mathrm{\&Sigma;\&phi;}}_{Td}$

${\mathrm{\&theta;}}_{X1}=\frac{1}{H}{\mathrm{\&Sigma;\&theta;}}_{Td}$

Then β _x1=| 90 ° of-θ _x1| (0 °≤β _x1≤ 90 °) be structural plane X ₁inclination angle;

Because tendency angular range is 0 ° ~ 360 °, be accurate Calculation tendency direction, need to utilize X ₁normal vector and Y, Z axis cosine value judge X ₁normal vector point to octant.Concrete method of discrimination is as follows:

Suppose X ₁normal vector be straight line L is X ₁with the intersection of X-Y plane in coordinate system, vector for horizontal projection vector, according to tendency direction definition, for structural plane X ₁tendency; α _x1, γ _x1be respectively with the angle of X, Z axis; k _x1for with the angle of Y-axis.Then have:

${\mathrm{cosk}}_{X1}=\frac{b_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2}}$

${\mathrm{cos\&alpha;}}_{X1}=\frac{a_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2+{c_{X1}}^2}}$

${\mathrm{cos\&gamma;}}_{X1}=\frac{c_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2+{c_{X1}}^2}}$

Calculate be structural plane X ₁tendency;

Repeat above-mentioned steps, bring the parameter in remaining structure face into, just can calculate inclination angle and the tendency of all structural planes;

4. completion structural plane

The disappearance putting cloud, tri patch on a small quantity may be caused in the sequence of operations such as collection, process, affect final bandwagon effect, imported in the softwares such as Geomagic, Pointcloud by the tri patch obtained and carry out completion, namely the contoured surface corresponding to each bunch finally obtained can be considered rock mass discontinuity;

5th step: image data feature extraction

1) imagery optical is corrected

According to digital camera imaging mathematical model, correct inner parameter and external parameter that parameter is divided into digital camera, formula 5-1 is digital camera imaging mathematical model formula,

$Z_c\left[\begin{array}{c}U\\ V\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{d_x}& 0& u_0\\ 0& \frac{1}{d_y}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\\ 1\end{array}\right]---(5-1)$

Order:

$M_1=\left[\begin{array}{ccc}\frac{1}{d_x}& 0& u_0\\ 0& \frac{1}{d_y}& v_0\\ 0& 0& 1\end{array}\right]M_2=\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]M_3=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]$

Wherein, Z _cfor photo imaging plane is relative to the distance of camera lens; U and V is the coordinate of Photograph image in imaging plane under pixel coordinate system; X _w, Y _w, Z _wfor the coordinate of real world object under overall world coordinate system; M ₁, M ₂for the inner parameter of digital camera, inner parameter is only relevant with camera lens with camera internal structure, can obtain by the shop instructions of camera and camera lens and Technical Specification Sheet; M ₃for the external parameter of digital camera, external parameter can be tested by camera calibration and be obtained through the inverting of camera imaging mathematical model formula;

After obtaining correcting parameter, can correct image picture according to formula 5-1, ripe business software or the subsidiary kit software of three-dimensional laser scanner generally can be adopted to carry out.

2) image gray processing

Captured by digital camera, the original state of digital picture is coloured image, first needs to change into gray level image to it.Gray level image refers to the digital picture only containing monochrome information, and brightness value change continuously.High resolving power digital image is converted to gray-scale map and digitizing, in fact exactly image is converted into a gray-scale value matrix F (M, N).Namely show that image size is the expression matrix of M × N number of pixel by M × N rank, in matrix, each value is expressed as the gray-scale value of pixel cell.Formula 5-2 is the expression-form of gray-scale value matrix.

3) outline line extracts

The method that outline line extracts adopts the image segmentation in Computer Image Processing method.The present invention proposes a kind of the mixing overall situation and local threshold method of applicable rock joint outline line extraction.

The first step utilizes the gray-scale value in rock joint district to be generally the feature of region local minimum, adopt local threshold method, test through a series of threshold values, min (the x in the threshold method of local can be obtained, y) be with (x, y) the Local Minimum gray-scale value in 7 × 7 pixel graticule mesh neighborhoods centered by, can reduce the noise in joint district like this; A (x, y) is made to be average gray value in 70 × 70 pixel graticule mesh neighborhoods centered by (x, y).First, the pixel compartments site of the min (x, y) that satisfies condition >=a (x, y) is found out in original image by local threshold method.Second step makes f (x, y)=0 of the pixel compartments site that satisfies condition form middle transition image.3rd step adopts the Otsu method of overall threshold method to determine segmentation threshold T to middle transfer image acquisition _o, and then carry out Iamge Segmentation to extract rock joint outline line.Otsu method finds threshold value according to statistical theory, and the best threshold values of Otsu method is determined by the inter-class variance maximal value of background image and target image.If the summation of pixel is N in image, the pixel on gray level l adds up to n ₁, N and n ₁relation represent such as formula 5-3, the probability distribution p of image histogram pixel _lmeet formula 5-4, the optimal threshold T of Otsu method _ofinally obtained by formula 5-5.

$N={\mathrm{\&Sigma;}}_{l=0}^{l-1}{n}_l---(5-3)$

$p_l=\frac{n_l}{N},{\mathrm{\&Sigma;}}_{l=0}^{l-1}{p}_l=1---(5-4)$

$T_O=\underset{0<l<1}{\mathrm{argmax}}\left\{{\mathrm{\&sigma;}}_B^2\left(l\right)\right\}---(5-5)$

In formula, for inter-class variance.

4) image denoising,

Picture noise can produce very large interference to follow-up at rock joint skeleton and topological relation extraction, so the denoising of reply rock joint contour line picture.

Concrete grammar is as follows:

I, binary conversion treatment is carried out to rock joint contour line picture;

II, the blackspot in the rock joint district in rock joint contour line picture or the hickie in rocky zone is eliminated.Use dilation operation operator in digital image morphology.Original image A is expanded by structural element B and may be defined as: be still made up of after reflection translation X pixel of structural element B all some x of intersection point with A.Namely operational formula is:

The present invention uses line structure element to carry out dilation operation to joint diagram picture.

III, the joint district after expansion process has good connectedness, uses the Bwareaopen function in Matlab software, filters out unwanted small size part in image according to connection district size;

IV, after the IIIth step process, the edge of joint outline line has some irregular tiny burrs, has an impact to follow-up joint skeletal extraction, so need to carry out the smooth of the edge to irregular tiny burr.Median filtering method is used to carry out smooth treatment herein.Medium filtering principle is: D given numerical value { a ₁, a ₂... a _dordered arrangement by size, when D is odd number, that numerical value being positioned at centre position is called as the intermediate value of this D numerical value; When D is even number, the mean value for those two numerical value in centre position is the intermediate value of this D numerical value, is denoted as med{a ₁, a ₂... a _d, in neighborhood window, the gray scale intermediate value of all pixels is as the gray-scale value of pixel filtered in window.That is: image is [x (I, J)] _{m × N}matrix, field window is A _d, after medium filtering, image output y (I, J) of pixel x (I, J) is designated as:

$y(I,J)={\mathrm{med}}_{A_d(I,J)}\&lsqb;x(I,J)\&rsqb;---(5-7)$

The square window of 5 × 5 is used to carry out medium filtering to image herein,

5) joint trace skeletal extraction

1., image thinning joint outline line can not be expressed rock joint framework characteristic intuitively and flutter and be opened up structure, needs further to extract joint trace skeleton, and extraction work is divided into two steps:; 2., disconnected branch is removed;

1., image thinning is exactly the pixel in bianry image is carried out pixel layer by layer on the basis keeping original shape and connectedness peel off, until residual image skeleton, image framework is the axis of image, and after image thinning, image residual pixel width is l;

Stripping conditions is as follows:

★ judges that whether meeting white point number in point 8 neighborhoods in binary map is 2 to 6;

It is all continuous that ★ judges whether meet white point in point 8 neighborhoods in binary map;

On ★ judges whether the point in binary map meets, left, the right side is not white point entirely;

On ★ judges whether point in binary map meets, left, under be not white point entirely;

2., remove disconnected branch after image thinning, still to there is the burr branch that breaks exactly need rejecting.Remove disconnected branch to need to identify the end points of joint skeleton line and point of crossing;

End points identification: make black picture element point value be 0 in bianry image; White pixel point value is 1, adopt the square detection window iteration of 3 × 3 to travel through whole image, check that all pixel point values are the pixel of 1, to check in its neighborhood pixel point value in 8 consecutive point to be the some number of 1, if number is 1, then inspection obtains this pixel is end points;

Junction recognition: point of crossing exists three wire cross or four line crunodes.The first step adopt the square detection window inspection of 3 × 3 traversal pixel point value be in the neighborhood of pixel points of 0 in 8 consecutive point pixel point value be 1 some number be 3 or 4; It is 5 × 5 squares that second step expands detection window, adds up the some number that the outermost pixel point value of 3 × 3 window is 1, if with the pixel point value of 3 × 3 window inner side be 1 some number identical, then to obtain this pixel be point of crossing in inspection;

Remove disconnected branch: Iamge Segmentation is become connected region independent separately with point of crossing by application endpoints, calculates the size of each connected region, i.e. pixel number N in connected domain _i, pixel number threshold value N in connected domain is set _tif, N _i< N _t, be just used as disconnected branch, deleted;

6) joint trace mark

Joint trace Information Statistics need to identify each rock joint trace, and method is as follows:

I, the connectedness of searching image pixel point set, and with different colours mark, utilize each searching in image of the Bwlabel function of Matlab software to be communicated with object, and according to the order searched out different round valuess mark,

II, Bwlabel function return value is the mark matrix of double type, thus image display cannot be realized, next utilize the label2rgb function of Matlab software, the round values of specifying each to be communicated with district in mark matrix is defined as color matrix, shows with rgb value image;

6th step: occurrence acquisition of information and image part

According to trace mark and Registration of Measuring Data, obtain the coordinate of each mark trace origin coordinates point and arbitrary intermediate point, suppose that the origin coordinates point of trace and the coordinate of arbitrary intermediate point are respectively E (x _e, y _e, z _e), F (x _f, y _f, z _f), G (x _g, y _g, z _g), a triangle can be determined at 3, three corresponding interior angle θ can be obtained _e, θ _f, θ _g:

Note:

$\stackrel{\&RightArrow;}{EF}=\stackrel{\&RightArrow;}{m}=(x_F-x_E,y_F-y_E,z_F-z_E)=(a_G,b_G,c_G)$

$\stackrel{\&RightArrow;}{EG}=\stackrel{\&RightArrow;}{n}=(x_G-x_E,y_G-y_E,z_G-z_E)=(a_F,b_F,c_F)$

$\stackrel{\&RightArrow;}{FG}=\stackrel{\&RightArrow;}{p}=(x_G-x_F,y_G-y_F,z_G-z_F)=(a_E,b_E,c_E)$

By ${\mathrm{cos\&theta;}}_E=\frac{\stackrel{\&RightArrow;}{m}\&CenterDot;\stackrel{\&RightArrow;}{n}}{\left|\stackrel{\&RightArrow;}{m}\right|\left|\stackrel{\&RightArrow;}{n}\right|}$ :

${\mathrm{\&theta;}}_E=\arccos \frac{\stackrel{\&RightArrow;}{m}\&CenterDot;\stackrel{\&RightArrow;}{n}}{\left|\stackrel{\&RightArrow;}{m}\right|\left|\stackrel{\&RightArrow;}{n}\right|}=\arccos \frac{a_G{a}_F+b_G{b}_F+c_G{c}_F}{\sqrt{a_G^2+b_G^2+c_G^2}\sqrt{a_F^2+b_F^2+c_F^2}}---(6-1)$

Using formula 6-1 in like manner can obtain θ again _f, θ _g

Note θ _j=max{ θ _e, θ _f, θ _g}

If 1. θ _j≤ 120 °, then on this trace, 3 institute's structure triangles can be considered to show up, and use formula 6-2, can ask this normal vector that goes out to show up

$\left\{\begin{array}{c}{a}_E{a}_J+b_E{b}_J+c_E{c}_J=0\\ a_F{a}_J+b_F{b}_J+c_F{c}_J=0\\ a_G{a}_J+b_G{b}_J+c_G{c}_J=0\end{array}\right.---(6-2)$

Normal vector { angle of cut of 0,0}, the inclination angle namely going out to show up, is expressed as with X-Y plane normal vector

Because tendency angular range is 0 ° ~ 360 °, be accurate Calculation tendency direction, need to utilize normal vector with Y, the cosine value of Z axis judges the octant that the normal vector of structural plane points to; Concrete method of discrimination is as follows: suppose that plane P is a structural plane, straight line L is the intersection of this structural plane and surface level, vector for the normal vector of plane P horizontal projection vector, according to tendency direction definition, for the tendency of structural plane; α _j, γ _jbe respectively normal vector with X, the angle of Z axis; k _jfor vector with the angle of Y-axis.Then have

${\mathrm{cosk}}_J=\frac{b_J}{\sqrt{{a_J}^2+{b_J}^2}}$

${\mathrm{cos\&alpha;}}_J=\frac{a_J}{\sqrt{{a_J}^2+{b_J}^2+{c_J}^2}}$

${\mathrm{cos\&gamma;}}_J=\frac{c_J}{\sqrt{{a_J}^2+{b_J}^2+{c_J}^2}}$

If 2. θ > 120 °, be then considered as common trace by this trace;

Utilize E, G point coordinate can obtain line EG equation: z=Qx+Py, appoint a bit (a in line taking _l, b _l, c _l)

The angle of line EG and coordinate plane X-Y, is trace inclination angle, is expressed as

Because tendency angular range is 0 ° ~ 360 °, be accurate Calculation tendency direction, need to utilize vector with Y, the cosine value of Z axis judges the octant that the normal vector of structural plane points to; Concrete method of discrimination is as follows: suppose vector vector for horizontal projection vector, according to tendency direction definition, for the tendency of trace EG; α _l, γ _lbe respectively vector with X, the angle of Z axis; k _lfor vector angle with Y-axis, then have

${\mathrm{cosk}}_L=\frac{b_L}{\sqrt{{a_L}^2+{b_L}^2}}$

${\mathrm{cos\&alpha;}}_L=\frac{a_L}{\sqrt{{a_L}^2+{b_L}^2+{c_L}^2}}$

${\mathrm{cos\&gamma;}}_L=\frac{c_L}{\sqrt{{a_L}^2+{b_L}^2+{c_L}^2}}$

7th step: laser point cloud and digitized video merge

1) approximate fits characteristic plane

Obtain each bunch of tri patch in work above, according to the apex coordinate that ID extracts triangular facet, suppose X ₁in certain Atria apex coordinate be (x _r1, y _r1, z _r1), (x _r2, y _r2, z _r2), (x _r3, y _r3, z _r3), so can try to achieve this leg-of-mutton centre of form coordinate (x according to following formula _r, y _r, z _r) characterize this triangular facet by the centre of form:

$x_R=\frac{x_{R1}+x_{R2}+x_{R3}}{3}$

$y_R=\frac{y_{R1}+y_{R2}+y_{R3}}{3}---(7-1)$

$z_R=\frac{z_{R1}+z_{R2}+z_{R3}}{3}$

Then simulate approximate characteristic plane by the method for least square fitting after calculating the centre of form coordinate with bunch tri patch, the usual shape of equation of space plane is as A _rx+B _ry+C _rz+1=0, goes out this plane with H the triangular facet centre of form coordinate fitting of same bunch and can be expressed as following matrix form:

$\left[\begin{array}{ccc}{x}_{R1}& y_{R1}& z_{R1}\\ \begin{array}{c}\&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\end{array}& \begin{array}{c}\&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\end{array}& \begin{array}{c}\&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\end{array}\\ x_{RH}& y_{RH}& z_{RH}\end{array}\right]\left[\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right]=\left[\begin{array}{c}-1\\ -1\\ -1\end{array}\right]$

Following form can be changed into through computing:

$\left[\begin{array}{ccc}{\mathrm{\&Sigma;x}}_g^2& {\mathrm{\&Sigma;x}}_g{y}_g& {\mathrm{\&Sigma;x}}_g{z}_g\\ {\mathrm{\&Sigma;x}}_g{y}_g& {\mathrm{\&Sigma;y}}_g^2& {\mathrm{\&Sigma;y}}_g{z}_g\\ {\mathrm{\&Sigma;x}}_g{z}_g& {\mathrm{\&Sigma;y}}_g{z}_g& {\mathrm{\&Sigma;z}}_g^2\end{array}\right]\left[\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right]=\left[\begin{array}{c}-\mathrm{\&Sigma;}{x}_g\\ -\mathrm{\&Sigma;}{y}_g\\ -{\mathrm{\&Sigma;z}}_g\end{array}\right](g=1,\mathrm{2...}H)$

So obtain: $\left[\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right]={\left[\begin{array}{ccc}{\mathrm{\&Sigma;x}}_g^2& {\mathrm{\&Sigma;x}}_g{y}_g& {\mathrm{\&Sigma;x}}_g{z}_g\\ {\mathrm{\&Sigma;x}}_g{y}_g& {\mathrm{\&Sigma;y}}_g^2& {\mathrm{\&Sigma;y}}_g{z}_g\\ {\mathrm{\&Sigma;x}}_g{z}_g& {\mathrm{\&Sigma;y}}_g{z}_g& {\mathrm{\&Sigma;z}}_g^2\end{array}\right]}^{-1}\left[\begin{array}{c}-\mathrm{\&Sigma;}{x}_g\\ -{\mathrm{\&Sigma;y}}_g\\ -{\mathrm{\&Sigma;z}}_g\end{array}\right](g=1,\mathrm{2...}H),$ Coefficient A _r, B _r, C _rcan try to achieve, plane equation can draw;

2) trace and structural plane merger

From cloud data, extract the coordinate of each point corresponding to trace, the distance of the almost plane a little simulated to each bunch is obtained in pointwise, gets the position relationship that minor increment can describe this trace and each clustering architecture face,

Point is to the algorithm of plan range: $D=\frac{Ax+By+Cz+1}{\sqrt{A^2+B^2+C^2}}---(7-2)$

Each bunch around exists 1 or many traces nearest with it, and this 1 or many traces and this bunch of merger can be entered same group, so the occurrence information of this group is the occurrence information of corresponding region.If exist around certain trace and there is not close structural plane, then this trace is classified as separately one group, in order to characterize the occurrence information of region; If exist around certain structural plane and there is not close trace, then this structural plane is classified as separately one group, in order to characterize the occurrence information of region.

Claims (2)

1. one kind merges the ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method of laser point cloud and digitized video, it is characterized in that the extraction that the measurement content of this measuring method includes the cloud data extraction of ORIENTATION OF DISCONTINUITY IN ROCK MASS, the view data of ROCK MASS JOINT trace occurrence merges mutually with cloud data, its form of presentation adopts trace and structural plane matched packet, jointly shows

The method stated;

In the extraction of cloud data, first carry out denoising with filter method, again triangulation carried out to cloud data thus set up preliminary Three-dimensional Rock model, give full play to the feature that cloud data degree of accuracy is high afterwards, structural plane is tentatively divided into groups based on the DBSCAN clustering method in high density UNICOM region by constructing eigenmatrix to the calculating of triangular facet, triangular facet normal vector and adopting, connectedness subsequently by triangular facet differentiates the secondary grouping carrying out structural plane, finally calculates the occurrence of structural plane; Rock joint trace need extract from view data, first gray processing process is carried out to digitized video, the outline line of rock is extracted again by image blend thresholding method, peel off and remove through pixel the framework information that amputation extracts trace, finally trace be marked and carry out occurrence calculating; After the strike-dip survey of structural plane and trace completes, trace and structural plane are rationally divided into groups, carry out differentiation merger by the position relationship that calculates between trace and structural plane, characterize the occurrence information of rock mass discontinuity with the net result after merger.

2. the ORIENTATION OF DISCONTINUITY IN ROCK MASS measuring method merging laser point cloud and digitized video as claimed in claim 1, is characterized in that this measuring method mainly comprises the steps:

The first step: collection site data;

Second step: cloud data and image data are carried out mating, merging;

3rd step: carry out triangulation to cloud data, calculates the parameter of subdivision Vee formation shape method phasor, carries out cluster analysis to normal vector;

4th step: the result after cluster analysis is merged, the data after being combined carry out occurrence calculating;

5th step: to correct image data, gray processing, then extracts the outline line at joint and skeleton and mark skeleton in a cloud;

6th step: calculate any 3 of trace corresponding point cloud and form leg-of-mutton interior angle, occurrence differentiation is carried out to trace and carries out occurrence calculating;

7th step: the structural plane that approximate fits extracts the distance of computation structure face and each trace, then divides into groups structural plane and trace according to the position relationship of structural plane and trace;

8th step: surveying work completes.