CN103093466A - Building three-dimensional change detection method based on LiDAR point cloud and image - Google Patents

Abstract

The invention provides a building three-dimensional change detection method based on a LiDAR point cloud and an image. Through LiDAR data filter and interpolation, a high-precision digital surface model (DSM) is generated, difference is carried out on the DSM and a candidate change area is obtained, then through a relevant formula, the candidate change area is enabled to correspond to a two-dimensional image, and change areas on the image are extracted through edge extraction, linear feature matching, edge difference, and the like. The building three-dimensional change detection method has the advantages of being capable of making full use of elevation information which is provided by LiDAR point cloud data, at the same time, combining information of the two-dimensional image to conduct the three-dimensional change detection on a manual building, and significantly improving precision, breadth and depth of the change detection.

Description

Buildings three dimensional change detection method based on LiDAR point cloud and image

Technical field

The invention belongs to the Surveying Science and Technology field, relate to a kind of buildings three dimensional change detection method based on LiDAR point cloud and image data, be mainly used in the fields such as artificial structure's extracting change information and processing.

Background technology

Change the data source that detects and remote sensing image arranged, orthography, map vector and LiDAR(laser radar) cloud data etc.Wherein, LiDAR is not subjected to the restriction of sunshine and weather condition, the round-the-clock earth observation of energy, and quick obtaining is low-cost, the three-dimensional coordinate on high-precision earth's surface, provides good data source for accurately extracting building feature.Especially, than additive method, the LiDAR cloud data can provide the artificial structure high-precision elevation information.But because airborne laser radar point cloud data lacks texture information, packing density skewness and be scrambling, discontinuously arranged in three dimensions simultaneously.Therefore directly utilize the LiDAR cloud data carry out artificial structure's three dimensional change detect more difficult, especially complicated at building shape, when the buildings edge details is more, be difficult to buildings is changed detection accurately and rapidly.The aviation image data space resolution of obtaining by aeroplane photography is high, and data distribute continuously simultaneously, have and extremely enrich texture information, but image data can not provide high-precision elevation information.Therefore LiDAR cloud data and image data are detected the mutual supplement with each other's advantages that can reach the different pieces of information source in conjunction with the three dimensional change of carrying out buildings.The control fusion of two kinds of data sources changes detection for artificial structure's high-level efficiency, high accuracy three-dimensional, to satisfy practical application, very important researching value is arranged for the technical development of LiDAR Data Post.Existing change detecting method has: (1) comprises the methods such as image difference, image ratio, image homing method, image vegetation index difference, Change vector Analysis and background subtracting based on the change detecting method of algebraic operation; (2) based on the change detecting method of image conversion, comprise the methods such as principal component analytical method (PCA), tasseled cap transformation method (KT), canonical correlation analysis; (3) based on the change detecting method of Images Classification, comprise the rear comparative approach of classification, multi-temporal image Direct Classification method (also claiming spectrum/time to classify mutually); (4) object-based change detecting method; (5) based on the change detecting method of statistical model; (6) based on the change detecting method of wavelet transformation.Although it is different that the variation that above-mentioned various change detecting method adopts detects elementary cell, the tactful difference of change detection is very large, but they all only use the plane information of bidimensional image, fully do not take elevation into account and change, and very important variation of artificial structure is with regard in performance variation in height.Therefore in the urgent need to seeking to consider the change detecting method of elevation variation.

Summary of the invention

The purpose of this invention is to provide a kind of artificial structure's three dimensional change detection method that elevation information changes of taking into account, it can overcome the deficiency of above-mentioned existing change detecting method technology, satisfies artificial structure's three dimensional change and detects the demand of using.

Technical scheme of the present invention is a kind of buildings three dimensional change detection method based on LiDAR point cloud and image, comprises the following steps,

Step 1 is carried out respectively filtering and interpolation processing to the LiDAR cloud data of two different times, generates the DSM of two different times;

Step 2 is asked for the difference of the DSM of two different times of step 1 gained, obtains DSM difference image;

Step 3 to the image processing of step 2 gained DSM difference, is obtained the candidate region of variation on DSM difference image;

Step 4 according to DSM and image coordinate corresponding relation, calculates the corresponding candidate change zone of candidate region of variation difference on two corresponding original aviation images of different times on step 3 gained DSM difference image;

Step 5 is carried out respectively Edge Gradient Feature to candidate change zone on two corresponding aviation images of different times of step 4 gained, obtains two breadths edge images;

Step 6 is carried out respectively extraction of straight line and coupling to step 5 gained two breadths edge images, the edge image after obtaining two width extraction of straight line and mating;

Step 7 is carried out the difference computing to the edge image after step 6 gained two width extraction of straight line and coupling, rejects simultaneously the edge outside unmatched straight line and linearity region, obtains the modified-image in the candidate change zone.

And, the described DSM of step 4 and image coordinate corresponding relation as shown in the formula,

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

Wherein, a ₁, a ₂, a ₃, b ₁, b ₂, b ₃, c ₁, c ₂, c ₃9 rotation matrix elements for the generation of image elements of exterior orientation; (X _s, Y _s, Z _s) be three outer orientation line elements of image; (X, Y, Z) is the ground point object coordinates on DSM; (x, y) is the picture side's coordinate on image.

And step 6 is described carries out respectively extraction of straight line and coupling to step 5 gained two breadths edge images, and the coupling implementation is as follows,

If two different times are designated as respectively T1 period and T2 period, the straight line that is located in the some candidate change zone of T2 image in period is expressed as S set _n, the straight line in the same candidate change zone on T1 image in period is expressed as S set _o,

For S set _nIn any straight line (ρ, θ), adopt following formula to calculate S set _oIn the distance of all straight lines,

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}$

Wherein, D represents S set _nUpper certain straight line (ρ, θ) and S set _oDistance between upper certain straight line (ρ ', θ ');

If calculate reckling in the gained distance less than certain given threshold value T, think that corresponding two straight lines of this distance mate, otherwise think S set _oIn do not have corresponding straight line and straight line (ρ, θ) to mate.

Advantage of the present invention is not only to take full advantage of the characteristics (straight line coupling) that images match detects, and takes full advantage of this important feature of artificial structure's height change.By the Treatment Analysis to the LiDAR data, utilize the DSM image that generates to detect more exactly artificial structure's variation.The method does not need manually to intervene in whole variation testing process, can satisfy in time the requirement of the fast database of renewal speed and Geographic Information System, and can carry out quantitative analysis, also can carry out to detected region of variation the judgement of qualitative change simultaneously.

Description of drawings

Fig. 1 is the buildings three dimensional change overhaul flow chart based on LiDAR data and image of the embodiment of the present invention.

Embodiment

Describe technical solution of the present invention in detail below in conjunction with drawings and Examples.

Referring to Fig. 1, the invention provides artificial structure's three dimensional change detection method, generate high precision DSM by LiDAR data filtering, interpolation, then DSM is carried out difference and obtain the candidate change zone, then by correlation formula, the candidate change zone is corresponded on bidimensional image, then by the region of variation on the extraction images such as edge extracting, linear feature coupling, edge difference.The input data are two phase airborne laser point clouds and two phase aviation images.Can adopt computer software technology to realize operation automatically during concrete enforcement.

The concrete methods of realizing of embodiment comprises the following steps:

Step 1, high precision DSM generates: the LiDAR cloud data to two different times carries out respectively filtering and interpolation processing, generates two width DSM(digital surface elevation models).

Embodiment processes the original LiDAR cloud data of new and old two different times, is designated as respectively T1 LiDAR in period data and T2 LiDAR in period data.At first two phase original LiDAR point clouds are carried out filtering with excluding gross error point (high point, utmost point low spot and noise spot), then the cloud data of excluding gross error point carried out triangulation network interpolation (TIN interpolation), adopt simultaneously certain scheme of colour, different colours represents different height value, and then generates the high precision DSM of two different times.The corresponding DSM of T1 LiDAR in period data and T2 LiDAR in period data is designated as respectively T1 DSM in period and T2 DSM in period.

Step 2 generates DSM difference image: ask for the difference of the DSM of two different times according to elevation, obtain DSM difference image.

Embodiment utilizes the variation of elevation, and the two phase high precision DSM that generate are carried out the difference computing, obtains DSM difference image.

Step 3, DSM candidate change zone generation: DSM difference image is carried out the processing means such as binary conversion treatment, edge tracking, obtain the candidate region of variation on DSM difference image.

Embodiment arranges suitable threshold value, DSM difference image is carried out binary conversion treatment obtain bianry image, then this bianry image is carried out the morphology border and follows the tracks of the candidate change zone that obtains on DSM difference image.

Step 4, aviation image candidate change territory generation: according to DSM and image coordinate corresponding relation, calculate the corresponding candidate change zone of candidate region of variation difference on two corresponding original aviation images of different times on step 3 gained DSM difference image.

Transform mathematical formulae as the formula (1) for ease of DSM and the image coordinate of implementing reference, provide to relate to:

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}---\left(1\right)$

Wherein, a ₁, a ₂, a ₃, b ₁, b ₂, b ₃, c ₁, c ₂, c ₃9 rotation matrix elements for the generation of image elements of exterior orientation; (X _s, Y _s, Z _s) be three outer orientation line elements of image; (X, Y, Z) is the ground point object coordinates on DSM; (x, y) is the picture side's coordinate on image.

Embodiment is by DSM in formula (1) and image coordinate projection conversion relation, calculate candidate region of variation corresponding zone on corresponding original aviation image of two phases in DSM, can obtain the candidate change zone of T1 image data in period and T2 image data in period.

Step 5, aviation image candidate change edges of regions is extracted: Edge Gradient Feature is carried out respectively in candidate change zone on two corresponding aviation images of different times of step 4 gained, obtain two breadths edge images.

Because the Canny operator can and suppress in accuracy of detection obtain good balance aspect noise, therefore embodiment adopts the Canny operator to carry out Edge Gradient Feature to two phase aviation image candidate change zones, obtain T1 edge image in period and T2 edge image in period.Canny operator specific implementation is prior art.

Step 6, edge image line feature extracting and matching: step 5 gained two breadths edge images are carried out respectively extraction of straight line and coupling, obtain the edge image (can referred to as coupling back edge image) after two width extraction of straight line and coupling.

Due to factors such as noise and uneven illuminations, the marginal point that has in edge image is discontinuous, must connect by the edge they are converted into significant edge.Embodiment uses the Hough conversion to extract the linear feature of edge image, then calculates two phase edge image rectilineal intervals to be matched from D according to formula (2), and the match is successful if D less than given threshold value T, thinks straight line, and it fails to match otherwise think straight line.Hough conversion specific implementation is prior art.

For the sake of ease of implementation, provide embodiment to use the Hough conversion to extract the linear feature in new and old image candidate change zone, and then carry out the linear feature coupling, concrete matching strategy is as follows:

The straight line that is located in the some candidate change zone of T2 image in period is expressed as S set _n:

s _n={(ρ ₀,θ ₀),(ρ ₁,θ ₁)......(ρ _n,θ _n)}

Wherein, ρ _i(i=0,1,2......n) in expression T2 image in period candidate change zone i bar straight line to the vertical range of initial point, θ _i(i=0,1,2 ... n) expression x axle is to the angle of the vertical line of i bar straight line, θ _i∈ [90 °, 90 °].

Straight line in same candidate change zone on T1 image in period is expressed as S set _o:

$S_o=\{({\mathrm{\ρ}}_0^{\′},{\mathrm{\θ}}_0^{\′}),({\mathrm{\ρ}}_1^{\′},{\mathrm{\θ}}_1^{\′})......({\mathrm{\ρ}}_n^{\′},{\mathrm{\θ}}_n^{\′})\}$

Wherein, ρ ' _i(i=0,1,2......n) in expression T1 image in period candidate change zone i bar straight line to the vertical range of initial point, θ ' _i(i=0,1,2......n) expression x axle is to the angle of the vertical line of i bar straight line, θ ' _i∈ [90 °, 90 °].

Set of computations S _nUpper certain straight line (ρ, θ) and S set _oDistance B between upper certain straight line (ρ ', θ '):

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}---\left(2\right)$

For S set _nIn any straight line, calculate S set _oIn the distance B of all straight lines, if the minor increment in all distances thinks that less than certain given threshold value T corresponding two straight lines of this distance mate, otherwise think S set _oIn do not have corresponding straight line and straight line (ρ, θ) to mate.

Step 7, the region of variation image generates: the edge image after step 6 gained two width extraction of straight line and coupling is carried out the difference computing, deletion the straight line that it fails to match and the edge outside the linearity region simultaneously obtain the modified-image in the candidate change zone, i.e. the region of variation image.

Specific embodiment described herein is only to the explanation for example of the present invention's spirit.Those skilled in the art can make various modifications or replenish or adopt similar mode to substitute described specific embodiment, but can't depart from spirit of the present invention or surmount the defined scope of appended claims.

Claims (3)

1. buildings three dimensional change detection method based on LiDAR point cloud and image is characterized in that: comprises the following steps,

Step 1 is carried out respectively filtering and interpolation processing to the LiDAR cloud data of two different times, generates the DSM of two different times;

Step 2 is asked for the difference of the DSM of two different times of step 1 gained, obtains DSM difference image;

Step 3 to the image processing of step 2 gained DSM difference, is obtained the candidate region of variation on DSM difference image;

Step 4 according to DSM and image coordinate corresponding relation, calculates the corresponding candidate change zone of candidate region of variation difference on two corresponding original aviation images of different times on step 3 gained DSM difference image;

Step 5 is carried out respectively Edge Gradient Feature to candidate change zone on two corresponding aviation images of different times of step 4 gained, obtains two breadths edge images;

Step 6 is carried out respectively extraction of straight line and coupling to step 5 gained two breadths edge images, the edge image after obtaining two width extraction of straight line and mating;

Step 7 is carried out the difference computing to the edge image after step 6 gained two width extraction of straight line and coupling, rejects simultaneously the edge outside unmatched straight line and linearity region, obtains the modified-image in the candidate change zone.

2. according to claim 1 based on the buildings three dimensional change detection method of LiDAR point cloud and image, it is characterized in that: the described DSM of step 4 and image coordinate corresponding relation as shown in the formula,

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

Wherein, a ₁, a ₂, a ₃, b ₁, b ₂, b ₃, c ₁, c ₂, c ₃9 rotation matrix elements for the generation of image elements of exterior orientation; (X _s, Y _s, Z _s) be three outer orientation line elements of image; (X, Y, Z) is the ground point object coordinates on DSM; (x, y) is the picture side's coordinate on image.

3. described buildings three dimensional change detection method based on LiDAR point cloud and image according to claim 1 and 2, it is characterized in that: step 6 is described carries out respectively extraction of straight line and coupling to step 5 gained two breadths edge images, and the coupling implementation is as follows,

If two different times are designated as respectively T1 period and T2 period, the straight line that is located in the some candidate change zone of T2 image in period is expressed as S set _n, the straight line in the same candidate change zone on T1 image in period is expressed as S set _o,

For S set _nIn any straight line (ρ, θ), adopt following formula to calculate S set _oIn the distance of all straight lines,

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}$

Wherein, D represents S set _nUpper certain straight line (ρ, θ) and S set _oDistance between upper certain straight line (ρ ', θ ');

If calculate reckling in the gained distance less than certain given threshold value T, think that corresponding two straight lines of this distance mate, otherwise think S set _oIn do not have corresponding straight line and straight line (ρ, θ) to mate.

Images