CN104851113A - Urban vegetation automatic extraction method of multiple-spatial resolution remote sensing image - Google Patents

Abstract

The present invention discloses an urban vegetation automatic extraction method of multiple-spatial resolution remote sensing images. The method comprises the steps of extracting vegetation initial plaques based on high-spatial resolution multispectral images, calculating visual perception parameters for the high-spatial resolution multispectral images, calculating textural feature parameters for high-spatial resolution panchromatic Images, automatic growing of multi-feature comprehensive vegetation regions, and obtaining a vegetation region distribution map. With adoption of image sectional division treatment, the image processing speed is accelerated. According to normalized difference vegetation index (NDVI) features, automatic selection of maximal mathematical expectation algorithm adaptive dynamic thresholds is adopted. Panchromatic wave bands and multi-spectrum wave bands of the remote sensing images are fully utilized, the NDVI features, visual features and textural features of the panchromatic wave bands and multi-spectrum wave bands are synthesized, vegetation region judgment is performed on the vegetation, and the accuracy is raised.

Description

The urban vegetation extraction method of multiple-spatial resolution remote sensing image

Technical field

The present invention relates to a kind of urban vegetation extraction method of multiple-spatial resolution remote sensing image, belong to remote sensing image field.

Background technology

Urban vegetation has solid carbon, releases foster function, in maintenance urban ecological environment, have vital role.The urban vegetation extracted based on high-resolution remote sensing image covers, and is the important indicator of assessment ecocity, garden-like cities.The method that present stage utilizes high-resolution remote sensing image to extract urban vegetation mainly contains following several:

(1) artificial visual decipher interpretation.The method mainly with artificial interpretation vegetation area, vegetation border of then sketching.The feature that Main Basis vegetation shows in the picture in visual interpretation comprises shape, size, color and tint, shade, position, texture relation etc.Visual interpretation adopts the method for manual work screen digitization, and automaticity is lower, and subjectivity is strong, needs abundant image interpretation practical experience, and has good understanding and cognition to extracted thematic information.But due to various information in remote sensing image can be utilized fully, remain the important method of remote sensing image process.

(2) OO supervised classification method.The method mainly utilizes the color of vegetation area in image, spatial texture characteristic, on OO Iamge Segmentation basis, calculate the various features of each object, (it comprises: spectrum average to comprise spectrum, brightness, mean square deviation, ratio, the spectral signature relevant to field, spectral signature relevant to parent object etc.) shape (comprises area, girth, the length of side, tight ness rating, length breadth ratio and the shape facility relevant to subobject and parent object), texture is (according to the second moment (energy) that gray level co-occurrence matrixes defines, contrast, relevant, entropy, unfavourable balance square) and topological characteristic etc., finally select corresponding feature to carry out object-based classification to extract.

These extracting method above-mentioned all existing defects and problems, comprising:

(1) when extracting vegetation area, needing to select sample in advance, generating priori, but not based on image self-characteristic.Have man-machine interaction in this process, artificial influence factor is many.

(2) automaticity is not high.The yardstick split in assorting process, the threshold value of vegetative breakdown, feature all need man-machine interactively to select, and these pilot process have had a strong impact on the automaticity that vegetation is extracted.In addition, in cutting procedure, entire image adopts same threshold value, does not take into full account difference (as meadow, coniferous forest, the broad-leaf forest etc.) feature difference on image between different vegetation types.

(3) only make use of the multispectral data in high spatial resolution image, do not have to utilize the feature with more high spatial resolution panchromatic wave-band, and in reality, the remotely-sensed data PAN and multi-spectral image as high spatial resolutions such as WorldView, QuickBird bundles.Only make use of the multispectral data in high spatial resolution image, do not have to utilize the feature with more high spatial resolution panchromatic wave-band, and in reality, the remotely-sensed data PAN and multi-spectral image as high spatial resolutions such as WorldView, QuickBird bundles.

Summary of the invention

In order to solve the problem, the present invention is directed to IKONOS, QuickBird, WorldView etc. and there is high-resolution remote sensing image that is panchromatic and multi light spectrum hands, full-automatic extraction urban vegetation region, proposes a kind of urban vegetation extraction method of multiple-spatial resolution remote sensing image.

The invention provides a kind of urban vegetation extraction method of multiple-spatial resolution remote sensing image, mainly comprise the following steps:

Step one, to extract based on the initial patch of the multispectral image vegetation of high spatial resolution;

Step 2, multispectral image computation vision perceptual parameters to high spatial resolution;

Step 3, textural characteristics parameter is calculated to the panchromatic wave-band image of high spatial resolution;

The vegetation area automatic increase of step 4, combination features, obtains vegetation area distribution plan.

Preferably, above-mentioned steps one specifically comprises the following steps:

Step 1.1 calculates NDVI vegetation spectrum index, and computing formula is

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}$

Wherein NIR represents the pixel value of near-infrared band, and R represents the pixel value of red band;

Step 1.2 pair NDVI image carries out piecemeal, every block size 300*300 pixel, when image edge part is less than 300*300 pixel, is as the criterion with the size of reality;

Step 1.3 distinguishes vegetation and non-vegetation to the NDVI block data automatic acquisition segmentation threshold of each 300*300 pixel, obtains initial vegetation area;

Step 1.4 repeats step 1.3, obtains the vegetation prime area of all piecemeals, merges the vegetation area in all piecemeals, obtain initial vegetation area, carries out connectedness mark, obtain the pixel coordinate of each vegetative patch to whole image vegetation area.

Preferably, in above-mentioned steps 1.3, the acquisition of segmentation threshold obtains by the following method:

(1) corresponding several ground mulching type in each piecemeal NDVI, with ω _irepresent, suppose that Gaussian density function distribution is obeyed in the NDVI distribution of each atural object classification, then probability distribution function p (x) that in each piecemeal, NDVI is total can be expressed as:

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left({\mathrm{\ω}}_i\right)p\left(x\right|{\mathrm{\ω}}_i)$

Wherein, other total quantity of n representation class, ω _irepresent i-th atural object classification, average and variance use m respectively _i, represent, its probability density function is expressed as

$p\left(x\right|{\mathrm{\ω}}_i)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\σ}}_i^2}\}$

Wherein, p (ω _i) be the probability density of each atural object classification, its condition met is:

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left({\mathrm{\ω}}_i\right)=1$

Under these conditions, the threshold value solving vegetation and non-vegetation can be converted into the quantity n estimating atural object classification, the ω of each classification _iaverage and variance divide m _i, concrete steps are as follows:

(2) frequency histogram f (x) of NDVI in piecemeal image is added up, the maximal value NDVI of NDVI in search piecemeal image _maxwith minimum value NDVI _min, by NDVI _minto NDVI _maxbe equally divided into 256 parts, statistics NDVI value appears at every a interval pixel quantity, obtains histogram f (x) of NDVI;

(3) single order calculating histogram f (x) of NDVI is reciprocal, is maximum point, minimum point in the position searching NDVI frequency histogram of 0 according to single order inverse;

(4) add up the quantity of maximum point, determine ground mulching categorical measure n, using the NDVI value at maximum value place as the initial mean value m of each ground mulching classification _i, calculate initial σ with adjacent two minimum points _iformula as follows:

σ _i＝(NDVI ₂-NDVI ₁)/4

Wherein, NDVI ₂, NDVI ₁represent the NDVI value that adjacent two minimal value places are corresponding respectively;

$p\left({\mathrm{\ω}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

Wherein, wherein ncount is the number of all pixels of NDVI value between adjacent two minimal values, and n_total is total number of pixels of piecemeal;

(5) based on maximal mathematical expectation (EM) algorithm, loop iteration estimation p (ω is adopted _i), m _i, eM algorithm is by loop iteration, and each iteration forms by asking expectation value and expectation maximization two steps; The former is according to the currency of parameter to be estimated, the expectation value of direct estimation probability density from observation data, and the latter carrys out the estimator of undated parameter by maximizing this expectation, and this two step hockets successively in whole iterative process, until iteration converges; Specific formula for calculation is as follows:

$p^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}{I*I}$

$m_k^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\σ}}_k^2\right)}^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}$

What above-mentioned three formulas were estimated is prior probability, average and standard deviation respectively, in formula, k represents a kth atural object classification, t and t+1 represents estimated value current and that next iteration is used respectively, i, j represents line number and the columns of NDVI image respectively, and X (i, j) represents the NDVI value that in NDVI image, the i-th row j arranges, conditional probability p (X (i, j) | ω _i) calculating and the value of total probability p (X (i, j)) drawn by the respective formula in step (1); As the p (ω of adjacent twice iterative computation _i), m _iand σ _ivalue be less than given threshold epsilon (ε=10 ^-8) time iteration ends;

(6) determining the initial vegetation area of piecemeal, supposing that the NDVI average value of each classification obtained in upper step is by sorting from big to small, m ₁, m ₂, m _n, corresponding variance is σ _i, then initial vegetation distribution plan can be expressed as VI (i, j), and wherein 1 represents vegetation area, and 0 represents territory, nonvegetated area, then

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\≤T\end{array}\right.$

Wherein threshold value T is tried to achieve by following formula

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\σ}}_1& m_1-m_2>1.5({\mathrm{\σ}}_1+{\mathrm{\σ}}_2)\\ m_2+1.5{\mathrm{\σ}}_2& m_1-m_2\≤1.5({\mathrm{\σ}}_1+{\mathrm{\σ}}_2)\end{array}\right..$

Preferably, connectedness mark employing 4 neighborhood in above-mentioned steps 1.4 or 8 neighborhoods.

Preferably, above-mentioned steps two specifically comprises the following steps:

Step 2.1 pair multispectral image has RGB color space transformation to HSI color space, by the result of variations of the brightness of following formulae discovery vegetation, colourity, saturation degree

$\left\{\begin{array}{cc}H=\left\{\begin{array}{cc}\mathrm{\θ},& \begin{array}{cc}\mathrm{if}& G\≥B\end{array}\\ 2\mathrm{\π}-\mathrm{\θ},& \mathrm{otherwise}\end{array}\right.& \mathrm{\θ}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)& \\ I=(R+G+B)/3& \end{array}\right.;$

Step 2.2 computation vision characteristic perception parameter, according to initial vegetative patch as sample, calculates the average μ of the brightness of the vegetation of each patch, colourity, saturation degree _x(H, S, I), variances sigma _x(H, S, I), calculates normalization tone parameter of saturation

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k}& \\ 0,& \mathrm{otherwise}\end{array}\right.,\begin{array}{cc}\mathrm{if}& H_L<H(i,j)<H_U\end{array}$

Wherein, H (i, j) is tone, and (i, j) represents pixel coordinate, μ _kthe average of a kth vegetative patch, the variance of its correspondence is σ _k, H _l, H _ufor lower limit and the upper limit of vegetation tone

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right..$

Preferably, above-mentioned steps three is specially:

For panchromatic wave-band video generation gray level co-occurrence matrixes, calculate 6 parametric textures of gray level co-occurrence matrixes: contrast, diversity, angle energy, entropy, homogeneity, coefficient of autocorrelation, specific formula for calculation is as follows:

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)|i-j|$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}{p(i-j)}^2$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)\log \left(p(i,j)\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)=\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

The selection of textural characteristics computation window mouth and the resolution of panchromatic image and multispectral image have relation, and suppose that the spatial resolution of panchromatic wave-band is n times of multispectral image, then the window that texture calculates is 2n-l, calculates once every n pixel.

Preferably, above-mentioned steps four is specially: be " seed " with each initial vegetative patch, and the pixel adjacent to " seed " patch is attached on seed region by the method adopting region to increase.

Preferably, the criterion that above-mentioned zone increases adopts textural characteristics, visually-perceptible characteristic synthetic to adopt multiple features voting to determine.

Preferably, the criterion of above-mentioned zone growth is as follows:

(1) to each initial vegetative patch, add up brightness (I), normalization tone (H), saturation degree (S), contrast (Contrast), diversity (Dissimilarity), angle energy (Energy), entropy (Entropy) respectively, average, the variance of homogeneity (Homogeneity), coefficient of autocorrelation and above-mentioned characteristic parameter, suppose that its value is μ, σ; To 8 neighborhood territory pixels of initial vegetative patch, according to each feature, judge whether it is vegetation, and discriminant function is defined as follows

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}$

$V(i,j)=\left\{\begin{array}{cc}1,\mathrm{if}& M(i,j)>T_2\\ 0,\mathrm{otherwise}& \end{array}\right.$

Wherein, the characteristic parameter of M (i, j) for calculating, at T ₂represent judgment threshold, i, j represent the ranks number of pixel, and 1 represents vegetation area, and 0 represents territory, nonvegetated area, and then obtain above-mentioned 9 characteristic parameters Potential vegetation areal map corresponding respectively;

(2) multiple features voting determine final vegetation area, this 9 width Potential vegetation areal map is voted, namely on vegetation area figure each pixel when have 6 and above time feature be labeled as 1 time, think final vegetation area.

Relative to prior art, the urban vegetation extraction method of multiple-spatial resolution remote sensing image provided by the invention has following advantage:

(1) present invention employs the vegetation threshold value automatic selecting method based on maximal mathematical expectation, improve Data processing automaticity.Have employed piecemeal treatment technology in addition, view picture image is avoided to adopt same threshold decision, avoid the feature difference that different vegetation types (as meadow, coniferous forest, broad-leaf forest etc.) shows in remotely-sensed data, enhance the adaptivity of vegetation threshold value.

(2) the present invention adopts NDVI to propose initial vegetation area, combines the visually-perceptible feature of high-resolution multi-spectral data, the textural characteristics of panchromatic wave-band on this basis, and the method adopting multiple features to vote determines final vegetation area.This method considers vision, the textural characteristics of high spatial resolution, improves the accuracy that vegetation is extracted.

Accompanying drawing explanation

Fig. 1 is schematic flow sheet of the present invention.

Embodiment

Understand for the ease of those of ordinary skill in the art and implement the present invention, below in conjunction with the drawings and the specific embodiments, the present invention is described in further detail.

This programme has high-resolution remote sensing image that is panchromatic and multi light spectrum hands for IKONOS, QuickBird, WorldView etc., full-automatic extraction urban vegetation region.Concrete step comprises:

Step one, to extract based on the initial patch of the multispectral image vegetation of high spatial resolution.

Calculate NDVI vegetation spectrum index.Specific formula for calculation is shown in (1)

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}---\left(1\right)$

Wherein NIR represents the pixel value of near-infrared band, and R represents the pixel value of red band.

Piecemeal is carried out to NDVI image, every block size 300*300 pixel, when image edge part is less than 300*300 pixel, be as the criterion with the size of reality.

Vegetation and non-vegetation are distinguished to the NDVI block data automatic acquisition segmentation threshold of each 300*300 pixel, obtains initial vegetation area.The method obtaining segmentation threshold is as follows:

Because vegetation has strong absorption at red band, and there is strong reflection near infrared, enhanced this characteristic of vegetation by the NDVI index of the calculating of red, near-infrared band, the difference of vegetation and other types of ground objects can better be reacted.The codomain scope [-1,1] of NDVI, for vegetation, NDVI is generally larger.

In each piecemeal NDVI, corresponding several ground mulching type is (with ω _irepresent), suppose that Gaussian density function distribution is obeyed in the NDVI distribution of each atural object classification, then probability distribution function p (x) that in each piecemeal, NDVI is total can be expressed as

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)p\left(x\right|{\mathrm{\&omega;}}_i)---\left(2\right)$

Wherein, other total quantity of n representation class, ω _irepresent i-th atural object classification, average and variance use m respectively _i, represent, its probability density function is expressed as formula (3), p (ω _i) be the probability density of each atural object classification, its condition met is shown in formula (4).

$p\left(x\right|{\mathrm{\&omega;}}_i)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\&sigma;}}_i^2}\}---\left(3\right)$

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)=1---\left(4\right)$

Under these conditions, the threshold value solving vegetation and non-vegetation can be converted into the quantity n estimating atural object classification, the ω of each classification _iaverage and variance divide m _i, detailed process is as follows:

Frequency histogram f (x) of NDVI in statistics piecemeal image.The maximal value NDVI of NDVI in search piecemeal image _maxwith minimum value NDVI _min, by NDVI _minto NDVI _maxbe equally divided into 256 parts, statistics NDVI value appears at every a interval pixel quantity, obtains histogram f (x) of NDVI.

The single order calculating histogram f (x) of NDVI is reciprocal, is maximum point, minimum point in the position searching NDVI frequency histogram of 0 according to single order inverse.

The quantity of statistics maximum point, determines ground mulching categorical measure n.Using the NDVI value at maximum value place as the initial mean value m of each ground mulching classification _i, calculate initial σ with adjacent two minimum points by formula (5) _i.

σ _i＝(NDVI ₂-NDVI ₁)/4 (5)

NDVI ₂, NDVI ₁represent the NDVI value that adjacent two minimal value places are corresponding respectively.

$p\left({\mathrm{\&omega;}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

Wherein ncount is the number of all pixels of NDVI value between adjacent two minimal values, and n_total is total number of pixels of piecemeal.

Based on maximal mathematical expectation (EM) algorithm, adopt loop iteration estimation p (ω _i), m _i, eM algorithm is by loop iteration, and each iteration forms by asking expectation value and expectation maximization two steps; The former is according to the currency of parameter to be estimated, the expectation value of direct estimation probability density from observation data, and the latter carrys out the estimator of undated parameter by maximizing this expectation, and this two step hockets successively in whole iterative process, until iteration converges; Specific formula for calculation is as follows:

$p^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}{I*I}---\left(5\right)$

$m_k^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}---\left(6\right)$

${\left({\mathrm{\&sigma;}}_k^2\right)}^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}---\left(7\right)$

What three formulas were estimated above is prior probability, average and standard deviation respectively, in formula, k represents a kth atural object classification, t and t+1 represents estimated value current and that next iteration is used respectively, i, j represents line number and the columns of NDVI image respectively, and X (i, j) represents the NDVI value that in NDVI image, the i-th row j arranges, conditional probability p (X (i, j) | ω _i) calculating see (3), the value of total probability p (X (i, j)) is provided by formula (2); As the p (ω of adjacent twice iterative computation _i), m _iand σ _ivalue be less than given threshold epsilon (ε=10 ^-8) time iteration ends.

Determine the initial vegetation area of piecemeal.Suppose that the NDVI average value of each classification obtained in upper step is by sorting from big to small, m ₁, m ₂, m _n, corresponding variance is σ _i, then initial vegetation distribution plan can be expressed as VI (i, j), and wherein 1 represents vegetation area, and 0 represents territory, nonvegetated area.

Then

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\&le;T\end{array}\right.---\left(8\right)$

Wherein threshold value T is shown in formula

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\&sigma;}}_1& m_1-m_2>1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\\ m_2+1.5{\mathrm{\&sigma;}}_2& m_1-m_2\&le;1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\end{array}\right.---\left(9\right)$

Repeat 1.3 steps, obtain the vegetation prime area of all piecemeals.Merge the vegetation area in all piecemeals, obtain initial vegetation area.Then connectedness mark is carried out to whole image vegetation area, obtain the pixel coordinate of each vegetative patch.Connective mark can adopt 4 neighborhoods or 8 neighborhoods.

Step 2, multispectral image computation vision perceptual parameters to high spatial resolution.

2.1 pairs of multispectral images have RGB color space transformation to HSI color space.

The image of red, green, blue wave band corresponding in high-resolution multi-spectral image, calculates the result of variations of vegetation brightness, colourity, saturation degree according to formula (10).

$\left\{\begin{array}{cc}H=\left\{\begin{array}{cc}\mathrm{\&theta;},& \begin{array}{cc}\mathrm{if}& G\&GreaterEqual;B\end{array}\\ 2\mathrm{\&pi;}-\mathrm{\&theta;},& \mathrm{otherwise}\end{array}\right.& \mathrm{\&theta;}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)& \\ I=(R+G+B)/3& \end{array}\right.---\left(10\right)$

2.2 computation vision characteristic perception parameters.

According to initial vegetative patch as sample, calculate the average μ of the brightness of the vegetation of each patch, colourity, saturation degree _x(H, S, I), variances sigma _xstatistical parameter values such as (H, S, I).Calculate normalization tone parameter of saturation

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k}& \\ 0,& \mathrm{otherwise}\end{array}\right.,\begin{array}{cc}\mathrm{if}& H_L<H(i,j)<H_U---\left(11\right)\end{array}$

In formula, (i, j) represents pixel coordinate, μ _kthe average of a kth vegetative patch, the variance of its correspondence is σ k, H _l, H _ufor lower limit and the upper limit of vegetation tone, see formula (11)

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right.---\left(11\right)$

Step 3, textural characteristics parameter is calculated to the panchromatic wave-band image of high spatial resolution

For panchromatic wave-band video generation gray level co-occurrence matrixes, calculate 6 parametric textures of gray level co-occurrence matrixes: contrast (Con), diversity (Dis), angle energy (ASM), entropy (Ent), homogeneity (Hom), coefficient of autocorrelation (Cor), concrete formula is shown in (11-17).The selection of textural characteristics computation window mouth and the resolution of panchromatic image and multispectral image have relation, and the spatial resolution of false color panchromatic wave-band is n times of multispectral image, then the window that texture calculates is 2n-1, calculates once every n pixel.Such as the resolution of quickbird image panchromatic wave-band image is 0.6 1m, and the resolution of multispectral image is 2.4 4m, and the size calculating window during panchromatic wave-band texture is 7 × 7, calculates a textural characteristics every 4 pixels.So just can by the visually-perceptible feature of the textural characteristics of panchromatic image and multispectral image and NDVI corresponding.

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)---\left(11\right)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)|i-j|---\left(12\right)$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}{p(i-j)}^2---\left(13\right)$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)\log \left(p(i,j)\right)---\left(14\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}---\left(15\right)$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)=\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}---\left(16\right)$

The vegetation area automatic increase of step 4, combination features, obtains vegetation area distribution plan

Be " seed " with each initial vegetative patch, the pixel adjacent to " seed " patch is attached on seed region by the method adopting region to increase.The criterion that region increases adopts textural characteristics, visually-perceptible characteristic synthetic to adopt multiple features voting to determine.

The criterion that region increases is as follows

To each initial vegetative patch, add up brightness (I), normalization tone (H), saturation degree (S), contrast (Contrast), diversity (Dissimilarity), angle energy (Energy), entropy (Entropy) respectively, homogeneity (Homogeneity), coefficient of autocorrelation, the average of these 9 characteristic parameters, variance, suppose that its value is μ, σ.To 8 neighborhood territory pixels of initial vegetative patch, according to each feature, judge whether it is vegetation, and discriminant function is defined as follows

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}---\left(17\right)$

$V(i,j)=\left\{\begin{array}{cc}1,\mathrm{if}& M(i,j)>T_2\\ 0,\mathrm{otherwise}& \end{array}\right.---\left(18\right)$

Wherein, M (i, j) refers to the characteristic parameter calculated, be exactly the visual signature parameter of 3 described in (1), 6 textural characteristics parameters, have 9, i.e. brightness, normalization tone, saturation degree, contrast, diversity, angle energy, entropy, homogeneity, coefficient of autocorrelation.Herein with any one variable that M replaces wherein.Represent judgment threshold at T2, i, j represent the ranks number of pixel, and 1 represents vegetation area, and 0 represents territory, nonvegetated area.

Like this, can brightness, normalization tone, saturation degree, contrast (Contrast), diversity (Dissimilarity), angle energy (Energy), entropy (Entropy), homogeneity (Homogeneity), coefficient of autocorrelation, these 9 characteristic parameters can obtain a width Potential vegetation areal map respectively.

Final vegetation area is determined in multiple features voting.Voted by this 9 width Potential vegetation areal map, namely on vegetation area figure, each pixel, when having 6 and above feature is labeled as 1, thinks final vegetation area.

The urban vegetation extraction method of multiple-spatial resolution remote sensing image provided by the invention, have employed image block dividing processing, accelerates image processing speed.According to NDVI feature, adopt the automatic selection of the adaptive dynamic threshold of maximal mathematical expectation algorithm.The panchromatic wave-band of remote sensing image, multi light spectrum hands are made full use of, to the NDVI feature extracted from panchromatic wave-band and multi light spectrum hands, visual signature, textural characteristics comprehensively, vegetation area judgement is carried out to vegetation, improves accuracy.Method provided by the invention also has certain applicability, as the NDVI vegetation index in this method can change EVI index, SR index etc. into.

The above, only that specific embodiment of the invention case is described, and be not used to limit of the present invention can practical range, such as all equivalences that those skilled in the art complete under the spirit do not departed from indicated by the present invention and principle change or modify, and must be covered by the scope of the claims in the present invention.

Claims (9)

1. a urban vegetation extraction method for multiple-spatial resolution remote sensing image, is characterized in that described method mainly comprises the following steps:

Step one, to extract based on the initial patch of the multispectral image vegetation of high spatial resolution;

Step 2, multispectral image computation vision perceptual parameters to high spatial resolution;

Step 3, textural characteristics parameter is calculated to the panchromatic wave-band image of high spatial resolution;

The vegetation area automatic increase of step 4, combination features, obtains vegetation area distribution plan.

2. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 1, is characterized in that: described step one specifically comprises the following steps:

Step 1.1 calculates NDVI vegetation spectrum index, and computing formula is

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}$

Wherein NIR represents the pixel value of near-infrared band, and R represents the pixel value of red band;

Step 1.2 pair NDVI image carries out piecemeal, every block size 300*300 pixel, when image edge part is less than 300*300 pixel, is as the criterion with the size of reality;

Step 1.3 distinguishes vegetation and non-vegetation to the NDVI block data automatic acquisition segmentation threshold of each 300*300 pixel, obtains initial vegetation area;

Step 1.4 repeats step 1.3, obtains the vegetation prime area of all piecemeals, merges the vegetation area in all piecemeals, obtain initial vegetation area, carries out connectedness mark, obtain the pixel coordinate of each vegetative patch to whole image vegetation area.

3. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 2, is characterized in that: in described step 1.3, the acquisition of segmentation threshold obtains by the following method:

(1) corresponding several ground mulching type in each piecemeal NDVI, with ω _irepresent, suppose that Gaussian density function distribution is obeyed in the NDVI distribution of each atural object classification, then probability distribution function p (x) that in each piecemeal, NDVI is total can be expressed as:

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)p\left(x\right|{\mathrm{\&omega;}}_i)$

Wherein, other total quantity of n representation class, ω _irepresent i-th atural object classification, average and variance use m respectively _i, represent, its probability density function is expressed as

$p\left(x\right|{\mathrm{\&omega;}}_i)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\&sigma;}}_i^2}\}$

Wherein, p (ω _i) be the probability density of each atural object classification, its condition met is:

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)=1$

Under these conditions, the threshold value solving vegetation and non-vegetation can be converted into the quantity n estimating atural object classification, the ω of each classification _iaverage and variance divide m _i, concrete steps are as follows:

(2) frequency histogram f (x) of NDVI in piecemeal image is added up, the maximal value NDVI of NDVI in search piecemeal image _maxwith minimum value NDVI _min, by NDVI _minto NDVI _maxbe equally divided into 256 parts, statistics NDVI value appears at every a interval pixel quantity, obtains histogram f (x) of NDVI;

(3) single order calculating histogram f (x) of NDVI is reciprocal, is maximum point, minimum point in the position searching NDVI frequency histogram of 0 according to single order inverse;

(4) add up the quantity of maximum point, determine ground mulching categorical measure n, using the NDVI value at maximum value place as the initial mean value m of each ground mulching classification _i, calculate initial σ with adjacent two minimum points _iformula as follows:

σ _i＝(NDVI ₂-NDVI ₁)/4

Wherein, NDVI ₂, NDVI ₁represent the NDVI value that adjacent two minimal value places are corresponding respectively;

$p\left({\mathrm{\&omega;}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

Wherein, wherein ncount is the number of all pixels of NDVI value between adjacent two minimal values, and n_total is total number of pixels of piecemeal;

(5) based on maximal mathematical expectation (EM) algorithm, loop iteration estimation p (ω is adopted _i), m _i, by loop iteration, each iteration forms by asking expectation value and expectation maximization two steps; Ask the step of expectation value according to the currency of parameter to be estimated, the expectation value of direct estimation probability density from observation data, the step of expectation maximization carrys out the estimator of undated parameter by maximizing this expectation, this two step hockets successively in whole iterative process, until iteration converges; Specific formula for calculation is as follows:

$p^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}{I*I}$

$m_k^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\&sigma;}}_k^2\right)}^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}$

What above-mentioned three formulas were estimated is prior probability, average and standard deviation respectively, in formula, k represents a kth atural object classification, t and t+1 represents estimated value current and that next iteration is used respectively, i, j represents line number and the columns of NDVI image respectively, and X (i, j) represents the NDVI value that in NDVI image, the i-th row j arranges, conditional probability p (X (i, j) | ω _i) calculating and the value of total probability p (X (i, j)) drawn by the respective formula in step (1); As the p (ω of adjacent twice iterative computation _i), m _iand σ _ivalue be less than given threshold epsilon (ε=10 ^-8) time iteration ends;

(6) determining the initial vegetation area of piecemeal, supposing that the NDVI average value of each classification obtained in upper step is by sorting from big to small, m ₁, m ₂, m _n, corresponding variance is σ _i, then initial vegetation distribution plan can be expressed as VI (i, j), and wherein 1 represents vegetation area, and 0 represents territory, nonvegetated area, then

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\&le;T\end{array}\right.$

Wherein threshold value T is tried to achieve by following formula

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\&sigma;}}_1& m_1-m_2>1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\\ m_2+1.5{\mathrm{\&sigma;}}_2& m_1-m_2\&le;1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\end{array}\right..$

4. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 2, is characterized in that: connectedness mark employing 4 neighborhood in described step 1.4 or 8 neighborhoods.

5. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 1, is characterized in that: described step 2 specifically comprises the following steps:

Step 2.1 pair multispectral image has RGB color space transformation to HSI color space, by the result of variations of the brightness of following formulae discovery vegetation, colourity, saturation degree

$\left\{\begin{array}{c}H=\left\{\begin{array}{cc}\mathrm{\&theta;},& \mathrm{ifG}\&GreaterEqual;B\\ 2\mathrm{\&pi;}-\mathrm{\&theta;},& \mathrm{otherwise}\end{array}\right.\mathrm{\&theta;}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)\\ I=(R+G+B)/3\end{array}\right.;$

Step 2.2 computation vision characteristic perception parameter, according to initial vegetative patch as sample, calculates the average μ of the brightness of the vegetation of each patch, colourity, saturation degree _x(H, S, I), variances sigma _x(H, S, I), calculates normalization tone parameter of saturation

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k},& \mathrm{if}{H}_L<H(i,j)<H_U\\ 0,& \mathrm{otherwise}\end{array}\right.$

Wherein, H (i, j) is tone, and (i, j) represents pixel coordinate, μ _kthe average of a kth vegetative patch, the variance of its correspondence is σ _k, H _l, H _ufor lower limit and the upper limit of vegetation tone

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right..$

6. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 1, is characterized in that described step 3 is specially:

For panchromatic wave-band video generation gray level co-occurrence matrixes, calculate 6 parametric textures of gray level co-occurrence matrixes: contrast, diversity, angle energy, entropy, homogeneity, coefficient of autocorrelation, specific formula for calculation is as follows:

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)|i-j|$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}p{(i,j)}^2$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)\log \left(p(i,j)\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

The selection of textural characteristics computation window mouth and the resolution of panchromatic image and multispectral image have relation, and suppose that the spatial resolution of panchromatic wave-band is n times of multispectral image, then the window that texture calculates is 2n-1, calculates once every n pixel.

7. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 1, it is characterized in that described step 4 is specially: with each initial vegetative patch for " seed ", the pixel adjacent to " seed " patch is attached on seed region by the method adopting region to increase.

8. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 7, is characterized in that the criterion that described region increases adopts textural characteristics, visually-perceptible characteristic synthetic to adopt multiple features voting to determine.

9. the urban vegetation extraction method of multiple-spatial resolution remote sensing image according to claim 8, is characterized in that the criterion that described region increases is as follows:

(1) to each initial vegetative patch, add up brightness, normalization tone, saturation degree, contrast, diversity, angle energy, entropy respectively, average, the variance of homogeneity, coefficient of autocorrelation and above-mentioned characteristic parameter, suppose that its value is μ, σ; To 8 neighborhood territory pixels of initial vegetative patch, according to each feature, judge whether it is vegetation, and discriminant function is defined as follows

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}$

$V(i,j)=\left\{\begin{array}{cc}1,& \mathrm{ifM}(i,j)>T_2\\ 0,& \mathrm{otherwise}\end{array}\right.$

Wherein, the characteristic parameter of M (i, j) for calculating, at T ₂represent judgment threshold, i, j represent the ranks number of pixel, and 1 represents vegetation area, and 0 represents territory, nonvegetated area, and then obtain above-mentioned 9 characteristic parameters Potential vegetation areal map corresponding respectively;

(2) multiple features voting determine final vegetation area, this 9 width Potential vegetation areal map is voted, namely on vegetation area figure each pixel when have 6 and above time feature be labeled as 1 time, think final vegetation area.