CN104484668A - Unmanned aerial vehicle multi-overlapped-remote-sensing-image method for extracting building contour line - Google Patents

Abstract

The invention discloses an unmanned aerial vehicle multi-overlapped-remote-sensing-image method for extracting a building contour line. The unmanned aerial vehicle multi-overlapped-remote-sensing-image method includes the steps that three-dimensional point cloud is generated with an aerial-triangulation and dense-matching combined method and filtered, and a building is detected from the point cloud; after the walls of the detected building are canceled, the general contour of the building is extracted from building top face information; the general contour of the building serves as a buffering area to be overlapped on spliced images and serves as shape prior information, evolution is carried out in the buffering area with a level set algorithm, and finally an accurate contour of the building is obtained. By means of the unmanned aerial vehicle multi-overlapped-remote-sensing-image method, as point cloud three-dimensional information generated by the multiple overlapped images is sufficiently used, and meanwhile the high-accuracy geometrical information of the high-resolution remote sensing images is used in a combined mode, the building contour extracting accuracy is remarkably improved, and the complexity of the method is lowered.

Description

Building contour line extraction method for unmanned aerial vehicle multi-overlapping remote sensing image

Technical Field

The invention relates to the technical field of remote sensing image application, in particular to a method for extracting a building contour line of an unmanned aerial vehicle remote sensing image.

Background

The building is an important geographic space element in a city, and has important significance in the application fields of city planning and management, city development and change, disaster detection and evaluation and the like. Building contour extraction is an important step in the establishment and updating of urban basic geographic information systems. As a novel remote sensing monitoring platform, the unmanned aerial vehicle has high flight operation intelligent degree, can autonomously fly and shoot according to a preset air route, provides remote sensing monitoring data and low-altitude video monitoring in real time, and has the characteristics of strong maneuverability, convenience, low cost and the like.

The high-resolution remote sensing image contains a large amount of abundant information, and the building outline extraction is often interfered by various other ground objects, such as building and non-building distinction, shielding of trees around the building, influence of road borders and the like. Therefore, the extraction of the building outline is performed only by a single image, and the technical difficulty is great. The extraction of the building outline not only needs to rely on the extraction and analysis of remote sensing two-dimensional information, but also needs to be combined with the building three-dimensional information, so that the two-dimensional information and the three-dimensional information are mutually fused and supplemented, and the extraction of the building outline in the remote sensing image is more facilitated. The typical method for extracting the building outline by using the high-resolution remote sensing image comprises the following steps: 1) and extracting the contour line of the building based on the single high-resolution remote sensing image. Although the remote sensing image with high resolution has clear building outline information, artificial buildings and non-buildings are difficult to distinguish, and in addition, tree sheltering around the buildings also generates certain interference on the outlines of the buildings, so the method has certain limitations. 2) Building contour extraction based on shadow assistance. Although building contour extraction with the aid of shadows indirectly utilizes height information of buildings, shadow extraction does not have certain universality, and parameters related to the requirement of obtaining the height of the buildings by utilizing the shadows are more, so that the method is difficult to meet the actual requirement. 3) And extracting the contour line of the building based on the Lidar and the remote sensing image. Although the method not only utilizes three-dimensional information of Lidar, but also integrates high-precision geometric outline information of images, and extracts the outline information of the building through mutual complementation of the advantages and disadvantages of two kinds of data. However, this method has difficulty in high-precision registration of Lidar and remote sensing images, and the Lidar data acquisition cost is expensive. 4) And extracting the contour line of the building based on the three-dimensional aerial image. Although this type of method obtains three-dimensional information by stereo matching, it extracts building contour information by the complementation of two types of information, while using high-precision two-dimensional information of images. However, the problem of this method is that the stereo image has a small size, which has a certain influence on extracting the contour of the buildings in the large area. Therefore, a method which is easy to obtain data, has a high automation degree of an extraction method, has a relatively high extraction result and meets the actual production needs is urgently sought.

Disclosure of Invention

The invention fully utilizes the point cloud three-dimensional information generated by the multiple overlapped images, and simultaneously combines the high-precision information of the remote sensing image, thereby obviously improving the precision of extracting the building outline.

The technical scheme of the invention is a building contour line extraction method of unmanned aerial vehicle multi-overlapping remote sensing images, which comprises the following steps:

the method comprises the following steps that firstly, adjustment is carried out on three pairs of unmanned aerial vehicle remote sensing images by utilizing an aerial vehicle, and meanwhile, the images are subjected to dense matching by utilizing a PMVS algorithm after GPU acceleration, and finally dense color point clouds with high precision are obtained;

splicing the remote sensing images of the unmanned aerial vehicle after the leveling;

step three, filtering the color point cloud; firstly, an improved morphological filtering algorithm is utilized to separate the ground from the non-ground, then, color invariants are utilized to filter vegetation in ground points, and finally, elevations and areas are utilized as thresholds to filter non-buildings;

detecting buildings in the point cloud by using a region growing algorithm;

deleting the wall surface of the building, and fitting the boundary of the top surface to finally obtain the rough contour information of the building;

step six, the coarse building contour obtained in the step three is used as a buffer area for building contour extraction and is superposed on the spliced image;

and step seven, simultaneously, utilizing the shape of the building rough contour as prior information, and evolving the building precise contour in the buffer area by using a level set algorithm.

The first step comprises the following steps:

(1.1) preprocessing the remote sensing image of the multi-view overlapping unmanned aerial vehicle by using prior information:

(1.2) carrying out aerial three-dimensional photogrammetry on the basis of the step (1.1), solving external orientation elements of each image by utilizing an aerial three-dimensional mesh, and carrying out integral adjustment of a beam method;

and (1.3) according to the image groups, on the basis of the step (1..2), performing rapid dense matching by using a PMVS (graphics processing unit) algorithm accelerated by a GPU (graphics processing unit) in the prior art to generate dense three-dimensional point cloud, and using the reconstructed point cloud as three-dimensional elevation data.

The second step comprises the following steps:

(2.1) feature extraction: extracting the features of the image by using SIFT;

(2.2) image registration: firstly, coarse registration is carried out, and matched feature points are searched by using a k-d tree; then, fine registration is carried out, and rough registration often has wrong matching points, so that the RANSAC algorithm is used for eliminating the wrong matching points; obtaining a transformation matrix between the images through the registration of the images;

(2.3) image stitching: splicing images through the transformation matrix obtained in the step (2.2);

(2.4) fusion of images: and after splicing, fusing the images by using a bilinear interpolation algorithm.

The third step comprises the following steps:

(3.1) separating the ground points and the non-ground points of the point cloud by using improved morphological filtering;

(3.2) filtering the vegetation by using the color invariant in the ground points in the step (3.1);

and (3.3) filtering out non-building target points by using a threshold value based on the characteristics of the building.

The step (3.1) comprises the following processes:

firstly, any point and neighborhood points thereof are taken to form a window with a fixed size, the lowest point in the window is detected through morphological opening operation, if the difference between the elevation value of the point in the window and the elevation of the lowest point is in a threshold range, the point is represented as a ground point, and all points in the point cloud are taken out for filtering;

then according to y2 xw^k+1, obtaining the size of a window required by next filtering, wherein the size of the window is smaller than the maximum value of a preset filtering window, and performing morphological filtering again; finally, when the window is larger than the preset window, finishing filtering; wherein k is the iteration number, and w is the size of the previous filtering window.

The step (3.2) comprises the following processes:

because the point cloud generated by image dense matching has color information, the green vegetation is filtered by utilizing the color invariant theory; let the coordinates of each point in the point cloud be (x, y, z), the color three channels be (R, G, B), and the threshold value of the color invariant for vegetation be T_gThe color invariant formula defined using the green and blue channels is:

<math> <mrow> <msub> <mi>&psi;</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>4</mn> <mi>&pi;</mi> </mfrac> <mo>&times;</mo> <mi>arctan</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein, I_g(x,y,z)、I_b(x, y, z) represents the green and blue component values of the point cloud at the (x, y, z) point;

ψ_g(x, y, z) represents the average of the values in (x,y, z) color invariant value of the dot; when psi_g＜T_gWhen the point is a vegetation point, the point is represented as a vegetation point; when psi_g＞T_gWhen it is, it means that the point is a non-vegetation point.

Fourthly, the point cloud in the third step is segmented by using a region growing method based on plane fitting to obtain a point cloud region V of each building_BiThe method specifically comprises the following steps:

(4.1) the points in the point cloud are divided into two types: if one point in the neighborhood of the point is default, the point belongs to the boundary point; otherwise such points belong to interior points;

(4.2) setting a certain interior point p₀(x₀,y₀,z₀),p₀The eight neighborhood points of (a): { p₁,p_2,…,p₈And fitting a plane equation (2) to the 9 points by using a least square method, and specifically realizing the following steps:

z＝ax+by+c (2)

for point p₀And 8 points within the neighborhood: (x)_i,y_i,z_i) I is 0,1,2, …, 7; fitting calculation formula (2) mean

A surface equation that minimizes equation (3);

<math> <mrow> <mi>S</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>7</mn> </munderover> <msup> <mrow> <mo>(</mo> <mi>ax</mi> <mo>+</mo> <mi>by</mi> <mo>+</mo> <mi>c</mi> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>

to minimize S, one should satisfy:

from this the following normal system of equations can be derived:

<math> <mrow> <mfenced open='|' close='|'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>x</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>y</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>n</mi> </mtd> </mtr> </mtable> </mfenced> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>b</mi> </mtd> </mtr> <mtr> <mtd> <mi>c</mi> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mi>i</mi> <mo>=</mo> <mn>0,1</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mn>7</mn> <mo>;</mo> <mi>n</mi> <mo>=</mo> <mn>8</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>

(4.3) finding the SSD and the variance of any interior point in the 8 neighborhood according to the formula (4):

<math> <mrow> <mi>SSD</mi> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>&Element;</mo> <mi>M</mi> </mrow> </munder> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein M is p₀And its set of eight neighborhood points, h_kAnd z_kRespectively obtaining an observation elevation value and a fitting plane elevation value;

(4.4) sequencing each point in the point cloud in an ascending order according to the SSD value, and taking the point with the smallest SSD value as a seed point;

(4.5) calculating the distance h from each point to the plane where the seed points are located in the neighborhood of the seed points, if h is<h_T,h_TIf the elevation threshold is the elevation threshold, merging the elevation threshold into the same area, and marking the area; otherwise, this is a non-patch point;

and (4.6) when all the neighborhood points are traversed, finding unprocessed points from the new SSD value as new seed points, and repeating the operation of (4.5) until all the points are traversed.

Step five, the point cloud area V of each building in the step four_BiDeleting wall points of the building, and finally obtaining rough contour information of the building by solving top surface boundary points, wherein the method specifically comprises the following steps:

(5.1) according to the divided surface patches in the step four, if the normal vector of the surface patch is parallel to the ground, the surface patch is indicated as a wall surface patch; if the normal vector of the patch is vertical to the ground, the patch is a roof patch; therefore, deleting the wall surface patches according to the normal vector;

(5.2) obtaining the boundary of the point cloud by utilizing an alpha-shape algorithm according to the top surface point cloud obtained in the step (5.1); the specific method comprises the following steps:

a) extracting an arbitrary point p1 from the top surface point cloud, searching a point set p2 with a distance of less than or equal to 2 alpha from the remaining points, wherein a is the radius of a circle, and setting p2 to { p21, p22, p23, … and p2n };

b) taking out a point p2i from any point p2, and finding out the circle center p0 of the points p1 and p2i by using the formulas (6) and (7); knowing two points (x1, y1), (x2, y2) and the radius α of the circle, the equation for the center of the circle (x0, y0) is as follows:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>1</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>1</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>2</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>2</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>

it is difficult to directly solve the equation, so the distance intersection algorithm in mapping is used to obtain:

$\left\{\begin{array}{c}x0=x1+\frac{x2-x1}{2}+H(y2-y1)\\ y0=y1+\frac{y2-y1}{2}+H(x2-x1)\end{array}\right.---\left(7\right)$

wherein,S²＝(x1-x2)²+(y1-y2)²；

c) finding the distances l from all points to p0 from the point set p2, if l > alpha, then p1 and p2i are boundary points; if l < alpha, the next step d) is carried out;

d) repeating the three steps a) b) c) for other points in p2 until all points in p2 are judged to be finished.

And step six, superposing the top surface point cloud contour points of each building obtained in the step five on the image to obtain a building buffer region omega on the image_Bi(ii) a The specific process is as follows:

setting a projection matrix P_3×4Computing P as a known inside and outside orientation element representing a single view geometry_3×4Matrix:

$P=\left[\begin{array}{ccc}-f& 0& x_0\\ 0& f& y_0\\ 0& 0& 1\end{array}\right]R^T\left[\begin{array}{cccc}1& 0& 0& -X_s\\ 0& 1& 0& -Y_s\\ 0& 0& 1& -Z_s\end{array}\right]---\left(8\right)$

wherein P is the projection matrix P_3×4F is the focal length of the image, x₀And y₀The eccentricity of the optical axis from the optical center in the horizontal direction and the vertical direction; x_s、Y_s、Z_sAs the coordinates of the camera center in the world coordinate system, R^TRepresents a 3 × 3 rotation matrix;

the projection formula is:

$\left[\begin{array}{c}\mathrm{xz}\\ \mathrm{yz}\\ z\end{array}\right]=P\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(9\right)$

calculating the projection of the point on the point cloud on the image by using the formula (8) and the formula (9); x, Y and Z represent the coordinates of the object point in the camera coordinate system, and X, Y and Z represent the coordinates of the object point in the world coordinate system.

Step seven, performing boundary connection on the points which are superposed on the image in the step six, superposing the obtained building rough contour on the spliced image as a buffer area, and simultaneously, utilizing the shape of the building rough contour as prior information and evolving the building precise contour in the buffer area by using a level set algorithm; the method comprises the following specific steps:

(7.1) carrying out boundary tracking on the points superposed on the image to obtain a closed contour; expanding the contour to obtain a buffer region omega for extracting the contour of the building on the image_Bi；

(7.2) performing local level set evolution in the remote sensing image in the buffer area by taking the contour of the image area obtained according to the step (7.1) as an initial level set and taking the contour shape of the image area obtained according to the step (7.1) as prior information to obtain a fine contour of a building of the remote sensing image; the specific implementation is as follows:

first, the buffer area of each building on the image is preset to be omega_BiThe shape of the contour obtained in (7.1) is [ phi ]_Pi，Ф_PiAs a prior shape, phi_SiAs a target to segment the contour, basisThe energy functional for shape constraints is defined as follows:

E_to＝E(c₁,c₂,φ_Si)+E_sh(φ_Si,φ_P1) (10)

the initial contour and the prior shape of the target are expressed by utilizing a moving target area:

<math> <mrow> <msub> <mi>E</mi> <mi>sh</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mo>,</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mo>&Integral;</mo> <msub> <mi>&Omega;</mi> <mi>Bi</mi> </msub> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>dxdy</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow> </math>

the energy function of the prior art based on the level set C-V model is:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <mi>E</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>,</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <mi>&mu;</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> <mi>v</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mi>dxdy</mi> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow> </math>

where H (φ) is a Heaviside function of the form: <math> <mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>1</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&GreaterEqual;</mo> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&lt;</mo> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> </math>

(phi) is a Dirac function of the formu₀(x, y) is the gray value of a certain point in the image area to be processed,modulus of the gradient of the current point, coefficient lambda₁，λ₂More than 0, mu, v ≥ 0 is a fixed parameter, and is generally selected as lambda₁＝λ₂1, μ ═ 0, v ═ 1, and (12) partial differential equations corresponding to the following equations:

<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>[</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> </mrow> </mfrac> <mo>-</mo> <mi>v</mi> <mo>-</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <msup> <mrow> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>]</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>]</mo> </mrow> <mn>2</mn> </msup> <mo>]</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein the parameter c₁、c₂Obtained according to the following formula,

<math> <mfenced open='' close=''> <mtable> <mtr> <mtd> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> </mfrac> </mtd> <mtd> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> </math>

according to the formula (13), a C-V level set evolution method is adopted to segment and extract the building region which is used as an image

Object v of building segmentation_iSegThe C-V level set is prior art.

Compared with the prior art, the invention has the following advantages and beneficial effects:

(1) although the data source is a single data source, two types of data information can be generated, and the effect of extracting the building outline by the two types of data sources is achieved. The complexity of the method is reduced, and the production cost is also saved.

(2) The three-dimensional information of the three-dimensional point cloud generated by the image is utilized to obtain the outline of the rough building, and meanwhile, the high-precision outline information of the high-resolution image is utilized, and the high-precision outline information and the high-resolution image are mutually supplemented and fused, so that the automation degree and the precision of the extraction of the outline of the building are improved.

(3) The building rough contour extracted by the three-dimensional point cloud information is used as the prior shape information of level set evolution and the buffer area for extracting the building contour, so that the speed and the precision of the level set evolution building contour are ensured.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 is a group chart of an image of point A in an oscillogram in step one.

Detailed Description

The invention provides a building contour line extraction method of an unmanned aerial vehicle remote sensing image, which comprises the steps of firstly utilizing a space-three algorithm in photogrammetry to combine with multi-view geometric stereo reconstruction in computer vision to quickly generate three-dimensional point cloud with geographic coordinates, then extracting rough contour information of a building through three-dimensional point cloud segmentation, establishing a buffer area for the rough contour information of the building to be superposed on a high-resolution image, then utilizing the rough contour as prior shape information, and utilizing a level set to carry out evolution iteration in the buffer area to obtain the precise geometric contour of the building. Since the three-dimensional point cloud is generated by image dense matching, there is no registration difficulty between the point cloud and the image. The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments, wherein a flow chart is shown in fig. 1, and the technical solution flow of the embodiments includes the following steps:

the method comprises the steps of firstly, carrying out adjustment on remote sensing images of the unmanned aerial vehicle by utilizing the air-to-three pairs, and carrying out dense matching on the images by utilizing a PMVS algorithm after GPU acceleration to finally obtain dense color point clouds with high precision and geographic coordinates. The specific implementation is as follows:

(1) preprocessing the remote sensing image of the multi-view overlapping unmanned aerial vehicle by using prior information:

there is certain overlap degree between the adjacent image of unmanned aerial vehicle aerial photograph. Because the data volume of taking photo by plane is very big, directly carry out three-dimensional reconstruction, on the one hand can not obtain better reconstruction effect, on the other hand can make the calculated amount of rebuilding big, the reconstruction time is longer on the other hand. Therefore, the images are grouped by using the existing POS information and the air belt prior information. Since the course overlapping degree of the unmanned aerial vehicle images is 80% in implementation and the side overlapping degree is 35%, 4 continuous images of the same flight band and two continuous images between the flight bands are divided into a group for a certain image. As shown in fig. 2, the image of point a is grouped in the space band diagram, and the black rectangular dotted line frame part is the image grouped in the same group as the image a.

(2) And (3) performing aerial three-dimensional photogrammetry on the basis of the step (1), solving external orientation elements of each image by utilizing an aerial three-dimensional mesh, and performing integral adjustment of a beam method. The step can be implemented by the prior art, and the invention is not described in detail.

(3) And (3) according to image grouping, performing rapid dense matching by using a PMVS algorithm accelerated by a GPU in the prior art on the basis of the step (2), generating dense three-dimensional point cloud, and using the reconstructed point cloud as three-dimensional elevation data. And step two, splicing the remote sensing images of the unmanned aerial vehicle after the adjustment. The specific process is as follows:

(1) feature extraction: and (5) utilizing SIFT to extract the features of the image.

(2) Image registration: firstly, coarse registration is carried out, and matched feature points are searched by using a k-d tree; and then carrying out fine registration, wherein the rough registration often has wrong matching points, so that the RANSAC algorithm is used for removing the wrong matching points. By registration of the images, a transformation matrix between the images can be obtained.

(3) And (3) splicing images: and (3) splicing images through the transformation matrix obtained in the step (2).

(4) Fusing images: and after splicing, fusing the images by using a bilinear interpolation algorithm.

And step three, filtering the color point cloud. The specific implementation is as follows:

(1) ground points and non-ground points of the point cloud are separated using improved morphological filtering. The operation process is as follows: firstly, any point and neighborhood points thereof are taken to form a window with a fixed size, the lowest point in the window is detected through morphological opening operation, if the difference between the elevation value of the point in the window and the elevation of the lowest point is in a threshold range, the point is represented as a ground point, and all points in the point cloud are taken out for filtering; then according to y2 xw^k+1 (where k is the number of iterations, and w is the size of the previous filtering window) to obtain the size of the window required by the next filtering, and the size of the window is smaller than the preset maximum value of the filtering window, and morphological filtering is performed once; and finally, when the window is larger than the preset window, finishing filtering.

(2) And (3) filtering the vegetation by utilizing the color invariant in the ground points in the step (1). Because the point cloud generated by image dense matching has color information, the green vegetation can be filtered by using a color invariant theory. Let the coordinates of each point in the point cloud be (x, y, z), the color three channels be (R, G, B), and the threshold value of the color invariant for vegetation be T_gThe color invariant formula defined using the green and blue channels is:

<math> <mrow> <msub> <mi>&psi;</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>4</mn> <mi>&pi;</mi> </mfrac> <mo>&times;</mo> <mi>arctan</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein, I_g(x,y,z)、I_b(x, y, z) represents the green and blue component values of the point cloud at the (x, y, z) point. Psi_g(x, y, z) represents a color invariant value at the (x, y, z) point. When psi_g＜T_gWhen the point is a vegetation point, the point is represented as a vegetation point; when psi_g＞T_gWhen it is, it means that the point is a non-vegetation point.

(3) Based on the characteristics of the building, for example, the height should not be less than 2 meters, and the area should not be less than 35 square meters. Non-building target points are filtered out using a threshold.

Step four, utilizing a region growing method based on plane fitting to segment the point cloud in the step three to obtain a point cloud region V of each building_Bi. The specific process is as follows:

(1) points in the point cloud are divided into two categories: if one point in the neighborhood of the point is default, the point belongs to the boundary point; otherwise such points belong to interior points.

(2) Let a certain interior point p₀(x₀,y₀,z₀),p₀The eight neighborhood points of (a): { p₁,p_2,…,p₈}. And fitting the plane equation (2) to the 9 points by using a least square method, wherein the method is specifically realized as follows:

z＝ax+by+c (2)

wherein a, b and c represent parameters of a plane equation of formula (2); for point p₀And 8 points within the neighborhood:

(x_i,y_i,z_i) I is 0,1,2, …, 7; fitting the calculated equation (2) to the plane equation minimizes equation (3).

<math> <mrow> <mi>S</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>7</mn> </munderover> <msup> <mrow> <mo>(</mo> <mi>ax</mi> <mo>+</mo> <mi>by</mi> <mo>+</mo> <mi>c</mi> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>

To minimize S, one should satisfy:

from this the following normal system of equations can be derived:

<math> <mrow> <mfenced open='|' close='|'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>x</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>y</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>n</mi> </mtd> </mtr> </mtable> </mfenced> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>b</mi> </mtd> </mtr> <mtr> <mtd> <mi>c</mi> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mi>i</mi> <mo>=</mo> <mn>0,1</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mn>7</mn> <mo>;</mo> <mi>n</mi> <mo>=</mo> <mn>8</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>

(3) the SSD and the variance of any interior point in the 8 neighborhood are found according to equation (4):

<math> <mrow> <mi>SSD</mi> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>&Element;</mo> <mi>M</mi> </mrow> </munder> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein M is p₀And its set of eight neighborhood points, h_kAnd z_kRespectively, an observation elevation value and a fitting plane elevation value.

(4) And sequencing each point in the point cloud in an ascending order according to the SSD value, and taking the point with the minimum SSD value as a seed point.

(5) Calculating the distance h from each point to the plane of the seed point in the neighborhood of the seed point, if h<h_TAnd merging into the same area and marking the area. Otherwise, this is a non-patch point. The corresponding threshold value can be preset by the person skilled in the art.

(6) And (5) when all the neighborhood points are traversed, finding unprocessed points from the new SSD value as new seed points, and repeating the operation in the step (5) until all the points are traversed.

Step five, carrying out point cloud area V on each building in the step four_BiAnd deleting wall points of the building, and finally obtaining the rough contour information of the building by solving the boundary points of the top surface. The specific implementation process is as follows:

(1) according to the divided surface patches in the step four, if the normal vector of the surface patch is parallel to the ground, the description is a wall surface patch; if the normal vector of the patch is perpendicular to the ground, it is described as a roof patch. The wall surface patches are deleted according to the normal vectors.

(2) And (4) obtaining the boundary of the point cloud by utilizing an alpha-shape algorithm according to the top surface point cloud obtained in the step (1). The specific method comprises the following steps: a) extracting an arbitrary point p1 from the top surface point cloud, searching a point set p2 with a distance of less than or equal to 2 alpha (alpha is the radius of a circle) from the rest points, and setting p2 to { p21, p22, p23, … and p2n }; b) taking out a point p2i from any point p2, and finding out the circle center p0 of the points p1 and p2i by using the formulas (6) and (7); knowing two points (x1, y1), (x2, y2) and the radius α of the circle, the equation for the center of the circle (x0, y0) is as follows:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>1</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>1</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>2</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>2</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>

it is difficult to directly solve the equation, so the distance intersection algorithm in mapping is used to obtain:

$\left\{\begin{array}{c}x0=x1+\frac{x2-x1}{2}+H(y2-y1)\\ y0=y1+\frac{y2-y1}{2}+H(x2-x1)\end{array}\right.---\left(7\right)$

wherein,S²＝(x1-x2)²+(y1-y2)²。

c) finding the distances l from all points to p0 from the point set p2, if l > alpha, then p1 and p2i are boundary points; if l < alpha, the next step d) is carried out; d) repeating the three steps a) b) c) for other points in p2 until all points in p2 are judged to be finished.

Step six, the top surface point cloud contour points of each building obtained in the step five are superposed on the image to obtain a building buffer area omega on the image_Bi. Since the boundary points of the top surface of the building in the step five are obtained by image multi-view stereo matching reconstruction, the points are re-projected back to the image only by utilizing a two-dimensional to three-dimensional projection matrix for inverse operation. The specific process is as follows:

setting a projection matrix P_3×4Computing P as a known inside and outside orientation element representing a single view geometry_3×4Matrix:

$P=\left[\begin{array}{ccc}-f& 0& x_0\\ 0& f& y_0\\ 0& 0& 1\end{array}\right]R^T\left[\begin{array}{cccc}1& 0& 0& -X_s\\ 0& 1& 0& -Y_s\\ 0& 0& 1& -Z_s\end{array}\right]---\left(8\right)$

wherein P is the projection matrix P_3×4F is the focal length of the image, x₀And y₀The eccentricity of the optical axis from the optical center in the horizontal direction and the vertical direction. X_s、Y_s、Z_sAs the coordinates of the camera center in the world coordinate system, R^TRepresenting a 3 x 3 rotation matrix.

The projection formula is:

$\left[\begin{array}{c}\mathrm{xz}\\ \mathrm{yz}\\ z\end{array}\right]=P\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(9\right)$

the projection of the point on the point cloud onto the image is calculated using equations (8) and (9). X, Y and Z represent the coordinates of the object point in the camera coordinate system, and X, Y and Z represent the coordinates of the object point in the world coordinate system.

And step seven, performing boundary connection on the points which are superposed on the image in the step six, superposing the obtained building rough contour on the spliced image as a buffer area, and simultaneously, utilizing the shape of the building rough contour as prior information and evolving the building precise contour in the buffer area by using a level set algorithm. The specific process is as follows:

(1) and carrying out boundary tracking on the points superposed on the image to obtain a closed contour. Expanding the contour to obtain a buffer region omega for extracting the contour of the building on the image_Bi。

(2) And (3) obtaining the contour of the image area as an initial level set according to the step (1), and performing local level set evolution in the remote sensing image in the buffer area by using the contour shape of the image area obtained according to the step (1) as prior information to obtain the fine contour of the building of the remote sensing image. The specific implementation is as follows:

first, the buffer area of each building on the image is preset to be omega_BiThe shape of the contour obtained in (1) is [ ]_Pi，Ф_PiAs a prior shape, phi_SiAs a target segmentation contour, an energy functional based on shape constraint is defined as follows:

E_to＝E(c₁,c₂,φ_Si)+E_sh(φ_Si,φ_P1) (10)

the initial contour and the prior shape of the target can be expressed by a moving target region:

<math> <mrow> <msub> <mi>E</mi> <mi>sh</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mo>,</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mo>&Integral;</mo> <msub> <mi>&Omega;</mi> <mi>Bi</mi> </msub> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>dxdy</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow> </math>

the energy function of the prior art based on the level set C-V model is:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <mi>E</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>,</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <mi>&mu;</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> <mi>v</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mi>dxdy</mi> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow> </math>

where H (φ) is a Heaviside function of the form: <math> <mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>1</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&GreaterEqual;</mo> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&lt;</mo> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> </math>

(phi) is a Dirac function of the formu₀(x, y) is the gray value of a certain point in the image area to be processed,modulus of the gradient of the current point, coefficient lambda₁，λ₂More than 0, mu, v ≧ 0 are fixed parameters, and it is generally recommended to take lambda₁＝λ₂1, μ ═ 0, v ═ 1, and (12) partial differential equations corresponding to the following equations:

<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>[</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> </mrow> </mfrac> <mo>-</mo> <mi>v</mi> <mo>-</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <msup> <mrow> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>]</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>]</mo> </mrow> <mn>2</mn> </msup> <mo>]</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein the parameter c₁、c₂Obtained according to the following formula,

<math> <mfenced open='' close=''> <mtable> <mtr> <mtd> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> </mfrac> </mtd> <mtd> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> </math>

according to the formula (13), a C-V level set evolution method is adopted to segment and extract the building region which is used as an image

Object v of building segmentation_iSegThe C-V level set is prior art and is not described in detail herein.

Is a distance d_Buffer1V is the corresponding position on the remote sensing image after the disaster_iEstablishing a local buffer v_iBufBuffer v_iBufThe image area within the outline is the current area omega to be processed. d_Buffer1Can be preset by the technicians in the field according to the situation.

Claims (10)

1. A building contour line extraction method of unmanned aerial vehicle multi-overlapping remote sensing images is characterized by comprising the following steps:

the method comprises the following steps that firstly, adjustment is carried out on three pairs of unmanned aerial vehicle remote sensing images by utilizing an aerial vehicle, and meanwhile, the images are subjected to dense matching by utilizing a PMVS algorithm after GPU acceleration, and finally dense color point clouds with high precision are obtained;

splicing the remote sensing images of the unmanned aerial vehicle after the leveling;

step three, filtering the color point cloud;

firstly, an improved morphological filtering algorithm is utilized to separate the ground from the non-ground, then, color invariants are utilized to filter vegetation in ground points, and finally, elevations and areas are utilized as thresholds to filter non-buildings;

detecting buildings in the point cloud by using a region growing algorithm;

deleting the wall surface of the building, and fitting the boundary of the top surface to finally obtain the rough contour information of the building;

step six, the coarse building contour obtained in the step three is used as a buffer area for building contour extraction and is superposed on the spliced image;

and step seven, simultaneously, utilizing the shape of the building rough contour as prior information, and evolving the building precise contour in the buffer area by using a level set algorithm.

2. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein the first step comprises the following steps:

(1.1) preprocessing the remote sensing image of the multi-view overlapping unmanned aerial vehicle by using prior information:

(1.2) carrying out aerial three-dimensional photogrammetry on the basis of the step (1.1), solving external orientation elements of each image by utilizing an aerial three-dimensional mesh, and carrying out integral adjustment of a beam method;

and (1.3) according to the image groups, on the basis of the step (1..2), performing rapid dense matching by using a PMVS (graphics processing unit) algorithm accelerated by a GPU (graphics processing unit) in the prior art to generate dense three-dimensional point cloud, and using the reconstructed point cloud as three-dimensional elevation data.

3. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein the second step comprises the following steps:

(2.1) feature extraction: extracting the features of the image by using SIFT;

(2.2) image registration: firstly, coarse registration is carried out, and matched feature points are searched by using a k-d tree; then, fine registration is carried out, and rough registration often has wrong matching points, so that the RANSAC algorithm is used for eliminating the wrong matching points; obtaining a transformation matrix between the images through the registration of the images;

(2.3) image stitching: splicing images through the transformation matrix obtained in the step (2.2);

(2.4) fusion of images: and after splicing, fusing the images by using a bilinear interpolation algorithm.

4. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein the third step comprises the following steps:

(3.1) separating the ground points and the non-ground points of the point cloud by using improved morphological filtering;

(3.2) filtering the vegetation by using the color invariant in the ground points in the step (3.1);

and (3.3) filtering out non-building target points by using a threshold value based on the characteristics of the building.

5. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapped remote sensing image according to claim 4, wherein the step (3.1) comprises the following processes:

firstly, any point and neighborhood points thereof are taken to form a window with a fixed size, the lowest point in the window is detected through morphological opening operation, if the difference between the elevation value of the point in the window and the elevation of the lowest point is in a threshold range, the point is represented as a ground point, and all points in the point cloud are taken out for filtering;

then according to y2 xw^k+1, obtaining the size of a window required by next filtering, wherein the size of the window is smaller than the maximum value of a preset filtering window, and performing morphological filtering again; finally, when the window is larger than the preset window, finishing filtering; wherein k is the iteration number, and w is the size of the previous filtering window.

6. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapped remote sensing image according to claim 4, wherein the step (3.2) comprises the following processes:

because the point cloud generated by image dense matching has color information, the green vegetation is filtered by utilizing the color invariant theory; let the coordinates of each point in the point cloud be (x, y, z), the color three channels be (R, G, B), and the threshold value of the color invariant for vegetation be T_gThe color invariant formula defined using the green and blue channels is:

<math> <mrow> <msub> <mi>&psi;</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>4</mn> <mi>&pi;</mi> </mfrac> <mo>&times;</mo> <mi>arctan</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>I</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>b</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein, I_g(x,y,z)、I_b(x, y, z) represents the green and blue component values of the point cloud at the (x, y, z) point; psi_g(x, y, z) represents a color invariant value at the (x, y, z) point; when psi_g＜T_gWhen the point is a vegetation point, the point is represented as a vegetation point; when psi_g＞T_gWhen it is, it means that the point is a non-vegetation point.

7. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image as claimed in claim 1, wherein in the fourth step, the point cloud in the third step is segmented by using a region growing method based on plane fitting to obtain the point cloud region V of each building_BiThe method specifically comprises the following steps:

(4.1) the points in the point cloud are divided into two types: if one point in the neighborhood of the point is default, the point belongs to the boundary point; otherwise such points belong to interior points;

(4.2) setting a certain interior point p₀(x₀,y₀,z₀),p₀The eight neighborhood points of (a): { p₁,p₂,…,p₈And fitting a plane equation (2) to the 9 points by using a least square method, and specifically realizing the following steps:

z＝ax+by+c (2)

for point p₀And 8 points within the neighborhood: (x)_i,y_i,z_i) I is 0,1,2, …, 7; fitting calculation formula (2) mean

A surface equation that minimizes equation (3);

<math> <mrow> <mi>S</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>7</mn> </munderover> <msup> <mrow> <mo>(</mo> <mi>ax</mi> <mo>+</mo> <mi>by</mi> <mo>+</mo> <mi>c</mi> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>

to minimize S, one should satisfy:

from this the following normal system of equations can be derived:

<math> <mrow> <mfenced open='|' close='|'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>x</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msup> <msub> <mi>y</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> </mtd> <mtd> <mi>n</mi> </mtd> </mtr> </mtable> </mfenced> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>b</mi> </mtd> </mtr> <mtr> <mtd> <mi>c</mi> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>x</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>y</mi> <mi>i</mi> </msub> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mi>&Sigma;</mi> <msub> <mi>z</mi> <mi>i</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>0,1</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mn>7</mn> <mo>;</mo> <mi>n</mi> <mo>=</mo> <mn>8</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>

(4.3) finding the SSD and the variance of any interior point in the 8 neighborhood according to the formula (4):

<math> <mrow> <mi>SSD</mi> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>&Element;</mo> <mi>M</mi> </mrow> </munder> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein M is p₀And its set of eight neighborhood points, h_kAnd z_kRespectively obtaining an observation elevation value and a fitting plane elevation value;

(4.4) sequencing each point in the point cloud in an ascending order according to the SSD value, and taking the point with the smallest SSD value as a seed point;

(4.5) calculating the distance h from each point to the plane where the seed points are located in the neighborhood of the seed points, if h is<h_T,h_TIf the elevation threshold is the elevation threshold, merging the elevation threshold into the same area, and marking the area; otherwise, this is a non-patch point;

and (4.6) when all the neighborhood points are traversed, finding unprocessed points from the new SSD value as new seed points, and repeating the operation of (4.5) until all the points are traversed.

8. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein in the fifth step, the point cloud area V of each building in the fourth step_BiDeleting wall points of the building, and finally obtaining rough contour information of the building by solving top surface boundary points, wherein the method specifically comprises the following steps:

(5.1) according to the divided surface patches in the step four, if the normal vector of the surface patch is parallel to the ground, the surface patch is indicated as a wall surface patch; if the normal vector of the patch is vertical to the ground, the patch is a roof patch; therefore, deleting the wall surface patches according to the normal vector;

(5.2) obtaining the boundary of the point cloud by utilizing an alpha-shape algorithm according to the top surface point cloud obtained in the step (5.1); the specific method comprises the following steps:

a) taking an arbitrary point p1 from the top surface point cloud, searching a point set p2 with a distance less than or equal to 2 alpha from the remaining points, wherein a is the radius of a circle, and setting p2 to { p21, p22, p23,.., p2n };

b) taking out a point p2i from any point p2, and finding out the circle center p0 of the points p1 and p2i by using the formulas (6) and (7); knowing two points (x1, y1), (x2, y2) and the radius α of the circle, the equation for the center of the circle (x0, y0) is as follows:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>1</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>1</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>2</mn> <mo>-</mo> <mi>x</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>2</mn> <mo>-</mo> <mi>y</mi> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>

it is difficult to directly solve the equation, so the distance intersection algorithm in mapping is used to obtain:

$\left\{\begin{array}{c}x0=x1+\frac{x2-x1}{2}+H(y2-y1)\\ y0=y1+\frac{y2-y1}{2}+H(x2-x1)\end{array}\right.---\left(7\right)$

wherein, <math> <mrow> <mi>H</mi> <mo>=</mo> <msqrt> <mfrac> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <msup> <mi>s</mi> <mn>2</mn> </msup> </mfrac> <mo>-</mo> <mfrac> <mn>1</mn> <mn>4</mn> </mfrac> </msqrt> <mo>,</mo> <msup> <mi>S</mi> <mn>2</mn> </msup> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mn>1</mn> <mo>-</mo> <mi>x</mi> <mn>2</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>y</mi> <mn>1</mn> <mo>-</mo> <mi>y</mi> <mn>2</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>;</mo> </mrow> </math>

c) finding the distances l from all points to p0 from the point set p2, if l > alpha, then p1 and p2i are boundary points; if l < alpha, the next step d) is carried out;

d) repeating the three steps a) b) c) for other points in p2 until all points in p2 are judged to be finished.

9. The method for extracting the contour line of the building of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein in the sixth step, the top point cloud contour points of each building obtained in the fifth step are superposed on the image to obtain a building buffer omega on the image_Bi(ii) a The specific process is as follows:

setting a projection matrix P_3×4Representing single-view geometryKnowing the inside and outside orientation elements, calculate P_3×4Matrix:

$P=\left[\begin{array}{ccc}-f& 0& x_0\\ 0& f& y_0\\ 0& 0& 1\end{array}\right]R^T\left[\begin{array}{cccc}1& 0& 0& -X_s\\ 0& 1& 0& -Y_s\\ 0& 0& 1& -Z_s\end{array}\right]---\left(8\right)$

wherein P is the projection matrix P_3×4F is the focal length of the image, x₀And y₀The eccentricity of the optical axis from the optical center in the horizontal direction and the vertical direction; x_s、Y_s、Z_sAs the coordinates of the camera center in the world coordinate system, R^TRepresents a 3 × 3 rotation matrix;

the projection formula is:

$\left[\begin{array}{c}\mathrm{xz}\\ \mathrm{yz}\\ z\end{array}\right]=P\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]---\left(9\right)$

calculating the projection of the point on the point cloud on the image by using the formula (8) and the formula (9); x, Y and Z represent the coordinates of the object point in the camera coordinate system, and X, Y and Z represent the coordinates of the object point in the world coordinate system.

10. The method for extracting the building contour line of the unmanned aerial vehicle multi-overlapping remote sensing image according to claim 1, wherein in the seventh step, the points overlapped on the image in the sixth step are subjected to boundary connection, the obtained building rough contour is overlapped on the spliced image as a buffer area, and meanwhile, the accurate contour of the building is evolved in the buffer area by using a level set algorithm by using the shape of the building rough contour as prior information; the method comprises the following specific steps:

(7.1) carrying out boundary tracking on the points superposed on the image to obtain a closed contour; expanding the contour to obtain a buffer region omega for extracting the contour of the building on the image_Bi；

(7.2) performing local level set evolution in the remote sensing image in the buffer area by taking the contour of the image area obtained according to the step (7.1) as an initial level set and taking the contour shape of the image area obtained according to the step (7.1) as prior information to obtain a fine contour of a building of the remote sensing image; the specific implementation is as follows:

first, the buffer area of each building on the image is preset to be omega_BiThe shape of the contour obtained in (7.1) is [ phi ]_Pi，Ф_PiAs a prior shape, phi_SiAs a target segmentation contour, an energy functional based on shape constraint is defined as follows:

E_to＝E(c₁,c₂,φ_Si)+E_sh(φ_Si,φ_P1)(10)

the initial contour and the prior shape of the target are expressed by utilizing a moving target area:

<math> <mrow> <msub> <mi>E</mi> <mi>sh</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mo>,</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mo>&Integral;</mo> <msub> <mi>&Omega;</mi> <mi>Bi</mi> </msub> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>&phi;</mi> <mi>si</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>&phi;</mi> <mi>Pi</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>dxdy</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow> </math>

the energy function of the prior art based on the level set C-V model is:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <mi>E</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>,</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>u</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> </mtd> </mtr> <mtr> <mtd> <mi>&mu;</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>+</mo> <mi>&nu;</mi> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mi>dxdy</mi> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow> </math>

where H (φ) is a Heaviside function of the form: <math> <mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>1</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&GreaterEqual;</mo> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> <mo>,</mo> </mtd> <mtd> <mi>&phi;</mi> <mo>&lt;</mo> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> </math>

(phi) is a Dirac function of the formu₀(x, y) is the gray value of a certain point in the image area to be processed,modulus of the gradient of the current point, coefficient lambda₁，λ₂More than 0, mu, v is more than or equal to 0 as a fixed parameter, and generally, lambda is taken₁＝λ_２1, μ ═ 0, ν ═ 1, and (12) partial differential equations corresponding to the following equations:

<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>[</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> </mrow> </mfrac> <mo>-</mo> <mi>&nu;</mi> <mo>-</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <msup> <mo>]</mo> <mn>2</mn> </msup> <mo>+</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mo>[</mo> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mo>]</mo> <mn>2</mn> </msup> <mo>]</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein the parameter c₁、c₂Obtained according to the following formula,

<math> <mrow> <msub> <mi>c</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mi>dxdy</mi> </mrow> </mfrac> </mrow> </math> <math> <mrow> <msub> <mi>c</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msub> <mi>u</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> <mrow> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mo>[</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>]</mo> <mi>dxdy</mi> </mrow> </mfrac> </mrow> </math>

according to the formula (13), a C-V level set evolution method is adopted to segment and extract the building region, and the image is taken as a building segmentation object V_iSegThe C-V level set is prior art.