CN111582146B - High-resolution remote sensing image city function partitioning method based on multi-feature fusion - Google Patents

High-resolution remote sensing image city function partitioning method based on multi-feature fusion Download PDF

Info

Publication number
CN111582146B
CN111582146B CN202010372618.6A CN202010372618A CN111582146B CN 111582146 B CN111582146 B CN 111582146B CN 202010372618 A CN202010372618 A CN 202010372618A CN 111582146 B CN111582146 B CN 111582146B
Authority
CN
China
Prior art keywords
word
image
feature
point
characteristic
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.)
Active
Application number
CN202010372618.6A
Other languages
Chinese (zh)
Other versions
CN111582146A (en
Inventor
孙伟伟
高子为
杨刚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ningbo University
Original Assignee
Ningbo University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Ningbo University filed Critical Ningbo University
Priority to CN202010372618.6A priority Critical patent/CN111582146B/en
Publication of CN111582146A publication Critical patent/CN111582146A/en
Application granted granted Critical
Publication of CN111582146B publication Critical patent/CN111582146B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/13—Satellite images
    • 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
    • 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/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G06F18/24155—Bayesian classification
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/28—Determining representative reference patterns, e.g. by averaging or distorting; Generating dictionaries
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/25—Determination of region of interest [ROI] or a volume of interest [VOI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Image Analysis (AREA)

Abstract

The invention relates to a high-resolution remote sensing image city function partitioning method based on multi-feature fusion, which comprises the following steps: step 1, preprocessing an image; step 2, distributing the characteristic value in each image to the visual word which is most similar to the characteristic value, and counting the corresponding word frequency of each visual word to form the visual word characteristic; constructing a multi-feature BoW visual dictionary; step 3, constructing an LDA probability topic model, and mining a high-dimensional semantic vector of the image by using the LDA probability topic model; step 4, training an SVM classifier according to the high-dimensional semantic vector obtained in the step 3; and 5, carrying out city function partitioning on the test set by using an SVM classifier. The invention has the beneficial effects that: POI data are introduced, so that the error score of the remote sensing data caused by same-object different spectrums and same-spectrum foreign objects is reduced; by comprehensively utilizing various characteristics of the image, including local characteristics, spectral characteristics, texture characteristics, surface temperature characteristics, spatial three-dimensional characteristics and POI characteristics, higher classification accuracy can be obtained under the condition that single characteristics of the image are not obvious.

Description

High-resolution remote sensing image city function partitioning method based on multi-feature fusion
Technical Field
The invention relates to the field of remote sensing image classification, in particular to a high-resolution remote sensing image urban function partitioning method based on multi-feature fusion.
Background
With the development of remote sensing technology, the time and space resolution of remote sensing images are continuously improved, and the data volume of the remote sensing images is increased rapidly. When a large amount of remote sensing data is faced, a method of manual visual interpretation needs to spend a great deal of time and labor for interpreting the remote sensing image. Therefore, how to automatically interpret the remote sensing image by using a computer becomes a hot research problem in the field of remote sensing. Meanwhile, Chinese economy and cities in a new normal state enter a new development stage, and a traditional city development mode faces a plurality of problems, so that new requirements and challenges are provided for city planning concepts, strategies and construction. In order to enhance the reasonable planning of urban functional areas and determine the most reasonable spatial layout of the urban functional areas, on the basis of improving the urban land utilization efficiency, various urban industries are gathered and exert the maximum efficiency, the urban land utilization efficiency is improved to a certain extent, and the effective implementation of a novel town strategy is ensured.
In order to overcome the gap between the low-level visual features and the high-level semantics of the remote sensing images, a method for modeling and describing land utilization scene semantics based on the middle-level features gradually gets wide attention. Especially, a Bag of Words model (Bag of Words model) in recent years has achieved great success in the application of image analysis and image classification, becomes a new and effective research idea for image content expression, and achieves certain results in remote sensing image land use scene classification. It represents the image as some local image blocks, each image block is represented as a sentence like words with different proportions, and the set of all words constitutes a visual dictionary.
The bag-of-words model has the characteristics of simple calculation, robustness to noise, illumination and local shielding and the like. However, the remote sensing image has abundant texture information and more local feature points, and the visual words obtained directly by the clustering method cannot necessarily reflect the scene features, so that the overall classification accuracy is not high. And for different building types in the high-resolution remote sensing image, the performance of each characteristic is not consistent. For example, some classes may be most suitable for classification by using spectral features due to the significance of the spectral features, while some classes have abundant textures, so that local features are more obvious. When the remote sensing image is only used for city function partitioning, the situations of 'same-object different-spectrum' and 'same-spectrum foreign matter' can occur, and the classification precision is reduced. It is obvious that the method of classifying different types of ground objects using the same feature is not applicable.
In summary, it is very important to provide a high-resolution remote sensing image city function partitioning method based on multi-feature fusion.
Disclosure of Invention
The invention aims to overcome the defects in the prior art, and provides a high-resolution remote sensing image city function partitioning method based on multi-feature fusion, which specifically comprises the following steps:
step 1, preprocessing an image, segmenting the image by using a proper grid, and selecting a training set and a testing set; calculating the POI type quantity proportion of each grid to obtain POI characteristics; respectively extracting local features, spectral features, textural features, earth surface temperature features and three-dimensional space features in a training set of the remote sensing image:
step 1.1, extracting local features of the remote sensing image: selecting an SURF algorithm, and solving a Hessian matrix corresponding to each pixel point (x, y) by adopting a Hessian matrix determinant approximate value image:
Figure GDA0002949883120000021
in the above formula, (x, y) is the pixel coordinate, and f (x, y) is the gray value of the coordinate point;
when the discrimination of the Hessian matrix obtains a local maximum value, judging that the current point is a brighter point or a darker point than other points in the surrounding neighborhood, and positioning the position of a key point; the discriminant of the Hessian matrix is:
Figure GDA0002949883120000022
performing box filtering calculation on the integral image by using a box filtering template added with the gradual change information to construct a scale space; comparing the det (H) of each pixel point with the det (H) of all adjacent points in the adjacent scale space, and when the det (H) of each pixel point is more than or less than all the adjacent points, determining that the point is an extreme point, and recording the position of the extreme point as a feature point;
constructing a 60-degree fan-shaped sliding window in a circular area with the characteristic point as the center of a circle and 6s as the radius; the s is a scale space of the feature points; traversing the whole circular area in a rotating way by using the radian of 0.2, and calculating Haar wavelet characteristic values in a fan-shaped sliding window; selecting the direction with the maximum sum of the Haar wavelet characteristic values as the main direction of the characteristic point; the solving method of the Haar wavelet characteristic sum is to accumulate the Harr wavelet characteristic values dx and dy of the image to obtain a vector (m)ω,θω)：
mω＝∑ωdx+∑ωdy (3)
θω＝arctan(∑ωdx/∑ωdy) (4)
The main direction θ is a direction corresponding to the maximum Harr wavelet feature value of the image, namely:
θ＝θω|max{mω} (5)
taking a square frame with a theta direction at the characteristic point, wherein the side length of the frame is 20s, and s is a scale space of the characteristic point; dividing a square frame with a theta direction into 4 multiplied by 4 sub-regions, and counting the sum of the Harr wavelet characteristic values of each sub-region in the horizontal direction, the sum of the absolute values in the horizontal direction, the sum of the vertical direction and the sum of the absolute values in the vertical direction; the SURF feature descriptor is composed of 4 multiplied by 64-dimensional feature vectors;
step 1.2, extracting spectral characteristics of the remote sensing image, and calculating the mean value and standard deviation of each wave band, wherein the implementation mode is as follows:
Figure GDA0002949883120000031
Figure GDA0002949883120000032
in the above formulas (6) and (7), n is the total number of pixels in the grid, viThe gray value of the ith pixel in the wave band;
step 1.3, extracting texture features of the remote sensing image, selecting an LBP operator, defining the LBP operator in a window of 3 multiplied by 3 pixels:
Figure GDA0002949883120000033
in the above formula, (x)c，yc) Representing the central point of each 3 × 3 area, wherein p represents the p-th pixel point except the central pixel point in the 3 × 3 window; i (c) represents the gray value of the central pixel point, and I (p) represents the gray value of the p-th pixel point in the field;
taking the window center pixel as a threshold, comparing the center pixel value with the gray values of the adjacent 8 pixels: if the neighboring pixel value is greater than the center pixel value, then the position is marked as 1; otherwise, the flag is 0:
Figure GDA0002949883120000034
in the above formula, x is the central pixel value; obtaining an 8-bit binary number, and using the 8-bit binary number as an LBP value of a window center pixel point to reflect the texture information of the 3 x 3 pixel window;
step 1.4, extracting surface temperature characteristics, and calculating the surface specific index according to NDVI:
Figure GDA0002949883120000035
in the above formula, epsilon is the earth surface emissivity; NDVI is the normalized vegetation index; the surface temperature is calculated according to a single window algorithm, which has the following formula:
Ts＝{a·(1-C-D)[b·(1-C-D)+C+D]·T-D·Ta}/C (11)
in the above formula, TsIs the surface temperature; a and b are empirical coefficients, wherein a is-67.35535 and b is 0.458608; t is the brightness temperature; t isaIs the average temperature of action of the atmosphere; the calculation formula for C and D is as follows:
C＝τ·ε (12)
D＝(1-τ)·[1+τ·(1-ε)] (13)
in the above formulas (12) to (13), τ is the atmospheric transmittance, and ε is the surface emissivity;
step 1.5, extracting three-dimensional space characteristics: three-dimensional space point cloud data are obtained through a LiDAR technology, and then a digital earth surface model DSM characteristic is generated through a Krigin space interpolation algorithm;
step 2, distributing the characteristic value in each image to the visual word which is most similar to the characteristic value, and counting the corresponding word frequency of each visual word to form the visual word characteristic; constructing a multi-feature BoW visual dictionary:
step 2.1, constructing various feature vocabulary lists: respectively carrying out K-means clustering on the local features, the spectral features, the textural features, the earth surface temperature features and the three-dimensional space features extracted in the step 1; each clustering center is a word, and all clustering centers are converged into a word list;
step 2.2, by calculating the distance between the characteristics of each grid and each word, assigning the characteristics to the word with the closest distance in the word list;
step 2.3, counting word frequencies corresponding to the words to generate a K-dimensional feature vector, splicing the K-dimensional feature vectors in a stacking mode, and representing each grid by using a multi-dimensional vector:
fi＝{suri,spei,lbpi,lsti,dsmi,poii} (14)
in the above formula, suriFor local feature vectors, speiAs spectral feature vector, lbpiAs texture feature vectors, lstiAs surface temperature feature vector, dsmiFor three-dimensional spatial feature vectors, poiiAs feature vectors of POI
Each mesh is described as a document with a characteristic vocabulary:
Doci＝{wordsur,wordspe,wordlbp,wordlst,worddsm,wordpoi}i (15)
step 3, constructing an LDA probability topic model, mining high-dimensional semantic vectors of the image by using the LDA probability topic model, and distributing the probability of each feature vector by extracting the high-dimensional semantic vectors contained in each feature vector by the LDA probability topic model:
step 3.1, giving each document in the document set according to the form of probability distribution by using an LDA probability topic model:
P(w|d)＝P(w|t)×P(t|d) (16)
in the above formula, w is a word, d is a document, and t is a theme; with topic t as intermediate layer, passing through two vectors
Figure GDA0002949883120000041
Give P (w | t) and P (t | d), θ, respectivelydRepresenting the probability that for an individual document D in each aggregate set D corresponds to a different individual topic t,
Figure GDA0002949883120000042
representing the probability vectors that different words are generated for the individual subjects t in each aggregate set TA;
step 3.2, the learning process of the LDA probability topic model is as follows:
for all d and t, first give θ randomlyd,
Figure GDA0002949883120000043
Assigning; for a particular document dsThe ith word w iniIf the word wiCorresponding topic is tjThe formula (15) is rewritten as:
Pj(wi|ds)＝P(wi|tj)×P(tj|ds) (17)
enumerating the subject t in the total set TA to obtain all Pj(wi|ds) (ii) a According to Pj(wi|ds) Result of (a) is dsThe ith word w iniSelecting a topic t, wiCorresponding subject t gets Pj(wi|dS) Topic t with the highest probabilityj；
An iterative process is defined as: performing P (w | D) calculation once on all w in all the documents D in the document total set D, and reselecting a theme;
the iterative process is repeated until thetad,
Figure GDA0002949883120000044
Convergence and final output of the estimated parameter thetad,
Figure GDA0002949883120000045
Obtaining a theme of each word and a high-dimensional semantic vector of each theme;
step 4, training an SVM classifier according to the high-dimensional semantic vector obtained in the step 3;
and 5, carrying out city function partitioning on the test set by using an SVM classifier.
Preferably, the POI feature in step 1 is used to re-classify the POI data according to the location name information.
Preferably, the SURF algorithm in step 1.1 constructs a scale pyramid by changing the size of the box filter while keeping the image size unchanged.
Preferably, in step 1.4, the earth surface temperature features are extracted in a satellite remote sensing data inversion mode.
Preferably, the SURF descriptors of step 1.1 have scale and rotation invariance, as well as invariance to changes in illumination.
Preferably, the LDA probabilistic topic model in step 3 is an unsupervised bayesian model.
The invention has the beneficial effects that: according to the method, POI data are introduced, so that the misclassification of the remote sensing data caused by same-object different spectrums and same-spectrum foreign objects is reduced; the method comprehensively utilizes various characteristics of the image, including local characteristics, spectral characteristics, texture characteristics, surface temperature characteristics, space three-dimensional characteristics and POI characteristics, and can obtain higher classification accuracy under the condition that single characteristics of the image are not obvious.
Drawings
FIG. 1 is a technical flow chart of the present invention.
Detailed Description
The present invention will be further described with reference to the following examples. The following examples are set forth merely to aid in the understanding of the invention. It should be noted that, for a person skilled in the art, several modifications can be made to the invention without departing from the principle of the invention, and these modifications and modifications also fall within the protection scope of the claims of the present invention.
A high-resolution remote sensing image city function partitioning method based on multi-feature fusion is disclosed, a technical flow chart of which is shown in figure 1, and the method specifically comprises the following steps:
step 1, preprocessing an image, segmenting the image by using a proper grid, and selecting a training set and a testing set; calculating the POI type quantity proportion of each grid to obtain POI characteristics; respectively extracting local features, spectral features, textural features, earth surface temperature features and three-dimensional space features in a training set of the remote sensing image:
step 1.1, extracting local features of the remote sensing image: selecting an SURF algorithm, and solving a Hessian matrix corresponding to each pixel point (x, y) by adopting a Hessian matrix determinant approximate value image:
Figure GDA0002949883120000051
in the above formula, (x, y) is the pixel coordinate, and f (x, y) is the gray value of the coordinate point;
when the discrimination of the Hessian matrix obtains a local maximum value, judging that the current point is a brighter point or a darker point than other points in the surrounding neighborhood, and positioning the position of a key point; the discriminant of the Hessian matrix is:
Figure GDA0002949883120000061
performing box filtering calculation on the integral image by using a box filtering template added with the gradual change information to construct a scale space; comparing the det (H) of each pixel point with the det (H) of all adjacent points in the adjacent scale space, and when the det (H) of each pixel point is more than or less than all the adjacent points, determining that the point is an extreme point, and recording the position of the extreme point as a feature point;
constructing a 60-degree fan-shaped sliding window in a circular area with the characteristic point as the center of a circle and 6s as the radius; the s is a scale space of the feature points; traversing the whole circular area in a rotating way by using the radian of 0.2, and calculating Haar wavelet characteristic values in a fan-shaped sliding window; selecting the direction with the maximum sum of the Haar wavelet characteristic values as the main direction of the characteristic point; the solving method of the Haar wavelet characteristic sum is to accumulate the Harr wavelet characteristic values dx and dy of the image to obtain a vector (m)ω,θω)：
mω＝∑ωdx+∑ωdy (3)
θω＝arctan(∑ωdx/∑ωdy) (4)
The main direction θ is a direction corresponding to the maximum Harr wavelet feature value of the image, namely:
θ＝θω|max{mω} (5)
taking a square frame with a theta direction at the characteristic point, wherein the side length of the frame is 20s, and s is a scale space of the characteristic point; dividing a square frame with a theta direction into 4 multiplied by 4 sub-regions, and counting the sum of the Harr wavelet characteristic values of each sub-region in the horizontal direction, the sum of the absolute values in the horizontal direction, the sum of the vertical direction and the sum of the absolute values in the vertical direction; the SURF feature descriptor is composed of 4 multiplied by 64-dimensional feature vectors;
step 1.2, extracting spectral characteristics of the remote sensing image, and calculating the mean value and standard deviation of each wave band, wherein the implementation mode is as follows:
Figure GDA0002949883120000062
Figure GDA0002949883120000063
in the above formulas (6) and (7), n is the total number of pixels in the grid, viIs at the same timeThe gray value of the ith pixel in the wave band;
step 1.3, extracting texture features of the remote sensing image, selecting an LBP operator, defining the LBP operator in a window of 3 multiplied by 3 pixels:
Figure GDA0002949883120000064
in the above formula, (x)c，yc) Representing the central point of each 3 × 3 area, wherein p represents the p-th pixel point except the central pixel point in the 3 × 3 window; i (c) represents the gray value of the central pixel point, and I (p) represents the gray value of the p-th pixel point in the field;
taking the window center pixel as a threshold, comparing the center pixel value with the gray values of the adjacent 8 pixels: if the neighboring pixel value is greater than the center pixel value, then the position is marked as 1; otherwise, the flag is 0:
Figure GDA0002949883120000071
in the above formula, x is the central pixel value; obtaining an 8-bit binary number, and using the 8-bit binary number as an LBP value of a window center pixel point to reflect the texture information of the 3 x 3 pixel window;
step 1.4, extracting surface temperature characteristics, and calculating the surface specific index according to NDVI:
Figure GDA0002949883120000072
in the above formula, epsilon is the earth surface emissivity; NDVI is the normalized vegetation index; the surface temperature is calculated according to a single window algorithm, which has the following formula:
Ts＝{a·(1-C-D)[b·(1-C-D)+C+D]·T-D·Ta}/C (11)
in the above formula, TsIs the surface temperature; a and b are empirical coefficients, wherein a is-67.35535 and b is 0.458608; t is the brightness temperature; t isaIs the average temperature of action of the atmosphere;the calculation formula for C and D is as follows:
C＝τ·ε (12)
D＝(1-τ)·[1+τ·(1-ε)] (13)
in the above formulas (12) to (13), τ is the atmospheric transmittance, and ε is the surface emissivity;
step 1.5, extracting three-dimensional space characteristics: three-dimensional space point cloud data are obtained through a LiDAR technology, and then a digital earth surface model DSM characteristic is generated through a Krigin space interpolation algorithm;
step 2, distributing the characteristic value in each image to the visual word which is most similar to the characteristic value, and counting the corresponding word frequency of each visual word to form the visual word characteristic; constructing a multi-feature BoW visual dictionary:
step 2.1, constructing various feature vocabulary lists: respectively carrying out K-means clustering on the local features, the spectral features, the textural features, the earth surface temperature features and the three-dimensional space features extracted in the step 1; each clustering center is a word, and all clustering centers are converged into a word list;
step 2.2, by calculating the distance between the characteristics of each grid and each word, assigning the characteristics to the word with the closest distance in the word list;
step 2.3, counting word frequencies corresponding to the words to generate a K-dimensional feature vector, splicing the K-dimensional feature vectors in a stacking mode, and representing each grid by using a multi-dimensional vector:
fi＝{suri,spei,lbpi,lsti,dsmi,poii} (14)
in the above formula, suriFor local feature vectors, speiAs spectral feature vector, lbpiAs texture feature vectors, lstiAs surface temperature feature vector, dsmiFor three-dimensional spatial feature vectors, poiiAs feature vectors of POI
Each mesh is described as a document with a characteristic vocabulary:
Doci＝{wordsur,wordspe,wordlbp,wordlst,worddsm,wordpoi}i (15)
step 3, constructing an LDA probability topic model, mining high-dimensional semantic vectors of the image by using the LDA probability topic model, and distributing the probability of each feature vector by extracting the high-dimensional semantic vectors contained in each feature vector by the LDA probability topic model: in the semantic model classification, the probability of the same word appearing in different subject contexts is different. And extending to an image classification neighborhood, wherein the probability of each characteristic value appearing in different class topics is different, and introducing an LDA topic model in order to determine the corresponding relation between the appearance probability of each characteristic value and the class topics.
Step 3.1, giving each document in the document set according to the form of probability distribution by using an LDA probability topic model:
P(w|d)＝P(w|t)×P(t|d) (16)
in the above formula, w is a word, d is a document, and t is a theme; with topic t as intermediate layer, passing through two vectors
Figure GDA0002949883120000081
Give P (w | t) and P (t | d), θ, respectivelydRepresenting the probability that for an individual document D in each aggregate set D corresponds to a different individual topic t,
Figure GDA0002949883120000082
representing the probability vectors that different words are generated for the individual subjects t in each aggregate set TA;
step 3.2, the learning process of the LDA probability topic model is as follows:
for all d and t, first give θ randomlyd,
Figure GDA0002949883120000085
Assigning; for a particular document dsThe ith word w iniIf the word wiCorresponding topic is tjThe formula (15) is rewritten as:
Pj(wi|ds)＝P(wi|tj)×P(tj|ds) (17)
enumerating the subject t in the total set TA to obtain all Pj(wi|ds) (ii) a According to Pj(wi|ds) Result of (a) is dsThe ith word w iniSelecting a topic t, wiCorresponding subject t gets Pj(wi|ds) Topic t with the highest probabilityj；
An iterative process is defined as: performing P (w | D) calculation once on all w in all the documents D in the document total set D, and reselecting a theme;
the iterative process is repeated until thetad,
Figure GDA0002949883120000083
Convergence and final output of the estimated parameter thetad,
Figure GDA0002949883120000084
Obtaining a theme of each word and a high-dimensional semantic vector of each theme;
step 4, training an SVM classifier according to the high-dimensional semantic vector obtained in the step 3; and labeling each grid by visual interpretation according to the classification standard of urban land classification and planning construction land standard of the samples in the training set to obtain a label set, training an SVM classifier by combining potential semantic feature vectors contained in each training sample, and classifying the test samples.
And 5, carrying out city function partitioning on the test set by using an SVM classifier.
And the POI characteristics in the step 1 are used for reclassifying the POI data according to the place name information. The introduced POI data solves the problem of building functions that the remote sensing data can not be extracted due to 'same-object different spectrum' and 'same-spectrum foreign matter'.
In the step 1.1, the SURF algorithm constructs the scale pyramid by changing the size of the box filter while keeping the size of the image unchanged. The SURF algorithm is an improvement over the SIFT algorithm.
And in the step 1.4, the earth surface temperature characteristics are extracted in a satellite remote sensing data inversion mode.
The SURF feature descriptor described in step 1.1 has scale and rotation invariance, and also has invariance to changes in illumination.
The LDA probability topic model in the step 3 is an unsupervised Bayesian model.
In the experimental data research of the example, the region of Wanke city in the town sea area of Ningbo at 11 th of 9 th of 2019 is selected, and the data includes score 2 data, TM data, LiDAR data and Baidu POI data.
According to 11/1/2017, the city function is divided into the following parts: the method comprises the following steps of classifying test sets by using trained classifiers, wherein the test sets comprise eight types of ground object types including residential sites, public management and public service sites, commercial service facility sites, industrial sites, logistics storage sites, road and traffic facility sites, public facility sites and green and square site sites.

Claims (4)

1. A high-resolution remote sensing image city function partitioning method based on multi-feature fusion is characterized by comprising the following steps:
step 1, preprocessing an image, segmenting the image by using a proper grid, and selecting a training set and a testing set; calculating the POI type quantity proportion of each grid to obtain POI characteristics; respectively extracting local features, spectral features, textural features, surface temperature features and three-dimensional space features in a training set of the remote sensing image;
step 1.1, extracting local features of the remote sensing image: selecting an SURF algorithm, and solving a Hessian matrix corresponding to each pixel point (x, y) by adopting a Hessian matrix determinant approximate value image:
Figure FDA0002949883110000011
in the above formula, (x, y) is the pixel coordinate, and f (x, y) is the gray value of the coordinate point;
when the discrimination of the Hessian matrix obtains a local maximum value, judging that the current point is a brighter point or a darker point than other points in the surrounding neighborhood, and positioning the position of a key point; the discriminant of the Hessian matrix is:
Figure FDA0002949883110000012
performing box filtering calculation on the integral image by using a box filtering template added with the gradual change information to construct a scale space; comparing the det (H) of each pixel point with the det (H) of all adjacent points in the adjacent scale space, and when the det (H) of each pixel point is more than or less than all adjacent points, determining that the pixel points with the det (H) values more than or less than all the adjacent points are extreme points, and recording the positions of the extreme points as feature points;
constructing a 60-degree fan-shaped sliding window in a circular area with the characteristic point as the center of a circle and 6s as the radius; the s is a scale space of the feature points; traversing the whole circular area in a rotating way by using the radian of 0.2, and calculating Haar wavelet characteristic values in a fan-shaped sliding window; selecting the direction with the maximum sum of the Haar wavelet characteristic values as the main direction of the characteristic point; the solving method of the Haar wavelet characteristic sum is to accumulate the Harr wavelet characteristic values dx and dy of the image to obtain a vector (m)ω，θω)：
mω＝∑ωdx+∑ωdy (3)
θω＝arctan(∑ωdx/∑ωdy) (4)
The main direction θ is a direction corresponding to the maximum Harr wavelet feature value of the image, namely:
θ＝θω|max{mω} (5)
taking a square frame with a theta direction at the characteristic point, wherein the side length of the frame is 20s, and s is a scale space of the characteristic point; dividing a square frame with a theta direction into 4 multiplied by 4 sub-regions, and counting the sum of the Harr wavelet characteristic values of each sub-region in the horizontal direction, the sum of the absolute values in the horizontal direction, the sum of the vertical direction and the sum of the absolute values in the vertical direction; the SURF feature descriptor is composed of 4 multiplied by 64-dimensional feature vectors;
step 1.2, extracting spectral characteristics of the remote sensing image, and calculating the mean value and standard deviation of each wave band, wherein the implementation mode is as follows:
Figure FDA0002949883110000021
Figure FDA0002949883110000022
in the above formulas (6) and (7), n is the total number of pixels in the grid, viThe gray value of the ith pixel in the wave band;
step 1.3, extracting texture features of the remote sensing image, selecting an LBP operator, defining the LBP operator in a window of 3 multiplied by 3 pixels:
Figure FDA0002949883110000023
in the above formula, (x)c，yc) Representing the central point of each 3 × 3 area, wherein p represents the p-th pixel point except the central pixel point in the 3 × 3 window; i (c) represents the gray value of the central pixel point, and I (p) represents the gray value of the p-th pixel point in the field;
taking the window center pixel as a threshold, comparing the center pixel value with the gray values of the adjacent 8 pixels: if the neighboring pixel value is greater than the center pixel value, then the position is marked as 1; otherwise, the flag is 0:
Figure FDA0002949883110000024
in the above formula, x is the central pixel value; obtaining an 8-bit binary number, and using the 8-bit binary number as an LBP value of a window center pixel point to reflect the texture information of the 3 x 3 pixel window;
step 1.4, extracting surface temperature characteristics, and calculating the surface specific index according to NDVI:
Figure FDA0002949883110000025
in the above formula, epsilon is the earth surface emissivity; NDVI is the normalized vegetation index; the surface temperature is calculated according to a single window algorithm, which has the following formula:
Ts＝{a·(1-C-D)[b·(1-C-D)+C+D]·T-D·Ta}/C (11)
in the above formula, TsIs the surface temperature; a and b are empirical coefficients, wherein a is-67.35535 and b is 0.458608; t is the brightness temperature; t isaIs the average temperature of action of the atmosphere; the calculation formula for C and D is as follows:
C＝τ·ε (12)
D＝(1-τ)·[1+τ·(1-ε)] (13)
in the above formulas (12) to (13), τ is the atmospheric transmittance, and ε is the surface emissivity;
step 1.5, extracting three-dimensional space characteristics: three-dimensional space point cloud data are obtained through a LiDAR technology, and then a digital earth surface model DSM characteristic is generated through a Krigin space interpolation algorithm;
step 2, distributing the characteristic value in each image to the visual word which is most similar to the characteristic value, and counting the corresponding word frequency of each visual word to form the visual word characteristic; constructing a multi-feature BoW visual dictionary:
step 2.1, constructing various feature vocabulary lists: respectively carrying out K-means clustering on the local features, the spectral features, the textural features, the earth surface temperature features and the three-dimensional space features extracted in the step 1; each clustering center is a word, and all clustering centers are converged into a word list;
step 2.2, by calculating the distance between the characteristics of each grid and each word, assigning the characteristics to the word with the closest distance in the word list;
step 2.3, counting word frequencies corresponding to the words to generate a K-dimensional feature vector, splicing the K-dimensional feature vectors in a stacking mode, and representing each grid by using a multi-dimensional vector:
fi＝{suri，spei，lbpi，lsti，dsmi，poii} (14)
in the above formula, suriFor local feature vectors, speiAs spectral feature vector, lbpiAs texture feature vectors, lstiAs surface temperature feature vector, dsmiFor three-dimensional spatial feature vectors, poiiIs POI feature vector;
each mesh is described as a document with a characteristic vocabulary:
Doci＝{wordsur，wordspe，wordlbp，wordlst，worddsm，wordpoi}i (15)
step 3, constructing an LDA probability topic model, mining high-dimensional semantic vectors of the image by using the LDA probability topic model, and distributing the probability of each feature vector by extracting the high-dimensional semantic vectors contained in each feature vector by using the LDA probability topic model;
step 3.1, giving each document in the document set according to the form of probability distribution by using an LDA probability topic model:
P(w|d)＝P(w|t)×P(t|d) (16)
in the above formula, w is a word, d is a document, and t is a theme; with topic t as intermediate layer, passing through two vectors
Figure FDA0002949883110000031
Give P (w | t) and P (t | d), θ, respectivelydRepresenting the probability that for an individual document D in each aggregate set D corresponds to a different individual topic t,
Figure FDA0002949883110000032
representing the probability vectors that different words are generated for the individual subjects t in each aggregate set TA;
step 3.2, the learning process of the LDA probability topic model is as follows:
for all d and t, first give θ randomlyd，
Figure FDA0002949883110000033
Assigning; for a particular document dsThe ith word w iniIf the word wiCorresponding topic is tjThe formula (15) is rewritten as:
Pj(wi|ds)＝P(wi|tj)×P(tj|ds) (17)
enumerating the subject t in the total set TA to obtain all Pj(wi|ds) (ii) a According to Pj(wi|ds) Result of (a) is dsThe ith word w iniSelecting a topic t, wiCorresponding subject t gets Pj(wi|ds) Topic t with the highest probabilityj；
An iterative process is defined as: performing P (w | D) calculation once on all w in all the documents D in the document total set D, and reselecting a theme;
the iterative process is repeated until thetad，
Figure FDA0002949883110000041
Convergence and final output of the estimated parameter thetad，
Figure FDA0002949883110000042
Obtaining a theme of each word and a high-dimensional semantic vector of each theme;
step 4, training an SVM classifier according to the high-dimensional semantic vector obtained in the step 3;
and 5, carrying out city function partitioning on the test set by using an SVM classifier.
2. The urban functional partitioning method for the high-resolution remote sensing images based on the multi-feature fusion as claimed in claim 1, wherein: and in the step 1, POI data is combined to be used as a geographic space big data characteristic to participate in city functional partitioning.
3. The urban functional partitioning method for the high-resolution remote sensing images based on the multi-feature fusion as claimed in claim 1, wherein: and 2, constructing a BoW dictionary in a multi-feature library cooperation mode.
4. The urban functional partitioning method for the high-resolution remote sensing images based on the multi-feature fusion as claimed in claim 1, wherein: the LDA probability topic model in the step 3 is an unsupervised Bayesian model.
CN202010372618.6A 2020-05-06 2020-05-06 High-resolution remote sensing image city function partitioning method based on multi-feature fusion Active CN111582146B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010372618.6A CN111582146B (en) 2020-05-06 2020-05-06 High-resolution remote sensing image city function partitioning method based on multi-feature fusion

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010372618.6A CN111582146B (en) 2020-05-06 2020-05-06 High-resolution remote sensing image city function partitioning method based on multi-feature fusion

Publications (2)

Publication Number Publication Date
CN111582146A CN111582146A (en) 2020-08-25
CN111582146B true CN111582146B (en) 2021-04-13

Family

ID=72111986

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010372618.6A Active CN111582146B (en) 2020-05-06 2020-05-06 High-resolution remote sensing image city function partitioning method based on multi-feature fusion

Country Status (1)

Country Link
CN (1) CN111582146B (en)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112183879A (en) * 2020-10-14 2021-01-05 北京大学 Method and device for classifying urban functional areas, electronic equipment and storage medium
CN112364776A (en) * 2020-11-12 2021-02-12 中国地质大学（武汉） Urban functional area recognition model establishing and classifying method and device
CN112884635A (en) * 2021-01-25 2021-06-01 中交广州航道局有限公司 Submarine environment visualization method and device based on ROV carrying dual-frequency forward-looking sonar
CN113447771A (en) * 2021-06-09 2021-09-28 上海交通大学 Partial discharge pattern recognition method based on SIFT-LDA characteristics
CN113657324B (en) * 2021-08-24 2024-06-21 速度科技股份有限公司 Urban functional area identification method based on remote sensing image ground object classification
CN113988198B (en) * 2021-11-03 2024-05-14 中国地质大学（武汉） Multi-scale city function classification method based on landmark constraint
CN114399685B (en) * 2022-03-25 2022-06-21 航天宏图信息技术股份有限公司 Remote sensing monitoring and evaluating method and device for forest pest and disease damage
CN115299937B (en) * 2022-09-16 2023-08-11 东联信息技术有限公司 Intelligence detection platform that falls down

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102622607A (en) * 2012-02-24 2012-08-01 河海大学 Remote sensing image classification method based on multi-feature fusion
US8369593B2 (en) * 2007-12-21 2013-02-05 Siemens Medical Solutions Usa, Inc. Systems and methods for robust learning based annotation of medical radiographs
CN106897417A (en) * 2017-02-22 2017-06-27 东南大学 A kind of construction method of the city space holographic map based on the fusion of multi-source big data
CN107247938A (en) * 2017-06-08 2017-10-13 中国科学院遥感与数字地球研究所 A kind of method of high-resolution remote sensing image City Building function classification
CN108665149A (en) * 2018-04-23 2018-10-16 深圳市城市规划设计研究院有限公司 A kind of measuring method and system based on city function degree of mixing
CN108764193A (en) * 2018-06-04 2018-11-06 北京师范大学 Merge the city function limited region dividing method of POI and remote sensing image
CN109325085A (en) * 2018-08-08 2019-02-12 中南大学 A kind of urban land identification of function and change detecting method

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8369593B2 (en) * 2007-12-21 2013-02-05 Siemens Medical Solutions Usa, Inc. Systems and methods for robust learning based annotation of medical radiographs
CN102622607A (en) * 2012-02-24 2012-08-01 河海大学 Remote sensing image classification method based on multi-feature fusion
CN106897417A (en) * 2017-02-22 2017-06-27 东南大学 A kind of construction method of the city space holographic map based on the fusion of multi-source big data
CN107247938A (en) * 2017-06-08 2017-10-13 中国科学院遥感与数字地球研究所 A kind of method of high-resolution remote sensing image City Building function classification
CN108665149A (en) * 2018-04-23 2018-10-16 深圳市城市规划设计研究院有限公司 A kind of measuring method and system based on city function degree of mixing
CN108764193A (en) * 2018-06-04 2018-11-06 北京师范大学 Merge the city function limited region dividing method of POI and remote sensing image
CN109325085A (en) * 2018-08-08 2019-02-12 中南大学 A kind of urban land identification of function and change detecting method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于高分一号影像和POI数据的土地功能分类方法研究;高晨;《中国博士学位论文全文数据库.基础科学辑》;20190515(第05期);第7页第1段，第3.1节，第28页第3段 *
面向高分辨率遥感影像场景语义理解的概率主题模型研究;朱祺琪;《中国博士学位论文全文数据库.基础科学辑》;20190615(第06期);摘要，第14页第1段，第19页第2.2.2节 *

Also Published As

Publication number Publication date
CN111582146A (en) 2020-08-25

Similar Documents

Publication Publication Date Title
CN111582146B (en) High-resolution remote sensing image city function partitioning method based on multi-feature fusion
Zhang et al. An object-based convolutional neural network (OCNN) for urban land use classification
Nahhas et al. Deep learning approach for building detection using lidar–orthophoto fusion
Du et al. Semantic classification of urban buildings combining VHR image and GIS data: An improved random forest approach
Zhu et al. Scene classification based on the sparse homogeneous–heterogeneous topic feature model
CN104915636B (en) Remote sensing image road recognition methods based on multistage frame significant characteristics
US11854244B2 (en) Labeling techniques for a modified panoptic labeling neural network
Luo et al. Cross-spatiotemporal land-cover classification from VHR remote sensing images with deep learning based domain adaptation
CN110399819A (en) A kind of remote sensing image residential block extraction method based on deep learning
CN108427919B (en) Unsupervised oil tank target detection method based on shape-guided saliency model
Jiao et al. A survey of road feature extraction methods from raster maps
Zhang et al. A dual attention neural network for airborne LiDAR point cloud semantic segmentation
Zhu et al. High-resolution remote sensing image scene understanding: A review
Cai et al. A comparative study of deep learning approaches to rooftop detection in aerial images
Khan et al. A comparative study of graphic symbol recognition methods
Zhou et al. Stratified Object‐Oriented Image Classification Based on Remote Sensing Image Scene Division
CN115019163A (en) City factor identification method based on multi-source big data
Azeez et al. Urban tree classification using discrete-return LiDAR and an object-level local binary pattern algorithm
Liu et al. Second-order texton feature extraction and pattern recognition of building polygon cluster using CNN network
Guo et al. Object discovery in high-resolution remote sensing images: a semantic perspective
Chu et al. Scene classification based on SIFT combined with GIST
Sun et al. The recognition framework of deep kernel learning for enclosed remote sensing objects
Kakooei et al. VHR semantic labeling by random forest classification and fusion of spectral and spatial features on Google Earth Engine
Norelyaqine et al. Architecture of Deep Convolutional Encoder‐Decoder Networks for Building Footprint Semantic Segmentation
Vohra et al. Spatial shape feature descriptors in classification of engineered objects using high spatial resolution remote sensing data

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