CN110070545A - A kind of method that textural characteristics density in cities and towns automatically extracts cities and towns built-up areas - Google Patents

Abstract

The present invention claims to protect a method for automatically extracting urban built-up area from urban texture feature density. The invention belongs to the field of remote sensing image information extraction. The method uses multi-scale simple linear iterative algorithm to construct remote sensing image urban texture feature density for high-resolution remote sensing data. , to automatically extract urban built-up areas. The method first uses the J‑M distance to analyze the sensitivity of different texture features of the gray level co-occurrence matrix in town extraction, and uses the texture difference features as the classification basis; then uses the Simple Linear Iterative Clustering (SLIC) to perform multiple scale Using superpixels as the range for calculating the local feature density, the feature density images after multiple segmentations are calculated and superimposed, and the feature density images of the images to be classified are constructed. The results extracted by this method have better completeness and correctness, especially in the aspect of completeness, the town edge extracted by this method is highly consistent with the actual edge, which can achieve the effect of manual detection.

Description

A method for automatic extraction of urban built-up area by urban texture feature density

technical field

The invention belongs to a method in the field of remote sensing image information extraction, and specifically relates to a method for high-resolution remote sensing image (10m and higher resolution), based on multi-layer superpixel segmentation, constructing urban texture feature density, and realizing automatic extraction of image data. method of urban built-up area.

Background technique

Urban information extraction based on remote sensing images refers to identifying and marking urban built-up areas in remote sensing satellite images, aerial images or images captured by drones. The scope, distribution and changes of urban built-up areas are of great significance to urban planning, land use resource allocation and environmental monitoring. They are also the data basis for many research directions, and have attracted the attention of a large number of researchers. However, complex ground information and diverse urban settlements make town extraction a very challenging task. With the rapid development of the global aerospace industry, the quality and resolution of acquired image data have improved, and high-resolution images have become easier to obtain. High-resolution remote sensing images have clear urban texture information, which is based on remote sensing image information. Extract important data sources provided by town information.

The extraction methods of urban built-up areas are mainly divided into two categories, spectrum-based methods and texture-based methods. Spectral-based methods refer to analyzing the differences in spectral curves between urban and non-urban cities through multi-spectral information of low- and medium-resolution data, and are mostly used to extract large areas of towns and detect changes. However, low- and medium-resolution data are limited by resolution, making it difficult to obtain high-precision urban areas. Texture-based methods refer to distinguishing between urban and non-urban texture differences in high-resolution data. Texture-based methods can be further divided into: edge density-based methods, line detection-based methods, and texture detection model-based methods. The first two of these rely too much on the quality of the input data and require the data to have a resolution higher than 2 meters. The method based on texture detection model is widely used. Gabor filtering and gray co-occurrence matrix are usually used to construct texture features. Among them, Gabor filtering has higher requirements on the resolution of data, while the method based on gray co-occurrence matrix can be applied to a variety of resolutions. data. Existing texture-based methods have high requirements on data quality, usually requiring data with a resolution higher than 5 meters. However, high-resolution data acquisition is expensive, which is not conducive to the extraction of large-scale urban built-up areas. In addition, the existing methods rely too much on the quality of the data itself and the accuracy of the extracted features. However, the phenomenon of "different objects with the same spectrum" in remote sensing images is common, and there are a large number of natural features similar to town textures, resulting in a large number of wrong extractions.

Remote sensing satellite data with a spatial resolution of 10 meters not only has clear texture information, but also has smoother texture information in urban built-up areas, and the data can be obtained for free, but so far it has been rarely used in urban extraction. According to the characteristics of remote sensing satellite image data with a resolution of 10 meters, the present invention explores the texture features suitable for extracting towns, and proposes a method for constructing feature density by using superpixel segmentation to improve the accuracy of town extraction.

SUMMARY OF THE INVENTION

The present invention aims to solve the above problems of the prior art. A method is proposed, which has better completeness and correctness of the extracted results, especially in the aspect of completeness. The urban edge extracted by this method is highly consistent with the actual edge, which can achieve the effect of manual detection. The technical scheme of the present invention is as follows:

A method for automatically extracting urban built-up area from town texture feature density, which includes the following steps:

1) Obtain a remote sensing image in the visible light band with a resolution of 10 m, and perform cropping, splicing, and cloud removal operations to obtain an area to be extracted;

2) Constructing a grayscale image of the preprocessed remote sensing image, filtering out low-frequency textures based on Fourier transform on the grayscale image, and obtaining a grayscale image of high-frequency components;

3) for the grayscale image of the high frequency component obtained in step 2), calculate the grayscale co-occurrence matrix difference feature of several directions, and extract the rotational invariance texture feature;

4) the image that the rotation invariance texture feature obtained in step 3) is formed is mapped to the grayscale image of 0 to 255, and uses the Otsu threshold method OTSU to extract the feature area;

5) Multi-scale segmentation is performed on the image of the region to be extracted obtained in step 1) by a simple linear iterative clustering SLIC algorithm to obtain a multi-layer superpixel distribution image, and the superpixel is used as a unit to obtain a feature density image at each scale, Overlay all feature density images to build a refined feature density image.

Further, described step 1) is to carry out cropping, splicing, and cloud removal processing to the obtained image, specifically including: a remote sensing image of one scene usually cannot completely cover the area to be extracted, and it is necessary to download adjacent images of multiple scenes, and use the MapBased of ENVI software. The Mosaic function stitches the images to obtain a completely covered image, and then uses the Subset Data from ROIs function of the ENVI software to crop out the part outside the extraction area according to the range of the area to be extracted (vector file .shp). It is very common for ground objects in remote sensing images to be occluded by clouds, and the part covered by clouds in the image to be processed cannot obtain a good extraction effect, so cloud removal processing is required. The commonly used cloud removal method is as follows: manually identify the cloud area, look for images of other time periods without clouds in the area, crop out the image of this area, and stitch it into the image to be processed to obtain a clean image without cloud coverage. Preprocessing can be done using ENVI software.

Further, the step 2) filters out the low-frequency texture based on the Fourier transform of the grayscale image, and obtains the grayscale image of the high-frequency component, which specifically includes: using the two-dimensional fast Fourier transform FFT to convert the spatial domain image to frequency. The frequency domain image is filtered with a Gaussian low-pass filter, and the filtering result is converted to the spatial domain with the inverse fast Fourier transform to obtain the background image, and the original grayscale image and the background image are subtracted to obtain the suppressed low-frequency components. grayscale image.

Further, using the two-dimensional fast Fourier transform FFT to convert the spatial domain image to the frequency domain, the formula is: N represents the width and height of the image, respectively, (x, y) represents the spatial coordinates of the image, (u, v) represents the coordinates of the frequency image, f(x, y) is the input image, and F(u, v) is the converted image. The Fourier frequency-domain image of , filtered with a Gaussian low-pass filter in the frequency-domain image: F'(u,v)=F(u,v)·H(u,v) where H(x,y) is the Gaussian filter weight function, D represents the Gaussian filter window midpoint (u,v) to the Gaussian window center The Euclidean distance of , σ represents the degree of expansion of the Gaussian weight function, F(x, y) is the frequency domain image, F′(x, y) is the filtered frequency domain image, and the filtering result is converted by inverse fast Fourier transform Go to the spatial domain to obtain the background image f'(x,y), and subtract the original grayscale image f(x,y) from the background image f'(x,y) to obtain a grayscale image that suppresses low-frequency components: G(x,y)=f(x,y)-f'(x,y), G(x,y) is the finally obtained high-frequency component image.

Further, described step 3) uses the high-frequency component image obtained in the previous step, builds the gray level co-occurrence matrix difference feature of several directions, and extracts the rotation invariant texture feature, and the construction process is as follows;

The gray-level co-occurrence matrix is a matrix of size n×n, where n is the gray level of the image, and the gray-level co-occurrence matrix is defined as follows: M(i,j)=[m(i,j,Δx,Δy); i,j∈(1,2,…n)] The grayscale co-occurrence matrix M is an n×n matrix, and the element value of the i-th row and the j-th column is: within a certain size of the image in the window range of w×w Satisfy that the position offset is (Δx, Δy), and the gray value is the number of occurrences of the pixel pair of i and j, respectively, the difference formula is in m _(i,j) is the value of the elements in the i-th row and the j-th column of the gray-level co-occurrence matrix M, and constructs n1 kinds of co-occurrence matrices in n1 directions, so as to calculate the difference features of n1 kinds of directions, and select n1 kinds of directions The minimum value on the rotation-invariant texture feature is constructed. The n1 directions are (Δx,Δy)∈{(1,0),(2,0),(2,1),(1,1),(1,2),(-2,1),( -1,1),(1,2), (0,2),(0,1)}.

Further, the step 4) uses the Otsu threshold method OTSU to extract the feature area, which specifically includes:

1) the feature image is linearly mapped to the grayscale image of 0-255, and the grayscale histogram is counted, that is, the number of pixels of each grayscale;

2) Normalized histogram. That is, the number of pixels per gray level divided by the total number of pixels;

3) Calculate the maximum between-class variance. Let t be the grayscale threshold, count the proportion w ₀ of the pixels from 0 to t gray level in the whole image, count the average gray value u ₀ of the 0 to t gray level, and count the proportion of the pixels from t to 255 gray levels. The proportion w ₁ of the whole image, count the average gray value u ₁ from t to 255 gray levels, calculate g=w ₀ ×w ₁ × (u ₀ -u ₁ ) ² , set t to 0-255 in turn, The t value that maximizes g is the threshold value, the part larger than t is the feature area, and the part smaller than t is the other area.

Further, in the step 5), the superpixels obtained by the linear iterative clustering algorithm SLIC are uniformly arranged in a cell shape, which specifically includes: dividing the visible light band of the image n times with different scales based on the SLIC algorithm, and the scales are {N _i ,i∈1,2,3…n}, N _i represents the average number of pixels in the i-th segmentation superpixel, where N ₁ is the minimum segmentation scale, N _n is the maximum segmentation scale, from N ₁ to N _n in order The same step size S is increased, and the n superpixel layers obtained by segmentation are denoted as {L _i , i∈1,2,3…n}, and {P _ij ,j∈1,2,3…} denote the L _i The jth superpixel, Num _ij represents the number of pixels in P _ij ;

In step 3), based on the difference feature image with a window size of 23, the feature map obtained by the OSTU method is denoted as, and the feature density D _ij of all P _ij is calculated based on L _f : In the formula, Num _f represents the number of feature points in the range of P _ij , in the corresponding range of L _f , to obtain the feature density L' _i of each _Li , and finally accumulate L' _i to obtain the final feature density distribution L _d : It means that the pixels corresponding to each feature density image are added to obtain the total feature density L _d , the feature density is obtained by the above method, and finally the feature density image L _d is obtained by stacking each layer, and an appropriate threshold is set for the density image L _d , The part with higher density is the result of town extraction. The set threshold is between 0 and 1. The set threshold will be multiplied by the number of layers to be divided, and the part larger than the threshold will be extracted as towns.

The advantages and beneficial effects of the present invention are as follows:

Based on the advantages of frequency domain filtering, the present invention proposes a filtering method for remote sensing images, which can filter out the interference of low-frequency information in the original image. The difference between the intensity value and the gray value of the natural texture increases, so that the OSTU algorithm is used to obtain a more accurate feature area.

The present invention proposes to use superpixels as the calculation range of the most local feature density, and introduces the boundary line of ground objects as a limit, and different ground objects do not calculate the feature density in the same local range, and obtain a better feature density distribution. The multi-layer feature density increases the density difference inside and outside the urban built-up area, and at the same time refines the urban fringe. To get more accurate town information, using the empirical threshold to extract can get better results.

Through the extraction experiment of the MSIL1C visible light waveband remote sensing image of Sentinel-2, the extraction results of the urban built-up area extraction method proposed by the present invention have greater completeness and accuracy than the existing mainstream methods.

Description of drawings

Figure 1 is the obtained Sentinel-2MSIL1C brightness image

Figure 2 is the image after filtering the image in the frequency domain

Figure 3 shows the calculated differential rotation invariant features

Figure 4 shows the feature area of OSTU binarization

Figure 5 shows the feature density calculated by multi-layer superpixels

Figure 6 shows the results of town extraction by threshold method

Figure 7 method flow chart

Detailed ways

The technical solutions in the embodiments of the present invention will be described clearly and in detail below with reference to the accompanying drawings in the embodiments of the present invention. The described examples are partial examples of the invention.

The technical scheme that the present invention solves the above-mentioned technical problems is:

According to the characteristics of the ground objects in the 10-meter remote sensing image, the present invention selects the difference rotation invariance characteristics of 10 directions, and uses the OUST threshold method to obtain the characteristic area. In view of the problem that existing methods rely too much on feature extraction. The invention proposes a method for calculating local feature density. The method is based on the segmentation features of SLIC fitting the edge of objects, constructs multi-layer superpixels, calculates the feature density inside the superpixel, reduces the feature density outside the city, and realizes the realization of urban High-precision division of regional and non-urban areas.

The technical scheme of the method in this paper is as follows.

(1) Obtain remote sensing images in the visible light band with a resolution of 10 m, and perform relevant preprocessing;

(2) Constructing a grayscale image, filtering out low-frequency information based on Fourier transform, and obtaining a grayscale image of high-frequency components;

(3) for the high-frequency image obtained in step (2), calculate the gray level co-occurrence matrix difference feature of multiple directions, and extract the rotation invariant texture feature;

(4) the rotation-invariant texture feature image obtained in step (3) is mapped to the grayscale image of 0 to 255, and uses the Otsu threshold method (OTSU) to extract the feature area;

(5) Multi-scale segmentation of the true-color composite image of Sentinel-2 is performed by a simple linear iterative clustering (SLIC) algorithm to obtain a multi-layer superpixel distribution image. Taking superpixels as the unit, the feature density images at each scale are obtained, and all feature density images are superimposed to construct a refined feature density image.

(6) Finally, select the appropriate density threshold to extract the urban area.

The present invention is further described below:

(1) Crop, splicing, and cloud removal of the obtained image to obtain a clean and complete image in the visible light band.

(2) A grayscale image is constructed using the brightness information of the processed multi-band image, and then the image is filtered in the frequency domain. Convert the spatial domain image to the frequency domain using a two-dimensional fast Fourier transform (FFT), the formula is: N represents the width and height of the image, respectively, (x, y) represents the spatial coordinates of the image, (u, v) represents the coordinates of the frequency image, f(x, y) is the input image, and F(u, v) is the converted image. Fourier frequency domain image of . Filter the frequency domain image with a Gaussian low-pass filter: F'(u,v)=F(u,v)·H(u,v) where H(x,y) is the Gaussian filter weight function, D represents the Gaussian filter window midpoint (u,v) to the Gaussian window center The Euclidean distance of , σ represents the degree of expansion of the Gaussian weight function. F(x,y) is the frequency domain image, F'(x,y) is the filtered frequency domain image, and the filtering result is converted to the spatial domain by inverse fast Fourier transform, and the background image f'(x,y) is obtained ), subtract the original grayscale image f(x,y) from the background image f′(x,y) to obtain a grayscale image that suppresses low-frequency components: G(x,y)=f(x,y)-f'(x,y). G(x,y) is the final high-frequency component image.

(3) The grayscale co-occurrence matrix is a matrix of size n×n, where n is the grayscale level of the image (usually the grayscale level of a grayscale image is 256). The grayscale co-occurrence matrix is defined as follows: M(i,j)=[m(i,j,Δx,Δy); i,j∈(1,2,...n) The grayscale co-occurrence matrix M is an n×n matrix, The element value of the i-th row and the j-th column is: within a window range of a size w × w, the position offset is (Δx, Δy), and the gray value is the pixel pair of i and j, respectively. frequency. The present invention constructs a gray level co-occurrence matrix with 10 positional relationships, and the 10 kinds of positional relationships are (Δx, Δy)∈{(1,0),(2,0),(2,1),(1,1), (1,2),(-2,1),(-1,1),(1,2),(0,2),(0,1)}. The difference formula is in is the value of the element in the ith row and the jth column of the gray level co-occurrence matrix M. 10 co-occurrence matrices were constructed in 10 directions, and the difference features in 10 directions were calculated, and the minimum value in the 10 directions was selected to construct the rotation-invariant texture feature.

(4) Otsu thresholding method (Nobuyuki Otsu; OTSU) is a classic automatic thresholding algorithm, which divides the image into two parts, background and target, according to the grayscale characteristics of the image. The greater the inter-class variance between the background and the target, the greater the difference between the two parts of the image. When part of the target is mistakenly classified as the background or part of the background is mistakenly classified as the target, the difference between the two parts will become smaller. Therefore, the segmentation that maximizes the between-class variance means that the probability of misclassification is minimized. We map the rotation-invariant texture feature image obtained in step (3) to a grayscale image ranging from 0 to 255, and use OSTU to obtain two parts of the image. The part with large grayscale value is the feature region to be extracted.

(5) Simple linear iterative clustering (SLIC) is an efficient superpixel segmentation algorithm. The superpixels obtained by the SLIC algorithm are uniformly arranged in a cell shape. Superpixels can be adaptively deformed according to the boundary of the base map, and can be arranged according to the edge of the feature. The SLIC algorithm is calculated based on the CIELAB color space. The (l, a, b,) color of each pixel and the spatial coordinates (x, y) form a 5-dimensional vector, and the Euclidean distance is used to measure the distance between the two pixels. Similarity, the distance calculation formula is as follows: where i and j represent two pixels respectively, d _c represents the color distance, d _s represents the spatial distance, D represents the similarity of the two pixels, N _c is the desired maximum color distance within the class, and N _s is the desired intra-class color distance. maximum spatial distance. The distance D of the pixel from the cluster center determines the subordinate superpixel of the pixel. Among them, N _c and N _s are two important indicators that affect the clustering results. In general, N _s is determined by the number of cluster centers: N represents the total number of pixels in the image, and K is the number of cluster centers, which is the only input parameter of the algorithm. In the present invention, the average number of pixels of superpixels is used as a scale to measure the segmentation of superpixels.

Based on the SLIC algorithm, the present invention divides the visible light band of the image for n times with different scales, and the scales are {N _i , i∈1, 2, 3...n} respectively, and _Ni represents the average number of pixels in the i-th division of superpixels . Among them, N ₁ is the minimum segmentation scale, and N _n is the maximum segmentation scale. From N ₁ to N _n , the same step size S increases sequentially. The n superpixel layers obtained by segmentation are expressed as {L _i ,i∈1,2,3…n}, and {P _ij ,j∈1,2,3…} is used to denote the _jth superpixel of Li, Num _ij represents the number of pixels in P _ij . In step (3), based on the difference feature image with a window size of 23, the feature map obtained by the OSTU method is denoted as L _f . Calculate the eigendensity D _ij of all P _ij based on L _f : In the formula, Num _f represents the number of feature points in the corresponding range in L _f within the range of P _ij . The characteristic density L' _i of each _{Li is obtained, and finally L' i is} accumulated to obtain the final characteristic density distribution L _d : Indicates that the pixels corresponding to each feature density image are added to obtain the total feature density L _d . Set N ₁ to 300, S=10, n=5 to obtain 5-layer superpixel segmentation results, obtain feature densities through the above methods respectively, and finally superimpose each layer to obtain a feature density image L _d .

(6) Set an appropriate threshold based on the feature density image L _d obtained in step (5), and the part with higher density is the result of town extraction. Usually the threshold is set to 0.7×n to get better results, where n is the number of segmentation layers.

1. The example uses the image of the Sentinel-2A satellite launched by the European Space Agency of northern Beijing in September 2017, the serial number is "S2A_MSIL1C_20170922T025541_N0205_R032_T50". The present invention extracts the urban area based on the texture information of the image, and generally does not need to do radiometric calibration and atmospheric correction, but only needs to do cloud removal processing and cropping. The image shown in Figure 1 is obtained.

2. The obtained Sentinel-2A satellite data has multiple bands. We take the maximum value of the three bands of visible light to obtain the brightness image, as shown in Figure 2. After the luminance image is filtered in the frequency domain, a grayscale image with low-frequency information removed is obtained as shown in Figure 3. The value of σ in the example is set to M and N represent the width and height of the image, respectively.

3. Downgrade the grayscale of the image obtained in the previous step to 0-31, set the window size w to 32, calculate the texture dissimilarity in 10 directions, and take the minimum value to obtain the feature image, as shown in Figure 4 .

4. Map the feature image obtained in the previous step to the gray level of 0-255, select the threshold value through the Otsu threshold method, divide the image into two parts, and take the part with the larger gray value as the feature area. As shown in Figure 5

5. Through a simple linear iterative clustering algorithm, set the minimum segmentation scale to 300 pixels, set the step size to 10, and set the maximum segmentation scale to 340, and perform 5 segmentations in this way to obtain a 5-layer superpixel segmentation image. The feature density of each layer of superpixels is calculated based on superpixels, and the five layers of feature density images obtained by the calculation are superimposed to construct a total feature density image. The results are shown in Figure 6.

6. Finally, set the feature density threshold, the part larger than the threshold is urban area, and the threshold set for the instance is 3.5

The town extraction method proposed by the present invention can extract the precise town range. Compared with the results of human interpretation, the completeness rate of the above example is 0.96, and the correct rate is 0.92.

The above embodiments should be understood as only for illustrating the present invention and not for limiting the protection scope of the present invention. After reading the contents of the description of the present invention, the skilled person can make various changes or modifications to the present invention, and these equivalent changes and modifications also fall within the scope defined by the claims of the present invention.

Claims (7)

1. a kind of method that textural characteristics density in cities and towns automatically extracts cities and towns built-up areas, which comprises the following steps:

1) remote sensing image of 10m resolution ratio visible light wave range is obtained, and cut, splice, cloud is gone to operate, area to be extracted is obtained Domain；

2) grayscale image is constructed to pretreated remote sensing image, Fourier transformation is based on to grayscale image and filters out low frequency texture, Obtain the grayscale image of high fdrequency component；

3) for the grayscale image of high fdrequency component obtained in step 2), the gray level co-occurrence matrixes otherness in several directions is calculated Feature, and extract rotational invariance textural characteristics；

4) image of the rotational invariance textural characteristics composition obtained step 3), is mapped as 0 to 255 gray level image, uses Otsu threshold method OTSU extracts characteristic area；

5) multiple dimensioned point is carried out to area image to be extracted obtained in step 1) using simple linear iteration cluster SLIC algorithm It cuts, obtains multilayer super-pixel distributed image, as unit of super-pixel, calculate the characteristic density image under each scale, be superimposed institute There is characteristic density image, constructs the characteristic density image of fining.

2. the method that a kind of cities and towns textural characteristics density according to claim 1 automatically extracts cities and towns built-up areas, feature It is, the step 1) cuts the image of acquisition, splices, cloud removing, and specifically include: a scape remote sensing images are usually not Region to be extracted can be completely covered, needs to download the adjacent image of more scapes, uses the Map Based Mosaic function of ENVI software The image that energy stitching image is completely covered, then uses the Subset Data from ROIs function of ENVI software, according to The range (vector file .shp) in region to be extracted crops the part extracted outside region；Ground object is by cloud in remote sensing images The case where blocking is very universal, cannot obtain good extraction effect by the part that cloud covers in image to be processed, need to carry out Cloud removal processing, removes cloud method are as follows: manual identified cloud sector domain finds other cloudless period images of the region, cuts out the area Area image, and be spliced in image to be processed, obtain the clean image of cloudless covering；Pretreatment work is complete using ENVI software At.

3. the method that a kind of cities and towns textural characteristics density according to claim 1 automatically extracts cities and towns built-up areas, feature It is, the step 2) is based on Fourier transformation to grayscale image and filters out low frequency texture, obtains the grayscale image of high fdrequency component, has Body includes: that space area image is transformed into frequency domain using two-dimensional fast fourier transform FFT, with gauss low frequency filter to frequency Filter result is transformed into spatial domain with inverse fast fourier transform, obtains background video by rate domain images filter, by former ash degree shadow As doing additive operation with background video, the grayscale image for the low frequency component that is inhibited.

4. the method that a kind of cities and towns textural characteristics density according to claim 3 automatically extracts cities and towns built-up areas, feature It is, it is described that space area image is transformed into frequency domain, formula using two-dimensional fast fourier transform FFT are as follows:M, N respectively indicates the width and height of image, (x, y) Indicate the space coordinate of image, (u, v) indicates that the coordinate of frequency image, f (x, y) are input picture, and F (u, v) is that conversion obtains Fourier frequency area image, with gauss low frequency filter to frequency domain images filter: Wherein H (x, y) is gaussian filtering weight function to F ' (u, v)=F (u, v) H (u, v), and D indicates gaussian filtering window midpoint (u, v) To the Euclidean distance at Gauss window center, σ indicates the degree of Gauss weight function extension, and F (x, y) is frequency domain image, and F ' (x, y) is Filter result is transformed into spatial domain with inverse fast fourier transform by filtered frequency domain image, obtain background video f ' (x, Y), former grayscale image f (x, y) and background video f ' (x, y) are done into additive operation, the grayscale image for the low frequency component that is inhibited:G (x, y)=f (x, y)-f ' (x, y), G (x, y) is final The high fdrequency component image of acquisition.

5. the method that a kind of cities and towns textural characteristics density according to claim 1 automatically extracts cities and towns built-up areas, feature It is, the step 3) constructs the gray scale symbiosis in several directions using high fdrequency component grayscale image obtained in previous step Matrix difference feature, and rotational invariance textural characteristics are extracted, process is as follows；

The gray level co-occurrence matrixes are the matrixes that a size is n × n, and n is the gray level of image, and gray level co-occurrence matrixes define such as Under: M (i, j)=[m (i, j, Δ x, Δ y)；I, j ∈ (1,2 ... n)] matrix of the gray level co-occurrence matrixes M for n × n, the i-th row, the The element value of j column are as follows: meet positional shift in the window ranges that image some size be w × w for (Δ x, Δ y), and gray scale Value is respectively number of the pixel to appearance of i and j, and otherness formula is Whereinm_{(i, j)}For the i-th row of gray level co-occurrence matrixes M, the value of jth column element is constructed with n1 direction N1 kind gray level co-occurrence matrixes choose the minimum value building on n1 direction so that the otherness feature in n1 kind direction be calculated Rotational invariance textural characteristics.N1 direction be respectively (Δ x, Δ y) ∈ (1,0), (2,0), (2,1), (1,1), (1,2), (- 2,1), (- 1,1), (1,2), (0,2), (0,1) }.

6. the method that a kind of cities and towns textural characteristics density according to claim 1 automatically extracts cities and towns built-up areas, feature It is, the step 4) extracts characteristic area using Otsu threshold method OTSU, it specifically includes:

1) gray level image for being 0-255 by characteristic image Linear Mapping counts grey level histogram, i.e., the pixel of each gray level Number；

2) normalization histogram.The number of pixels of i.e. each gray level is divided by sum of all pixels；

3) maximum between-cluster variance is calculated.If t is gray threshold, statistics 0 accounts for the ratio w of whole sub-picture to the pixel of t gray level₀, system Meter 0 arrives the average gray value u of t gray level₀, the pixel of statistics t to 255 gray levels accounts for the ratio w of whole sub-picture₁, count t to 255 The average gray value u of gray level₁, calculate g=w₀×w₁×(u₀-u₁)², 0-255 successively is set by t, makes the maximum t value of g i.e. For threshold value, the part greater than t is characterized region, and the part less than t is other regions.

7. the method that a kind of cities and towns textural characteristics density according to claim 1 automatically extracts cities and towns built-up areas, feature It is, the super-pixel that the step 5) is obtained by linear iteraction clustering algorithm SLIC is specifically included at cellular evenly distributed: The segmentation of n times different scale is carried out to the visible light wave range of image based on SLIC algorithm, scale is respectively { N_i, i ∈ 1,2, 3...n }, N_iThe average pixel quantity of i-th segmentation super-pixel is represented, wherein N₁For smallest partition scale, N_nFor maximum fractionation ruler Degree, from N₁To N_nSuccessively with identical step-length S increase, the n super-pixel layer divided is expressed as { L_i, i ∈ 1,2,3...n }, With { P_ij, j ∈ 1,2,3... } and indicate L_iJ-th of super-pixel, Num_ijIndicate P_ijPixel quantity；

The otherness characteristic image for being 23 based on window size in step 3), obtains characteristic pattern with OSTU method and is denoted as, and is based on L_fMeter Calculate all P_ijCharacteristic density D_ij:Num in formula_fRepresent P_ijIn range, in L_fThe characteristic point of middle corresponding range Quantity, obtain each L_iCharacteristic density L '_i, finally add up L '_iObtain final characteristic density distribution L_d:It indicates for the corresponding pixel of each characteristic density image to be added, obtains total characteristic density L_d, by upper The method of stating obtains characteristic density, is finally superimposed each layer and obtains characteristic density image L_d, to density image L_dSuitable threshold value is set, The biggish part of density is cities and towns extraction as a result, set threshold value is between 0-1, and the threshold value of setting will be multiplied by segmentation The number of plies, the part greater than threshold value will be extracted as cities and towns.