CN107679441B - Method for extracting height of urban building based on multi-temporal remote sensing image shadow - Google Patents
Method for extracting height of urban building based on multi-temporal remote sensing image shadow Download PDFInfo
- Publication number
- CN107679441B CN107679441B CN201710077106.5A CN201710077106A CN107679441B CN 107679441 B CN107679441 B CN 107679441B CN 201710077106 A CN201710077106 A CN 201710077106A CN 107679441 B CN107679441 B CN 107679441B
- Authority
- CN
- China
- Prior art keywords
- shadow
- building
- image
- height
- length
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/176—Urban or other man-made structures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
- G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Multimedia (AREA)
- Remote Sensing (AREA)
- Computer Graphics (AREA)
- Image Processing (AREA)
Abstract
The invention discloses a method for extracting urban building height information based on multi-temporal remote sensing image shadows, which aims to form a systematic, quick and accurate method, wherein two remote sensing images with different time phases (the two images have different solar altitude angles and satellite azimuth angles) are used for carrying out batch analysis and calculation on target building shadows, so that the height of a target building in an extraction area can be quickly and accurately obtained, and the efficiency of work such as measurement, three-dimensional modeling and the like can be improved; the specific implementation is mainly that the height value of a target building in a region is calculated by calculating the shadow length of the building and selecting a corresponding mathematical model according to the included angle relationship between a mark building on an image and the shadow of the mark building; two groups of building height results obtained by calculating two remote sensing images are applied, a fitting equation is constructed based on a least square method, and the height of the corrected target building is obtained, so that the precision of the calculation result is improved.
Description
Technical Field
The invention relates to the technical field of remote sensing image processing, mathematical model design and three-dimensional construction, in particular to a method for extracting height information of urban buildings based on multi-temporal remote sensing image shadows.
Background
Along with the development of times, the urbanization expansion speed of China is faster and faster, and the urban expansion comprises two aspects of urban horizontal space expansion and urban vertical space expansion: horizontal spatial expansion manifests itself as a sharp increase in the area of the built-up zone; the vertical space mainly represents the increase of the height of the urban building, and the expansion of the longitudinal height of the city relieves the pressure caused by the expansion of the horizontal space on one hand, and can form a compact urban space form on the other hand, thereby facilitating the expansion of decisions such as urban planning, resource allocation and the like. The building is an important part of the urban composition, and the height of the building is increased to reflect the expansion of the vertical space of the city. The development trend of the urban horizontal space is combined with the expansion condition of the vertical space, and comprehensive analysis is performed, so that implementation of relevant urban decisions, simulation of urban landscape forms and construction of different urban models are facilitated.
The traditional building height extraction method is mainly divided into the following three types:
1) and the direct reading method is used for directly reading the number of windows of the building in the vertical direction on the remote sensing image to obtain the floor information and estimating the height of the building. The method is convenient and quick, but has the problems that the extraction precision is difficult to guarantee, and the method is easily influenced by differences of different displacement amplitudes of the target building on the image, different image scale and the like.
2) The projection method is characterized in that the projection difference of a target building on an aerial photo is measured and substituted into a formula:
h＝(δL/R)×H
h is the height of the target building, δ L is the projection difference of the building, R is the distance between the actual position point of the building and the main point of the image on the image, and H is the flight height of the image. The method has large limitation and is only suitable for buildings with high estimated actual height and position points far away from the image principal point.
3) The altitude difference calculation algorithm is characterized in that two points at the top and the bottom of a target building are determined on one picture by selecting a pair of stereopairs, then the corresponding image points with the same name are found on the other picture, the distance between the two points at the top and the bottom is measured to be used as the height h of the building, the result delta L obtained by subtracting the parallax of the image points with the same name is used as the parallax comparison, and the parallax comparison is substituted into a formula:
h＝[ΔL/(L+d)]×H
d is the average baseline length corresponding to the stereo pair, and H is the aerial height of the photo. The method has the problems of complex processing process, large calculated amount, difficult control of precision and poor effect.
At present, in practical application, methods such as vehicle-mounted radar and unmanned aerial vehicle scanning are commonly used to quickly acquire parameter information such as height of a target building, so as to construct a three-dimensional model. The methods have the advantages of high efficiency, high speed, convenience and quickness, but have the relative defect of being greatly influenced by the limitation of time and space ranges. In the remote sensing image, the shadow area of the building is rich in various information, and in contrast, the remote sensing image shadow is used for extracting the building height, so that the limitation of time and space can be completely eliminated. Firstly, the method can be used for estimating the building heights of the target area in different periods, and further reproducing the building models in the past period; in addition, under the condition of low required precision requirement, the construction of a building group model in a large-scale area can be quickly realized.
The building height extraction by using the remote sensing image shadow needs to be carried out through a reasonable remote sensing processing process of the system, required information is extracted from the building height, and particularly, the projection length of the building in the shadow direction is obtained to calculate the height of the corresponding building. The method has the main principle that the shadow length data and the corresponding mathematical model are combined to analyze and calculate the height of the building, the calculation results of images in different time phases are utilized to evaluate the height, the height value is calculated out by smaller errors, and the distribution condition of urban building groups is summarized, so that the method has important practical significance for the improvement of urban decision efficiency, the reasonable allocation of urban resources, the urban planning direction and the like.
Therefore, at present, a new method which is simple and convenient to operate, wide in application range, relatively accurate in calculation result and convenient to control in precision is urgently needed to be found for estimating the height of the building.
Disclosure of Invention
Aiming at the problem of extracting the height information of the building by using the shadow, the invention aims to form a systematic, quick and accurate method, and the shadow of the target building is analyzed and calculated in batch by using two remote sensing images with different time phases (the two images have different solar altitude angles and satellite azimuth angles), so that the height of the target building in an extraction area can be quickly and accurately obtained, and the efficiency of the work such as measurement, three-dimensional modeling and the like can be improved; the specific implementation is mainly that the height value of a target building in a region is calculated by calculating the shadow length of the building and selecting a corresponding mathematical model according to the included angle relationship between a mark building on an image and the shadow of the mark building; two groups of building height results obtained by calculating two remote sensing images are applied, a fitting equation is constructed based on a least square method, and the height of the corrected target building is obtained, so that the precision of the calculation result is improved.
The technical scheme mainly comprises the following steps:
the method comprises the following steps: acquiring an original high-resolution remote sensing image containing a target area, and carrying out preprocessing processes such as cutting, radiometric calibration, atmospheric correction, waveband fusion and the like on the image to obtain target area images of 2015 and 2016;
step two: selecting a proper training area sample by combining an image classification method and empirical knowledge, performing supervision and classification processing on a target image by using a support vector machine classifier, performing image classification post-processing (principal component analysis and cluster analysis) on a classification result, extracting shadows and buildings in the classification post-processing result, converting the shadows and the buildings into vector areas, and loading the shadow vectors to the target area image for display;
step three: selecting a result with a better classified post-processing effect, drawing equidistant parallel line segments covering a target area along a shadow projection direction (satellite azimuth of original image data), and performing intersection processing on the parallel line segments in the shadow direction and the shadow area to obtain an intersection line segment in the shadow direction; or the shadow vector area is divided by using equidistant parallel line segments, so that the shadow area corresponding to each building can be divided into a plurality of irregular graphs with equal width;
step four: opening the intersection line segments of the shadow directions obtained in the third step, obtaining the length of the intersection line segment corresponding to each shadow area one by one, and obtaining an average value after removing the minimum value and the maximum value to be used as the projection direction length of the shadow on the image; or opening an attribute table corresponding to each shadow region, counting the area of the small region of the divided irregular polygon, solving the average value U of the polygons belonging to the same shadow region, then taking the irregular polygon as a parallelogram along the shadow projection direction and the outline of the building at the other edge, and dividing the area U by the distance to obtain the projection length (shadow length) of the building along the shadow direction;
step five: three different relations exist among the sun, the satellite and the building shadow, and the relation among the sun, the satellite and the building shadow is determined by selecting regular mark buildings with vertexes, extracting corresponding shadow vector areas of the regular mark buildings and loading the regular mark buildings to an image for analysis. Wherein, the included angle between the shadow central axis and the mark building central axis is theta, if theta is more than 0 degree and less than 90 degrees, the sun and the satellite are positioned at the same side, and the length of the shadow on the image obtained in the fourth step is less than the length of the actual shadow; if theta is larger than or equal to 90 degrees and smaller than 180 degrees, the sun and the satellite are positioned on different sides, but the length of the shadow on the image obtained in the step four is larger than the length of the actual shadow; if theta is 180 degrees, the sun and the satellite are positioned on different sides, but the length of the shadow obtained on the image in the step four is the length of the actual shadow;
step six: determining the relation among the Sun, the satellite and the shadow of the building by combining the condition of the included angle theta in the step five, establishing a corresponding mathematical model by utilizing the Sun Elevation angle (Sun Elevation) and satellite Elevation angle (SatelliteElation) information in the original image, and calculating the height of the building;
step seven: constructing a linear regression equation (namely linear fitting) by using the heights of the buildings in the two years obtained in the step six, so as to obtain a result with more accurate height of the target building;
step eight: and constructing a corresponding building three-dimensional model by combining the building height information in the sixth step and the ArcScene10.2 function by using the building surface vectors obtained by image classification post-processing.
The original high resolution remote sensing image used in the first step is the data of a WorldView-2 image purchased. With four multispectral bands with a resolution of 1.8 meters and one panchromatic band with a resolution of 0.5 meters. Finding a target area on the image, and cutting each wave band by utilizing a vector diagram of the target area; then, respectively carrying out radiometric calibration and atmospheric correction on the cut multispectral wave band images and panchromatic wave band images, and eliminating the influence of atmospheric molecules and aerosol scattering on the image quality; and finally, performing band fusion on the processed multispectral band and the panchromatic band to obtain a target area synthetic image with the resolution of 0.5 meter, so as to facilitate the next operation.
In the second step, the selection of the training area sample firstly needs to classify all the ground features in the target area image by virtue of the experience of manual interpretation, and the ground feature classes are classified into six classes of water bodies, buildings, shadows, roads, bare lands and green lands. In the ENVI, after about 20 samples are selected from each category, the sample separability is evaluated, and the samples with the separability larger than 0.9 meet the requirement and are suitable samples, otherwise, the samples need to be selected again. After the classification is supervised, the main component analysis and the cluster analysis are carried out on the classification result, so that the classified shadows and buildings are gathered together, the plaques are removed, and the calculation of the length of the shadows is facilitated.
The satellite azimuth in step three, and the solar altitude and satellite altitude in step six are all searched from the original image data file used. The number of the small irregular areas divided in the step three is as large as possible, so that each divided small area graph can be approximately drawn into a parallelogram.
In the fourth step, no matter the length of the shadow projection direction on the image is obtained by calculating the parallel intersection line segment or the divided shadow area, the effect is all the same, the specific implementation needs to be considered according to the distribution condition of the target building, and the mode of adopting the parallel intersection line segment is simpler and more convenient when the distribution is dense; when the distribution is sparse, the calculation is carried out by utilizing the divided shadow area, and the effect is better. The embodiment of the invention calculates the projection length in the shadow direction by using the segmented shadow area.
In the fifth step, three different relations exist among the sun, the satellite and the building shadow, and different mathematical models are selected according to the condition of the used image data. The method mainly comprises a case that the sun and the satellite are on the same side and two cases that the sun and the satellite are on different sides. The specific distinction needs to find a mark building on the used image data to ensure that the building and the shadow thereof have stronger symmetry and the shape similar to a marker post is the best, namely, the mark building needs to have a vertex and regular shape, the symmetry is good, meanwhile, the shadow quality also needs to be kept good, the mark building can be distinguished from nearby ground objects obviously, and the whole shape is not influenced by other surrounding factors. The selection of the mark building directly determines the selection of a mathematical calculation model, and has great influence on the accuracy and precision of results. Fig. 8 shows an image of the present invention including a target area and its landmark buildings, wherein the selected landmark building is a landmark of the target area implemented by the present invention: the large corn meets the requirements of the above-mentioned sign buildings. The vertex of the mark building is set to be A, the top end of the corresponding shadow is set to be A ', and the selected mark building has good symmetry and regular shadow, so that the mark building can pass through the top end A' of the shadow to be used as a central axis of a shadow area, namely a symmetric line, and is intersected with the boundary line of the building and the shadow at a point O. Then OA' is the central axis of the shadow, the connection point O is the apex A of the landmark building, and OA is the central axis of the building. The included angle between the shadow central axis and the mark building central axis is theta, if the included angle is 0 degrees < theta <90 degrees, the sun and the satellite are positioned at the same side, and the shadow length on the image obtained in the fourth step is smaller than the actual shadow length; if theta is larger than or equal to 90 degrees and smaller than 180 degrees, the sun and the satellite are positioned at different sides, and the length of the shadow on the obtained image is larger than the length of the actual shadow; if θ is 180 °, the sun and the satellite are located on different sides, and the projection length in the shadow direction is the actual shadow length.
The linear regression line in the step seven adopts the principle of least square method, and the calculation result of the building height in 2015 is taken as an independent variable xi2016 year result as dependent variable yiThen the regression equation can be set to
Find the corresponding average value
When corresponding to the residual error
When the minimum value is the linear regression line, the linear regression line is respectively substituted into
The coefficients a and b are obtained to obtain corresponding linear regression equation, and the linear regression equation is substituted into xiAnd obtaining a correction result of the building height estimated value. Further, since building No. 5 has irregularities, the height of building No. 5 in the correction result is its maximum height at this time. When the model is constructed, the model is divided into an upper part and a lower part, and different height values are respectively given to the planar vectors of the two parts of buildings, wherein the height value of the upper part of the building is the building height correction result of the building No. 5 in the table of FIG. 11, and the height of the lower part of the building is a height correction value obtained by calculation by utilizing the length of the projection of the corresponding part: 28.42m, based on which a corresponding model of building No. 5 can be constructed.
After adopting such design, the invention has at least the following advantages:
1) the images of two different time phases are selected for calculation, comparison analysis can be carried out, a regression linear equation is further constructed, the result is corrected, and the influence of errors generated in the process of extracting the height information is reduced. The actual height of the buildings No. 1 to 8 corresponding to the research area used by the invention is obtained through application survey and actual measurement of starfish total station ATS-320M, and is shown in the following table: unit (rice)
Wherein the building height correction value is the value of the corresponding point of the regression line made by the building heights of 15 years and 16 years. Therefore, the value of the building height for 15 years, the value for 16 years and the deviation square sum of the corrected value and the corresponding actual value can be respectively calculated, and the calculation formula is as follows:
wherein, aiIs the actual height of building No. i, yiFor the building height of i number corresponding to each situation, the corresponding building height can be obtained
Therefore, the height correction value calculated by using the linear regression equation obtained by the height values in two years is more accurate and has higher precision.
2) The relation among the sun, the satellite and the shadow of the building is considered thoroughly, different mathematical models are constructed, corresponding models can be adopted for calculation according to different selected image data, and the accuracy and precision of the result are improved.
3) The process of extracting the building height information is more systematic, the difficulty of work such as three-dimensional modeling and measurement is greatly reduced, the limitations of time and space factors are ignored, building models in different periods and different areas can be constructed conveniently and quickly.
4) The shadow length is calculated by using the intersection line segment in the shadow direction or the divided shadow area, so that the error generated in the process of obtaining the shadow length can be reduced, the precision of the calculation result is higher, and the accuracy of the estimated building height is ensured.
Drawings
The foregoing is an overview of the technical solutions of the present invention, and the present invention is described in detail below with reference to the accompanying drawings and the detailed description.
Fig. 1 is an image of the pre-treatment results in 2015 of the study area according to the present invention.
FIG. 2 is an image of 2016. pretreatment results in the study area of the present invention.
Fig. 3 is an 2015 image map of the study area loaded with shadow area vectors according to the present invention.
Fig. 4 is a 2016-year image map of a study area loaded with shadow plane vectors.
Figure 5 is an 2015 corresponding image of a study area of the present invention loaded with parallel line segmented shadow vectors.
FIG. 6 is a 2016 corresponding image of a study area loaded with parallel line segmented shadow vectors according to the present invention.
Fig. 7 is a table of shadow length data of buildings corresponding to years 2015 and 2016 in the study area of the present invention.
FIG. 8 is an image of a target area and its landmark building according to the present invention.
FIG. 9 is a diagram showing the correspondence between the landmark buildings and their shadows in the image used in the present invention.
FIG. 10 is a diagram of a mathematical model for calculating building height based on shadow length according to the present invention.
Fig. 11 is a 2015 and 2016 building height data table calculated based on shadow length in accordance with the present invention.
Fig. 12 is a table of the building height correction value data based on the two-year height data according to the present invention.
Fig. 13 is a three-dimensional model of a target area building of the present invention.
Detailed Description
The following detailed description of specific embodiments of the invention is provided in connection with the accompanying drawings.
The research area of the invention selects a part of Zhengzhou city gold water area, the data is WorldView-2 image data of 2015 and 2016, and the processing of each step is carried out by combining ENVI5.1 and Arcgis10.2 software. Referring to fig. 1 to 11, a method for rapidly constructing a three-dimensional model of a building based on multi-temporal remote sensing image shadows includes the following steps:
the method comprises the following steps: and acquiring an original high-resolution remote sensing image containing a target area, and carrying out preprocessing processes such as cutting, radiometric calibration, atmospheric correction, waveband fusion and the like on the image to obtain target area images of 2015 and 2016.
Step two: by combining an image Classification method and experience knowledge, after selecting a proper training area sample, selecting a Support Vector Machine classifier (Support Vector Machine Classification) to supervise and classify the sample, and then performing image Classification post-processing on a Classification result, wherein the Classification post-processing comprises major/minor component Analysis (major/minor Analysis) and clustering processing (column Classes), so that a classified shadow area is more uniform and full, the shadow area is convenient to identify, and the next step is facilitated. The result with better processing effect after classification is selected for the next processing, as shown in fig. 3 and 4,
step three: a linear vector file is newly built in the ArcCatalog, the linear vector file is loaded into the Arcmap to be used as a template for drawing a line, a line segment with a certain length along the projection direction (satellite azimuth angle) is built by utilizing an element construction tool or a COGO tool, a Parallel Copy (Copy Parallel) option is selected in a pull-down menu of an Editor (Editor), the distance (translation distance) of the Parallel Copy option is set to be 5 m, a first Parallel line of the line segment is built, all the line segments are sequentially selected to be translated until all shadow parts of a target area are intersected with the Parallel line segment. Then selecting Analysis Tools in a System tool box (System Tools), and acquiring line segments of intersection parts of the parallel line segments and the shadow areas by utilizing an intersection function of an extraction Analysis tool; then, by using a "dividing plane" tool with an "advanced editing" function in archap, a shadow vector area is divided by using an intersection line segment, or a shadow is directly divided by using an equidistant parallel line, so that a shadow area corresponding to each building can be divided into a plurality of irregular graphs with equal width, as shown in fig. 5 and 6.
Step four: opening the attribute table corresponding to each shadow region obtained in the third step, counting the areas of the small regions of the divided irregular polygons one by one, summarizing and summarizing the polygons belonging to the same shadow region, sequencing the polygons from large to small, taking the maximum value as a starting point, counting and summarizing the polygons with errors within 50 square meters, calculating the average value U of the polygons, then taking the irregular polygons as a parallelogram along the shadow projection direction and along the outline of the building, and dividing the average area U by the distance, wherein the corresponding height of the side in the shadow projection direction is the distance of the drawn parallel lines (the distance in the experiment is 5 meters), and the obtained result is the projection length (the shadow length) of the building along the shadow direction, and is shown in a table in fig. 7.
Step five: three different relations exist among the sun, the satellite and the building shadow, and the relation among the sun, the satellite and the building shadow is determined by selecting regular mark buildings with vertexes, extracting corresponding shadow vector areas and loading the shadow vector areas to the image. As shown in fig. 8, the landmark building selected for the implementation of the technology of the present invention is the CBD landmark of zhengdong new district in zhengzhou city: big corn. Let A be the vertex of the landmark building, A 'be the shadow projection point of the building vertex, OA' be the central axis of the shadow, where O is the intersection point of the central axis of the shadow and the building bottom, OA is the central axis of the building on the image, and the included angle between the central axis of the shadow and the central axis of the landmark building is θ, as shown in FIG. 9. As can be seen from the above, if θ < 0 ° <90 °, the sun and the satellite are located on the same side, and the length of the shadow on the image obtained in step four is smaller than the actual shadow length; if theta is more than or equal to 90 degrees and less than 180 degrees, the sun and the satellite are positioned at different sides, and the length of the shadow on the obtained image is larger than the length of the actual shadow; if θ is 180 °, the sun and the satellite are still on different sides, and the length of the shadow on the obtained image is the actual shadow length.
Obviously, in the remote sensing image of the target area in the implementation of the technology of the present invention, the included angle θ between the central axis of the mark building and the central axis of the shadow thereof satisfies the following condition: theta is more than or equal to 90 degrees and less than 180 degrees, so that the condition that the sun and the satellite are positioned on different sides and the length of the shadow on the image is larger than the length of the actual shadow is determined.
Step six: and determining the relation among the Sun, the satellite and the building shadow by combining the condition of the included angle theta in the fifth step, determining a corresponding mathematical model by utilizing the Sun Elevation angle (Sun Elevation) and the satellite Elevation angle (SatelliteElation) information in the original image, and calculating to obtain the building height by combining the length of the shadow direction projection on the image obtained in the above step. Among them, as shown in fig. 10, it needs to be considered in three cases:
1) the sun altitude β and the satellite altitude angle α are in different directions, namely the sun and the satellite are on different sides, if the included angle theta between the central axis of the mark building and the shadow central axis meets the condition that theta is 180 degrees, the actual shadow length corresponding to the building is the shadow direction projection length L calculated in the step2Corresponding to a building height of
H＝L2tanβ
2) the solar altitude β and the satellite altitude angle alphaThe same direction, that is, the sun and the satellite are on different sides, if the included angle theta between the central axis of the mark building and the central axis of the shadow satisfies: theta is more than or equal to 90 degrees<180 degrees, the actual shadow length S of the building is smaller than the shadow length L obtained in the step four2The method comprises the following steps:
S＝L2-H/tanα
the corresponding building height is then:
H＝S×tanβ
finishing to obtain H ═ L2×[(tanα+tanβ)/(tanα×tanβ)]
3) the sun altitude beta and the satellite altitude α are in the same direction, namely the sun and the satellite are on the same side, and the included angle theta between the central axis of the sign building and the central axis of the shadow satisfies 0 DEG<θ<And 90 degrees, then the shadow direction projection L obtained in the step four is carried out2Only a partial projection of the shadow, the total shadow length is:
S＝L1+L2
as shown in case 3) of FIG. 10, L1Is the length of the partial shadow that is occluded due to satellite altitude. At this time, the building height is
H＝(L1+L2)tanβ
And the two formulae are combined and substituted into L2Height of building can be obtained
H＝L2×[(tanα-tanβ)/(tanα×tanβ)]
The image of the experimental area is analyzed, and obviously, the included angle theta between the central axis of the marked building and the central axis of the shadow meets the following requirements: theta is larger than or equal to 90 degrees and smaller than 180 degrees, calculation is carried out according to the model of the situation 2) in the figure 10, and the height H of the building in the target area is obtained and is shown in a table in figure 11.
Step seven: constructing a linear regression equation (namely linear fitting) by using the building heights of the two years obtained in the step six, taking the height calculation result of the year 2015 as an independent variable x, and taking the height calculation result of the year 2016 as a dependent variable y, and obtaining a linear fitting relationship between the height calculation result and the dependent variable y as follows: 1.0063 x-0.8272, R20.9852. And the more accurate result of the corresponding building height is obtained, as shown in the table of fig. 12.
Step eight: and adding the height information of the building corrected in the seventh step into a corresponding building planar vector attribute table, opening a planar vector corresponding to the target building in ArcScene10.2, and stretching and displaying according to the height attribute to obtain a corresponding building three-dimensional model, as shown in FIG. 13.
The present invention has been described in further detail with reference to the embodiments, and it is not to be construed that the embodiments are limited thereto. For those skilled in the art to which the present invention pertains and related technologies, the extension, operation method and data replacement should fall within the protection scope of the present invention based on the technical solution of the present invention.
Claims (7)
1. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadow is characterized by comprising the following steps of:
the method comprises the following steps: obtaining an original high-resolution remote sensing image containing a target area, and carrying out cutting, radiometric calibration, atmospheric correction and waveband fusion pretreatment on the original high-resolution remote sensing image to obtain a target area image which is continuous for two years;
step two: selecting a proper training area sample by combining an image classification method and experience knowledge, performing supervision and classification processing on a target image by using a support vector machine classifier, performing image classification post-processing on a classification result, extracting shadows and buildings in the classification post-processing result, converting the shadows and the buildings into vector areas, and loading the shadow vectors to the target area image for display;
step three: selecting a result with accurate processing effect after classification, drawing equidistant parallel line segments covering a target area along the shadow projection direction, namely the satellite azimuth angle of the original image data, and performing intersection processing on the parallel line segments in the shadow direction and the shadow area to obtain an intersection line segment in the shadow direction; or the shadow vector area is divided by using equidistant parallel line segments, and the shadow area corresponding to each building is divided into a plurality of irregular graphs with equal width;
step four: extracting intersection line segments in the shadow direction obtained in the third step, obtaining the length of the intersection line segment corresponding to each shadow area one by one, and obtaining an average value after removing the minimum value and the maximum value to be used as the projection direction length of the shadow on the image; or opening an attribute table corresponding to each shadow region, counting the area of the small region of the divided irregular polygon, solving the average value U of the polygons belonging to the same shadow region, then regarding the irregular polygon as a parallelogram along the shadow projection direction and along the outline of the building, and dividing the area average value U by the distance to obtain the projection length of the building along the shadow direction;
step five: three different relations exist among the sun, the satellite and the building shadow, regular mark buildings with vertexes are selected, corresponding shadow vector areas are extracted, and the shadow vector areas are loaded on an image for analysis to determine the relation among the three; wherein, the included angle between the shadow central axis and the mark building central axis is theta, if theta is more than or equal to 0 degree and less than 90 degrees, the sun and the satellite are positioned at the same side, and the length of the shadow direction on the image obtained in the fourth step is less than the length of the actual shadow; if theta is larger than or equal to 90 degrees and smaller than 180 degrees, the sun and the satellite are positioned at different sides, but the length of the shadow direction on the image obtained in the step four is larger than the length of the actual shadow; if theta is 180 degrees, the sun and the satellite are positioned on different sides, but the length of the shadow direction obtained on the image in the step four is the actual shadow length;
step six: determining the relation among the sun, the satellite and the building shadow by combining the condition of the included angle theta in the step five, establishing a corresponding mathematical model by utilizing the sun altitude angle and the satellite altitude angle information in the original image, and calculating the height of the building;
step seven: constructing a linear regression equation by using the heights of the buildings in the two years obtained in the step six, so as to obtain a result of accurate height of the target building;
step eight: and (5) building three-dimensional models are constructed through ArcScene10.2 functions by utilizing the building surface vectors obtained by image classification post-processing and combining the building height information in the step seven.
2. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadow according to claim 1, wherein the original high-resolution remote sensing image used in the first step is obtained WorldView image data, the original high-resolution remote sensing image has four multispectral wave bands with the resolution of 1.8 m and a panchromatic wave band with the resolution of 0.5 m, a target area is found on the image, and each wave band is cut by utilizing a vector diagram of the target area; then, respectively carrying out radiometric calibration and atmospheric correction on the cut multispectral wave band images and panchromatic wave band images, and eliminating the influence of atmospheric molecules and aerosol scattering on the image quality; and finally, performing band fusion on the processed multispectral band and the panchromatic band to obtain a target area synthetic image with the resolution of 0.5 meter, so as to facilitate the next operation.
3. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadow as claimed in claim 1, wherein in the second step, the selection of the training area sample firstly needs to classify all the land features in the target area image by virtue of the experience of manual interpretation, and the categories of the land features in the second step are divided into six categories of water, buildings, shadows, roads, bare land and green land; in ENVI, after 20 samples are selected in each category, the samples are evaluated for separability, the samples with separability larger than 0.9 meet the requirement and are suitable samples, otherwise, the samples need to be selected again; after the classification is supervised, the main component analysis and the cluster analysis are carried out on the classification result, so that classified shadows and buildings are gathered together, the plaques are removed, and the calculation of the shadow length is facilitated.
4. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadows according to claim 1, wherein the satellite azimuth in the third step, and the solar altitude and the satellite altitude in the sixth step are searched from the original image data file; the number of the small irregular areas divided in the step three is more, so that each divided small area graph can be approximately abstracted into a parallelogram.
5. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadow is characterized in that the length of the shadow projection direction on the image is obtained by using parallel intersection line segments or the divided shadow area, the length needs to be considered according to the distribution condition of a target building, and the parallel intersection line segments are adopted when the distribution is dense; when the distribution is sparse, the calculation is performed by using the divided shadow area.
6. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadows according to claim 1, wherein three different relations exist among the sun, the satellite and the building shadows in the fifth step, and different mathematical models are selected according to the condition of image data; three relations among the sun, the satellite and the shadow of the building are as follows: the sun and the satellite are on the same side, and the sun and the satellite are on different sides; the specific distinction firstly needs to find a mark building on the used image data, ensures that the building and the shadow thereof have symmetry, selects the mark building similar to the shape of a marker post, namely the mark building needs to have a vertex and regular shape, has good symmetry, simultaneously has good shadow quality requirement, can obviously distinguish from nearby ground objects, and has no influence of other surrounding factors on the whole shape; the selection of the mark building directly determines the selection of a mathematical calculation model, and has influence on the accuracy and precision of results; setting the top point of the mark building as A and the top end of the corresponding shadow as A ', wherein the selected mark building has good symmetry and regular shadow, and the axis passing through the top end A' of the shadow is taken as the central axis of the shadow area, namely a symmetric line, and is intersected with the boundary line of the building and the shadow at a point O; OA' is the central axis of the shadow, the connection point O is the vertex A of the mark building, and then OA is the central axis of the building; the included angle between the central axis of the shadow and the central axis of the mark building is theta, if theta is more than or equal to 0 degree and less than 90 degrees, the sun and the satellite are positioned at the same side, and the shadow length on the image obtained in the fourth step is smaller than the actual shadow length; if theta is larger than or equal to 90 degrees and smaller than 180 degrees, the sun and the satellite are positioned at different sides, and the length of the shadow on the obtained image is larger than the length of the actual shadow; if θ is 180 °, the sun and the satellite are located on different sides, and the projection length in the shadow direction is the actual shadow length.
7. The method for extracting the height of the urban building based on the multi-temporal remote sensing image shadow as claimed in claim 1, wherein the linear regression equation in the seventh step adopts the principle of least square method, and the calculation result of the building height in the first year is used as an independent variable xiThe result of the second year is the dependent variable yiThen set the regression equation to
Find the corresponding average value
When corresponding to the residual error
When the minimum value is the linear regression line, the coefficients a and b are respectively obtained and substituted into xiAnd obtaining a correction result of the building height estimated value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710077106.5A CN107679441B (en) | 2017-02-14 | 2017-02-14 | Method for extracting height of urban building based on multi-temporal remote sensing image shadow |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710077106.5A CN107679441B (en) | 2017-02-14 | 2017-02-14 | Method for extracting height of urban building based on multi-temporal remote sensing image shadow |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107679441A CN107679441A (en) | 2018-02-09 |
| CN107679441B true CN107679441B (en) | 2020-06-02 |
Family
ID=61134007
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710077106.5A Expired - Fee Related CN107679441B (en) | 2017-02-14 | 2017-02-14 | Method for extracting height of urban building based on multi-temporal remote sensing image shadow |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107679441B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111721267A (en) * | 2020-07-22 | 2020-09-29 | 河南大学 | Method for predicting building height by using satellite image |
Families Citing this family (19)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108765488B (en) * | 2018-03-29 | 2022-03-04 | 武汉大学 | Shadow-based high-resolution remote sensing image building height estimation method |
| CN108733077A (en) * | 2018-05-31 | 2018-11-02 | 合肥赛为智能有限公司 | A kind of agricultural unmanned machine operation Intelligentized regulating and controlling system |
| CN109190236B (en) * | 2018-08-28 | 2023-05-09 | 山东建筑大学 | Method for characterizing surface morphology features of machined workpiece |
| CN109389051A (en) * | 2018-09-20 | 2019-02-26 | 华南农业大学 | A kind of building remote sensing images recognition methods based on convolutional neural networks |
| CN110222586A (en) * | 2019-05-15 | 2019-09-10 | 清华大学 | A kind of calculating of depth of building and the method for building up of urban morphology parameter database |
| CN110736435B (en) * | 2019-08-29 | 2021-05-14 | 昆明理工大学 | Height measuring device and method based on solar geometric optics |
| CN110390715B (en) * | 2019-09-11 | 2023-10-13 | 桂林理工大学 | Method for simultaneously detecting shadows of building roof, building wall and ground |
| CN111104850B (en) * | 2019-10-30 | 2023-09-26 | 中国四维测绘技术有限公司 | Remote sensing image building automatic extraction method and system based on residual error network |
| CN111398958B (en) * | 2020-04-03 | 2022-06-14 | 兰州大学 | Method for determining correlation between ground settlement and building height of loess excavation area |
| CN111597496B (en) * | 2020-05-18 | 2022-06-24 | 中国科学院空天信息创新研究院 | Rational polynomial coefficient-based single satellite-borne remote sensing image target height calculation method |
| CN111666910B (en) * | 2020-06-12 | 2024-05-17 | 北京博能科技股份有限公司 | Airport clearance area obstacle detection method and device and electronic product |
| CN111721268B (en) * | 2020-07-22 | 2021-09-07 | 河南大学 | Method and device for accurately inverting building height |
| CN113203399B (en) * | 2021-04-16 | 2022-02-15 | 青岛地质工程勘察院(青岛地质勘查开发局) | Underground space resource quantity analysis method |
| CN113487634B (en) * | 2021-06-11 | 2023-06-30 | 中国联合网络通信集团有限公司 | Method and device for associating building height and area |
| CN113408419A (en) * | 2021-06-18 | 2021-09-17 | 亿景智联(北京)科技有限公司 | Method for measuring and calculating building height by utilizing multi-feature fusion |
| CN113886937B (en) * | 2021-12-08 | 2022-04-05 | 深圳小库科技有限公司 | Rapid projection calculation method based on irregular three-dimensional space object |
| CN114373009B (en) * | 2022-01-13 | 2023-02-03 | 中国科学院空天信息创新研究院 | Building shadow height measurement intelligent calculation method based on high-resolution remote sensing image |
| CN114061548B (en) * | 2022-01-14 | 2022-06-17 | 山东省地质测绘院 | Building surveying and mapping method and system based on unmanned aerial vehicle remote sensing |
| CN116977469B (en) * | 2023-08-02 | 2024-01-23 | 中国水利水电科学研究院 | Community scale city form data batch generation method based on random slicing |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102147936A (en) * | 2011-03-09 | 2011-08-10 | 浙江大学 | Cascade-based method for seamlessly superposing two-dimensional vectors on three-dimensional topography surface |
| CN104463868A (en) * | 2014-12-05 | 2015-03-25 | 北京师范大学 | Rapid building height obtaining method based on parameter-free high-resolution image |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2455359C (en) * | 2004-01-16 | 2013-01-08 | Geotango International Corp. | System, computer program and method for 3d object measurement, modeling and mapping from single imagery |
-
2017
- 2017-02-14 CN CN201710077106.5A patent/CN107679441B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102147936A (en) * | 2011-03-09 | 2011-08-10 | 浙江大学 | Cascade-based method for seamlessly superposing two-dimensional vectors on three-dimensional topography surface |
| CN104463868A (en) * | 2014-12-05 | 2015-03-25 | 北京师范大学 | Rapid building height obtaining method based on parameter-free high-resolution image |
Non-Patent Citations (3)
| Title |
|---|
| Satellite images analysis for shadow detection and building height estimation;Gregoris Liasis等;《ISPRS Journal of Photogrammetry and Remote Sensing》;20141128;第119卷;第437-540页 * |
| SHADOW ANALYSIS TECHNIQUE FOR EXTRACTION OF BUILDING HEIGHT USING HIGH RESOLUTION SATELLITE SINGLE IMAGE AND ACCURACY ASSESSMENT;P.L.N. Rajua等;《The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences》;20141128;第1185-1192页 * |
| 高分辨率遥感影像阴影与立体像对提取建筑物高度比较研究;陈亭等;《遥感科学与应用技术》;20160930;第18卷(第9期);第1267-1275页 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111721267A (en) * | 2020-07-22 | 2020-09-29 | 河南大学 | Method for predicting building height by using satellite image |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107679441A (en) | 2018-02-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107679441B (en) | Method for extracting height of urban building based on multi-temporal remote sensing image shadow | |
| CN110570428B (en) | Method and system for dividing building roof sheet from large-scale image dense matching point cloud | |
| CN110717983B (en) | Building elevation three-dimensional reconstruction method based on knapsack type three-dimensional laser point cloud data | |
| Chen et al. | Automatic extraction of blocks from 3D point clouds of fractured rock | |
| Matei et al. | Building segmentation for densely built urban regions using aerial lidar data | |
| CN102074047B (en) | High-fineness urban three-dimensional modeling method | |
| CN102147250B (en) | Digital line graph mapping method | |
| Wang et al. | Modeling indoor spaces using decomposition and reconstruction of structural elements | |
| CN109949326A (en) | Contour of building line drawing method based on Backpack type three-dimensional laser point cloud data | |
| CN112595258A (en) | Ground object contour extraction method based on ground laser point cloud | |
| CN110780307B (en) | Method for obtaining road cross section based on storage battery car-mounted laser point cloud mobile measurement system | |
| CN107844802A (en) | Water and soil conservation value method based on unmanned plane low-altitude remote sensing and object oriented classification | |
| CN103324916B (en) | Vehicle-mounted and aviation LiDAR data method for registering based on building profile | |
| CN101915570B (en) | Vanishing point based method for automatically extracting and classifying ground movement measurement image line segments | |
| CN114998536A (en) | Model generation method and device based on novel basic mapping and storage medium | |
| CN106338277B (en) | A kind of building change detecting method based on baseline | |
| Sun et al. | Building displacement measurement and analysis based on UAV images | |
| Liu et al. | Intelligent scanning for optimal rock discontinuity sets considering multiple parameters based on manifold learning combined with UAV photogrammetry | |
| CN106446306A (en) | Gauss sphere cluster based machine part reverse engineering modeling method | |
| Demir | Automated detection of 3D roof planes from Lidar data | |
| Shen et al. | A novel baseline-based method to detect local structural changes in masonry walls using dense terrestrial laser scanning point clouds | |
| Namouchi et al. | Piecewise horizontal 3d roof reconstruction from aerial lidar | |
| CN110298852A (en) | Geological boundary extraction method based on unmanned plane image chromatography | |
| CN115984721A (en) | Method for realizing country landscape management based on oblique photography and image recognition technology | |
| CN112365534B (en) | Large coal pile volume measurement method based on monocular camera three-dimensional reconstruction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200602 Termination date: 20210214 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |