CN104881867A

CN104881867A - Method for evaluating quality of remote sensing image based on character distribution

Info

Publication number: CN104881867A
Application number: CN201510243039.0A
Authority: CN
Inventors: 田岩; 张慧敏; 阮崇武; 许毅平
Original assignee: Huazhong University of Science and Technology
Current assignee: Huazhong University of Science and Technology
Priority date: 2015-05-13
Filing date: 2015-05-13
Publication date: 2015-09-02

Abstract

The invention discloses a remote sensing image quality evaluation method based on feature distribution, comprising the following steps: sampling each category in the input remote sensing image, extracting grayscale and features; calculating the mean value and variance of each category in the feature samples, as The initial value of the EM estimation (GMM estimation) of the Gaussian model is estimated by the Gaussian mixture model GMM; according to the variance and weight of each category feature obtained by the GMM estimation, the quality evaluation model based on the kappa coefficient is constructed to calculate the quality of the remote sensing image The kappa coefficient. The present invention makes full use of the image features of the remote sensing image when evaluating the quality evaluation of the remote sensing image, and uses the method of image classification to verify the evaluation method, so that the image quality can be effectively evaluated objectively. The fixed mean Gaussian mixture model GMM proposed by the invention has better convergence.

Description

A Quality Evaluation Method of Remote Sensing Image Based on Feature Distribution

技术领域 technical field

本发明属于遥感图像质量评价应用领域，具体涉及遥感图像纹理特征的计算方法、纹理分布的估计方法以及基于特征分布的遥感图像分类精度估计模型。 The invention belongs to the application field of quality evaluation of remote sensing images, and in particular relates to a calculation method of texture features of remote sensing images, an estimation method of texture distribution and a classification accuracy estimation model of remote sensing images based on feature distribution.

背景技术 Background technique

目前，各国遥感探测的发展方兴未艾，我国也已经启动高分辨率对地观测的计划。随着该计划的实施，我国将获取大量自主产权的遥感数据。为了利用这些数据大幅度提高我国自主对地观测信息的能力，包括图像分类等在内的图像解译是不可或缺的一个技术环节。 At present, the development of remote sensing detection in various countries is in the ascendant, and my country has also launched a high-resolution earth observation plan. With the implementation of the plan, my country will obtain a large amount of remote sensing data with independent property rights. In order to use these data to greatly improve my country's independent earth observation information capabilities, image interpretation, including image classification, is an indispensable technical link.

目前的图像质量评价方法主要从主观评价的角度进行的，就机器图像分类而言，这些方法无法反映机器分类的可分性。图像的主观评价较高，对人工判读会有极大的帮助，但是对于机器分类而言，无法度量用于遥感图像地物分类所用特征的可分性和稳定性。 The current image quality evaluation methods are mainly carried out from the perspective of subjective evaluation. As far as machine image classification is concerned, these methods cannot reflect the separability of machine classification. The subjective evaluation of images is high, which will be of great help to manual interpretation, but for machine classification, it is impossible to measure the separability and stability of features used for remote sensing image classification.

本发明提出基于特征分布的图像质量评价方法(其中特征包括灰度，能量，对比度，逆差矩，熵，相关性)，从遥感图像地物分类的角度研究卫星全色图像质量评估方法，这与传统的侧重于图像逼真度的质量评估具有迥然不同的含义，具有重要的实际应用价值。 The present invention proposes an image quality evaluation method based on feature distribution (wherein the feature includes gray scale, energy, contrast, inverse moment, entropy, correlation), and studies the satellite panchromatic image quality evaluation method from the perspective of remote sensing image ground object classification, which is similar to The traditional quality assessment that focuses on image fidelity has very different meanings and has important practical application value.

发明内容 Contents of the invention

针对遥感图像在机器分类方面是否具有良好的可分性和稳定性，即如何去评价其在机器分类方面的质量好坏，本发明提出了一种基于特征分布的图像质量评价方法，其具体步骤如下： Aiming at whether the remote sensing image has good separability and stability in machine classification, that is, how to evaluate its quality in machine classification, the present invention proposes an image quality evaluation method based on feature distribution, and its specific steps as follows:

S1：对图像中各个地物类别进行采样，提取灰度及纹理特征：设有遥感图像I，选取图像中各个类别的部分区域作为样本(根据)，以样本中的每个像元为中心点，开窗大小为K×K(K取5～21)，在此图像块的中统计其中心点的灰度共生矩阵，计算出纹理特征，过程如下： S1: Sampling each feature category in the image, and extracting grayscale and texture features: a remote sensing image I is set, and some areas of each category in the image are selected as samples (according to), and each pixel in the sample is used as the center point , the size of the window is K×K (K takes 5 to 21), count the gray level co-occurrence matrix of the center point in this image block, and calculate the texture feature, the process is as follows:

(1)灰度共生矩阵统计：取图像块(K×K)中任意一点〔x，y〕及偏离它的另一点〔x+a，y+b〕，设该点对的灰度值为〔g₁，g₂〕。令点〔x，y〕在整个画面上移动，则会得到各种〔g₁，g₂〕值，设灰度值的级数为k，则(g₁，g₂)的组合共有k²种。对于整个画面，统计出每一种〔g₁，g₂〕值出现的次数，然后排列成一个方阵，再用〔g₁，g₂〕出现的总次数将它们归一化为出现的概率P〔g₁，g₂〕，这样的方阵称为灰度共生矩阵。距离差分值〔a，b〕取不同的数值组合，得到不同数值组合下的联合概率矩阵。 (1) Gray level co-occurrence matrix statistics: Take any point [x, y] in the image block (K×K) and another point [x+a, y+b] that deviates from it, and set the gray value of the point pair as [g ₁ , g ₂ ]. Let the point [x, y] move on the whole screen, and various [g ₁ , g ₂ ] values will be obtained. If the number of levels of the gray value is k, then the combination of (g ₁ , g ₂ ) has a total of k ² kind. For the whole picture, count the number of occurrences of each [g ₁ , g ₂ ] value, and then arrange them into a square matrix, and then use the total number of occurrences of [g ₁ , g ₂ ] to normalize them to the probability of occurrence P[g ₁ , g ₂ ], such a square matrix is called a gray level co-occurrence matrix. The distance difference value [a, b] takes different numerical combinations to obtain the joint probability matrix under different numerical combinations.

(2)纹理特征计算：常用的描述纹理的特征有能量(Energy)、对比度(Contrast)、逆差矩(Inverse Difference Moment)、熵(Entropy)和相关(Correlation)，其中，L为灰度级，为归一化的灰度共生矩阵，能量用来衡量灰度分布的均匀性，其计算公式为： (2) Texture feature calculation: Commonly used features describing texture include Energy, Contrast, Inverse Difference Moment, Entropy, and Correlation, where L is the gray level, is a normalized gray-scale co-occurrence matrix, and energy is used to measure the uniformity of gray-scale distribution, and its calculation formula is:

${f f}_{11} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} \overset{^^}{p p} {((i i,, j j))}^{22}$

对比度反映了局部纹理的变化剧烈程度，其计算公式为： Contrast reflects the intensity of local texture changes, and its calculation formula is:

${f f}_{22} = = {Σ Σ}_{| | i i - - j j | | = = 00}^{L L - - 11} {((i i - - j j))}^{22} {{{Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} \overset{^^}{p p} ((i i,, j j))}}$

逆差距反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀，其计算公式为： The inverse gap reflects the homogeneity of the image texture and measures the local variation of the image texture. A large value indicates that there is a lack of change between different regions of the image texture, and the locality is very uniform. The calculation formula is:

${f f}_{33} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} \frac{p p ((i i,, j j))}{11 + + {((i i - - j j))}^{22}}$

熵是图像中所具有的信息量的度量，纹理的复杂度越高就意味着图像信息量越大，其熵越大，公式如下： Entropy is a measure of the amount of information in an image. The higher the complexity of the texture, the greater the amount of image information and the greater the entropy. The formula is as follows:

${f f}_{44} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} - - \overset{^^}{p p} ((i i,, j j)) ln ln ((\overset{^^}{p p} ((i i,, j j))))$

相关性指标度量灰度共生矩阵里的各值在行与列上的相似程度。因此，相关性的大小反映了图像中局部灰度分布的相关性。当矩阵元素值均匀相等时，相关就大；相反，如果矩阵像元值相差很大则相关值小。其中μ_x，μ_y为分别为灰度共生矩阵行方向和列方向上的均值，δ_x，δ_y分别为灰度共生矩阵行方向和列方向上的方差，其计算公式如下： The correlation index measures the similarity of each value in the gray level co-occurrence matrix on the row and column. Therefore, the magnitude of the correlation reflects the correlation of the local gray distribution in the image. When the matrix element values are evenly equal, the correlation is large; on the contrary, if the matrix cell values differ greatly, the correlation value is small. Among them, μ _x , μ _y are the mean values in the row direction and column direction of the gray level co-occurrence matrix respectively, δ _x , δ _y are the variances in the row direction and column direction of the gray level co-occurrence matrix respectively, and the calculation formula is as follows:

${f f}_{55} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} \frac{((i i * * j j)) \overset{^^}{p p} ((i i,, j j)) - - {μ μ}_{x x} {μ μ}_{y the y}}{{δ δ}_{x x} {δ δ}_{y the y}}$

式中： In the formula:

${μ μ}_{x x} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} i* i* \overset{^^}{p p} ((i i,, j j))$

${μ μ}_{y the y} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} i i * * \overset{^^}{p p} ((i i,, j j))$

${δ δ}_{x x} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} {((i i - - {μ μ}_{x x}))}^{22} * * \overset{^^}{p p} ((i i,, j j))$

${δ δ}_{y the y} = = {Σ Σ}_{i i = = 00}^{L L - - 11} {Σ Σ}_{j j = = 00}^{L L - - 11} {((j j - - {μ μ}_{y the y}))}^{22} * * \overset{^^}{p p} ((i i,, j j))$

S2：均值约束的GMM(高斯混合模型)参数估计：遥感影像中各类地物由于某些地物的特征分布极其相似，特征的概率密度曲线混叠较严重。所以本文基于遥感图像的特点提出添加均值约束的GMM算法对提取的样本特征进行计算。其计算过程如下： S2: Mean-constrained GMM (Gaussian Mixture Model) parameter estimation: In remote sensing images, the feature distribution of various ground features is very similar, and the probability density curves of the features are seriously aliased. Therefore, based on the characteristics of remote sensing images, this paper proposes the GMM algorithm with mean value constraints to calculate the extracted sample features. Its calculation process is as follows:

(1)设当前图像中当前特征的每类均值为m₁，...m_C，则令初始参数值为θ⁰＝{a₁，...，a_C，m₁，...m_C，δ₁，...δ_C}；其中a₁，...，a_C为每类的类别权重，该权重由用户根据图像地物先验分布确定，如无先验知识，就默认为1/C，即平均分布，满足δ₁，...δ_C为特征值方差，C为类别数； (1) Assuming that the mean value of each class of the current feature in the current image is m ₁ ,...m _C , then the initial parameter value is θ ⁰ ={a ₁ ,...,a _C ,m ₁ ,...m _C , δ ₁ , ... δ _C }; where a ₁ , ..., a _C is the category weight of each category, which is determined by the user according to the prior distribution of image features. If there is no prior knowledge, the default is is 1/C, that is, the average distribution, satisfying δ ₁ ,...δ _C is the variance of eigenvalues, and C is the number of categories;

(2)由θ⁰迭代t次得到θ^t，t为迭代次数，利用参数值θ^t计算当前像素属于第j类的后验概率β_j(x)，j＝1，2，...C，计算式如下： (2) Obtain θ ^t from θ ⁰ iteration t times, t is the number of iterations, use the parameter value θ ^t to calculate the posterior probability β _j (x) that the current pixel belongs to the jth class, j=1, 2,...C , the calculation formula is as follows:

${β β}_{j j} ((x x)) = = \frac{{a a}_{j j} g g ((x x,, {μ μ}_{j j},, {δ δ}_{j j}))}{{Σ Σ}_{k k = = 11}^{C C} g g ((x x,, {μ μ}_{k k},, {δ δ}_{k k}))}$

其中g(x，μ_j，δ_j)为均值为μ_j，方差为δ_j的高斯分布的概率密度函数，即 where g(x, μ _j , δ _j ) is the probability density function of Gaussian distribution with mean value μ _j and variance δ _j , namely

$g g ((x x,, {μ μ}_{j j},, {δ δ}_{j j})) = = \frac{11}{{((22 π π))}^{11 / / 22} {| | {δ δ}_{j j} | |}^{11 / / 22}} exp exp ((- - [[\frac{11}{22} {((x x - - {μ μ}_{j j}))}^{T T} {δ δ}_{j j}^{- - 11} ((x x - - {μ μ}_{j j}))]]));;$

(3)固定均值，以θ⁰＝{a₁，...，a_C，m₁，...m_C，δ₁，...δ_C}为初值，迭代计算各参数，包括地物类别权重、协方差，计算公式为： (3) Fixed mean value, with θ ⁰ ={a ₁ ,...,a _C ,m ₁ ,...m _C ,δ ₁ ,...δ _C } as the initial value, iteratively calculate each parameter, including ground The weight and covariance of the object category, the calculation formula is:

${β β}_{j j} ((x x)) = = \frac{{a a}_{j j} g g ((x x,, {μ μ}_{j j},, {δ δ}_{j j}))}{{Σ Σ}_{k k = = 11}^{C C} g g ((x x,, {μ μ}_{k k},, {δ δ}_{k k}))},, j j = = 1,2 1,2 . . . . C C,,$

${a a}_{j j}^{new new} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i})),, j j = = 1,2 1,2 . . . . C C,,$

${δ δ}_{j j}^{new new} = = \frac{{Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i})) (({x x}_{i i} - - {μ μ}_{j j})) {(({x x}_{i i} - - {μ μ}_{j j}))}^{T T}}{{Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i}))},, j j = = 1,2 1,2,, . . . . . . C C,,$

N为图像中像素个数； N is the number of pixels in the image;

重复(2)(3)，若||θ^new-θ||＜ζ，ζ为误差值(根据精度需求选取，本发明取小于10^-5)，迭代停止，则得到固定当前均值为m₁，...m_c时的权重和方差估计结果。 Repeat (2) (3), if ||θ ^new -θ||<ζ, ζ is the error value (selected according to the accuracy requirement, the present invention is less than 10 ^-5 ), the iteration stops, and the fixed current mean value is m ₁ , ... m _c weight and variance estimation results.

S3：构建基于kappa系数的质量评价模型： S3: Build a quality evaluation model based on the kappa coefficient:

kappa系数是一种评价整体精度的度量指标，它可以作为样本的分类结果与真实的地物类型的一致性检验。其公式如下： The kappa coefficient is a measure to evaluate the overall accuracy, and it can be used as a consistency test between the classification result of the sample and the real object type. Its formula is as follows:

$Kappa = \frac{N Σ_{i = 1}^{C} n_{ii} - Σ_{i = 1}^{C} N_{i +} N_{+ i}}{N^{2} - Σ_{i = 1}^{C} N_{i +} N_{+ i}} = \frac{\frac{Σ_{i = 1}^{C} n_{ii}}{N} - Σ_{i = 1}^{C} \frac{N_{i +} N_{+ i}}{N^{2}}}{1 - Σ_{i = 1}^{C} \frac{N_{i +} N_{+ i}}{N^{2}}}$ (式1) $Kappa = \frac{N Σ_{i = 1}^{C} {no}_{i} - Σ_{i = 1}^{C} N_{i +} N_{+ i}}{N^{2} - Σ_{i = 1}^{C} N_{i +} N_{+ i}} = \frac{\frac{Σ_{i = 1}^{C} {no}_{i}}{N} - Σ_{i = 1}^{C} \frac{N_{i +} N_{+ i}}{N^{2}}}{1 - Σ_{i = 1}^{C} \frac{N_{i +} N_{+ i}}{N^{2}}}$ (Formula 1)

令P_k是通过S2计算出的类别权重，作为先验概率，且有P_kk表示真实情况下属于第k类，也被正确归为第k类的像素比例；P_kadd代表实际不属于第k类，而被误分为第k类的像素比例。利用P_k、P_kk以及P_kadd，可将Kappa系数改写成为与分类结果直接相关的、各部分意义明确的新形式，如下所示： Let P _k be the category weights calculated by S2 as prior probabilities, and have P _kk represents the proportion of pixels that actually belong to the kth class and are correctly classified as the kth class; P _kadd represents the proportion of pixels that do not actually belong to the kth class but are mistakenly classified as the kth class. Using P _k , P _kk and P _kadd , the Kappa coefficient can be rewritten into a new form that is directly related to the classification result and has clear meanings for each part, as follows:

$kappa = 1 - \frac{Σ_{k = 1}^{C} P_{kadd}}{Σ_{k = 1}^{C} [P_{k} (P_{kk} + P_{kadd})]}$ (式2) $kappa = 1 - \frac{Σ_{k = 1}^{C} P_{kadd}}{Σ_{k = 1}^{C} [P_{k} (P_{kk} + P_{kadd})]}$ (Formula 2)

由上式可知，Kappa系数的预估依赖于P_k、P_kk以及P_kadd的求解。 It can be seen from the above formula that the estimation of Kappa coefficient depends on the solution of P _k , P _kk and P _kadd .

由于高斯分布是常见的分布形态，这里假设各类地物特征均服从多维高斯分布。对于一维特征，设μ₁，μ₂，...μ_C(μ₁＜μ₂＜…＜μ_C分别为各类的均值；σ₁，σ₂...σ_C为对应各类的标准差；P₁，P₂，..P_C为对应各类的先验概率；函数Φ(x)表示标准正态分布的累积分布函数；h_j(1≤j＜C)表示第j类与第j+1类间的分类决策面，则任选k＝1，2..C，P_kk和P_kadd可表示成： Since the Gaussian distribution is a common distribution form, it is assumed that all kinds of ground features obey the multidimensional Gaussian distribution. For one-dimensional features, set μ ₁ , μ ₂ , ...μ _C (μ ₁ <μ ₂ <...<μ _C are the mean values of each type; σ ₁ , σ ₂ ...σ _C are the corresponding values of each type Standard deviation; P ₁ , P ₂ , ..P _C is the prior probability corresponding to each category; function Φ(x) represents the cumulative distribution function of the standard normal distribution; h _j (1≤j<C) represents the jth class The classification decision surface between the j+1th category, then optional k=1, 2..C, P _kk and P _kadd can be expressed as:

$P_{kk} = \{\begin{matrix} P {x \leq h_{k} | x &Element; C_{k}}, k = 1 \\ P {h_{k - 1} < x \leq h_{k} | x &Element; C_{k}}, 2 \leq k < C \\ P {x > h_{k} | x &Element; C_{k}}, k = C \end{matrix}$ (式3) $P_{kk} = \{\begin{matrix} P {x \leq h_{k} | x &Element; C_{k}}, k = 1 \\ P {h_{k - 1} < x \leq h_{k} | x &Element; C_{k}}, 2 \leq k < C \\ P {x > h_{k} | x &Element; C_{k}}, k = C \end{matrix}$ (Formula 3)

$P_{kadd} = Σ_{j = 1, j! = k}^{C} P_{j} * P {h_{k - 1} \leq x < h_{k} | x &Element; C_{k}}$ (式4) $P_{kadd} = Σ_{j = 1, j! = k}^{C} P_{j} * P {h_{k - 1} \leq x < h_{k} | x &Element; C_{k}}$ (Formula 4)

因此，基于不同的决策方法得到不同的分界面h_j(1≤j＜C)，P_kk和P_kadd的值随之变化。本发明基于最小距离准则。 Therefore, different interface h _j (1≤j<C) is obtained based on different decision-making methods, and the values of P _kk and P _kadd change accordingly. The invention is based on the minimum distance criterion.

最小距离准则是一种常见的分类准则。它利用不同类别间地物属性相差较大，而同一类别内地物属性相差较小的原则，根据像素或对象到各类中心的距离进行分类。若像素或对象到某一类的类中心距离最小，则该像素或者对象则被标记为此类。例如，对于灰度图像而言，若x_i为任意一个像素，μ_j(j＝1，2，..C)为第j类的中心，且则x_i被分类第t类。 The minimum distance criterion is a common classification criterion. It makes use of the principle that the attributes of ground objects in different categories differ greatly, while the attributes of ground objects in the same category have small differences, and classify according to the distance between pixels or objects and the centers of various types. If the pixel or object has the smallest distance to the class center of a certain class, the pixel or object is marked as this class. For example, for a grayscale image, if x _i is any pixel, μ _j (j=1, 2, ..C) is the center of class j, and Then x _i is classified into the tth class.

基于上述正态分布的假设，根据最小距离分类准则下，第j类与第j+1类的分界面可表示为： Based on the above assumption of normal distribution, according to the minimum distance classification criterion, the interface between class j and class j+1 can be expressed as:

${h h}_{j j} = = \frac{{μ μ}_{j j} + + {μ μ}_{j j} + + 11}{22} ((11 \leq \leq j j < < C C))$

则对不同的k，利用概率理论，(式3)、(式4)可计算如下： Then for different k, using probability theory, (Equation 3) and (Equation 4) can be calculated as follows:

①若k＝1， ①If k=1,

$\{\begin{matrix} P_{kk} = P_{k} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \end{matrix}$ (式5) $\{\begin{matrix} P_{kk} = P_{k} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \end{matrix}$ (Formula 5)

②若1＜k＜C， ②If 1<k<C,

$\{\begin{matrix} P_{kk} = P_{k} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix}$ (式6) $\{\begin{matrix} P_{kk} = P_{k} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix}$ (Formula 6)

③若k＝C， ③If k=C,

$\{\begin{matrix} P_{kk} = P_{k} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 1}^{C - 1} P_{j} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix}$ (式7) $\{\begin{matrix} P_{kk} = P_{k} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 1}^{C - 1} P_{j} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix}$ (Formula 7)

(式2)即可作为面向图像分类的图像质量评价模型，(式5)～(式7)为各类均服从正态分布下的图像质量评价模型具体计算方法。当分类准则和分类特征确定时，一幅图像的分类性能完全可以根据(式5)计算Kappa系数进行推估. (Equation 2) can be used as an image quality evaluation model for image classification, and (Equation 5) to (Equation 7) are specific calculation methods for image quality evaluation models under normal distribution. When the classification criteria and classification features are determined, the classification performance of an image can be estimated by calculating the Kappa coefficient according to (Formula 5).

附图说明 Description of drawings

图1：本发明方法的基本流程图。 Figure 1: Basic flow diagram of the method of the present invention.

图2：实验数据 Figure 2: Experimental data

图3：德国KOB地区的数据，从中截取田地、林地和城区三种地物类型，构建200*600大小的图片 Figure 3: The data of the KOB region in Germany, from which three types of ground objects, namely fields, woodlands and urban areas, were intercepted to construct a 200*600 size picture

图4是本发明中基于均值约束GMM估计结果与GMM估计结果(灰度特征)比较图； Fig. 4 is a comparison diagram based on the mean value constraint GMM estimation result and the GMM estimation result (gray feature) in the present invention;

图5是本发明与基于最小距离分类的真实kappa系数计算结果比较图。 Fig. 5 is a comparison chart of the present invention and the real kappa coefficient calculation results based on minimum distance classification.

具体实施方式 Detailed ways

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及示例性实施例，对本发明进行进一步详细说明。如图1所示，本发明的具体过程为： In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be further described in detail below in conjunction with the accompanying drawings and exemplary embodiments. As shown in Figure 1, the concrete process of the present invention is:

(1)提取纹理特征：对输入的大小为600×200的遥感图像进行基于灰度共生矩阵的纹理提取，本实施例开窗大小为取21×21。 (1) Extract texture features: perform texture extraction based on the gray-level co-occurrence matrix for the input remote sensing image with a size of 600×200, and the window size in this embodiment is 21×21.

(2)获得均值初值：对原始遥感图像，以图像中地物类别为标准(本实施例图像类别C取3)，选取各个类别图像采样，统计出灰度共生矩阵，计算各个特征的均值和方差，选择灰度特征为本实施例的特征，作为GMM估计均值初值。 (2) Obtaining the initial value of the mean value: for the original remote sensing image, take the object category in the image as the standard (the image category C in this embodiment is taken as 3), select image samples of each category, and calculate the gray level co-occurrence matrix, and calculate the mean value of each feature and variance, the gray feature is selected as the feature of this embodiment, as the initial value of the GMM estimated mean.

(3)固定均值的GMM参数估计：初始参数值θ⁰＝{a₁，a₂，a₃，m₁，m₂，m₃，δ₁，δ₂，δ₃}＝{0.33，0.33，0.33，65.6，14.8，117.75，20.4，7.34，61.39}，利用参数值θ^t计算后验概率β_j(x)，j＝1，2，..C，t为迭代次数。固定均值，计算各参数(权重、协方差)： (3) GMM parameter estimation with fixed mean: initial parameter value θ ⁰ = {a ₁ , a ₂ , a ₃ , m ₁ , m ₂ , m ₃ , δ ₁ , δ ₂ , δ ₃ }={0.33, 0.33, 0.33, 65.6, 14.8, 117.75, 20.4, 7.34, 61.39}, use the parameter value θ ^t to calculate the posterior probability β _j (x), j=1, 2, .. C, t is the number of iterations. Fixed mean, calculate each parameter (weight, covariance):

均值为μ_j，方差为δ_j的高斯分布的概率密度函数为： The probability density function of a Gaussian distribution with mean μ _j and variance δ _j is:

$g g ((x x,, {μ μ}_{j j},, {δ δ}_{j j})) = = \frac{11}{{((22 π π))}^{11 / / 22} {| | {δ δ}_{j j} | |}^{11 / / 22}} exp exp ((- - [[\frac{11}{22} {((x x - - {μ μ}_{j j}))}^{T T} {δ δ}_{j j}^{- - 11} ((x x - - {μ μ}_{j j}))]]))$

代入公式： Into the formula:

${β β}_{j j} ((x x)) = = \frac{{a a}_{j j} g g ((x x,, {μ μ}_{j j},, {δ δ}_{j j}))}{{Σ Σ}_{k k = = 11}^{C C} g g ((x x,, {μ μ}_{k k},, {δ δ}_{k k}))},, j j = = 1,2 1,2,, 33,,$

${a a}_{j j}^{new new} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i})),, j j = = 1,2 1,2,, 33,,$

${δ δ}_{j j}^{new new} = = \frac{{Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i})) (({x x}_{i i} - - {μ μ}_{j j})) {(({x x}_{i i} - - {μ μ}_{j j}))}^{T T}}{{Σ Σ}_{i i = = 11}^{N N} {β β}_{j j} (({x x}_{i i}))},, j j = = 1,2 1,2,, 33$

重复计算β_j〔x〕和各参数，若||θ^new-θ||＜ζ，ζ为误差值(根据精度需求选取，本实施例取于10^-5)，迭代停止，则得到当前均值为{m₁，m₂，m₃}＝{65.6，14.8，117.75}时的其他参数估计结果{a₁，a₂，a₃，δ₁，δ₂，δ₃}＝{0.2992，0.3022，0.3985，21.45，7.82，63.00} Repeat the calculation of β _j [x] and each parameter, if ||θ ^new -θ||<ζ, ζ is the error value (selected according to the accuracy requirement, this embodiment takes 10 ^-5 ), the iteration stops, and the current mean value is obtained Estimation results of other parameters { _{a 1} _, a ₂ , a ₃ , δ ₁ , δ ₂ , δ ₃ _} ₌ {0.2992, 0.3022, 0.3985, 21.45, 7.82, 63.00}

(4)图像分类精度指标kappa系数计算：由得到的相应均值的权重和方差计算出相应的kappa系数， (4) Calculation of the kappa coefficient of the image classification accuracy index: calculate the corresponding kappa coefficient from the weight and variance of the corresponding mean value obtained,

$kappa kappa = = 11 - - \frac{{Σ Σ}_{k k = = 11}^{C C} {P P}_{kadd kadd}}{{Σ Σ}_{k k = = 11}^{C C} [[{P P}_{k k} (({P P}_{kk kk} + + {P P}_{kadd kadd}))]]}$

计算得到当前图像的kappa系数为0.56。 The calculated kappa coefficient of the current image is 0.56.

Claims

1. a Remote Sensing Image Quality evaluation method for feature based distribution, is characterized in that comprising the steps:

(1) obtain remote sensing images, intuitively obtained the atural object classification number C of image by image naked eyes;

(2) atural object all kinds of in image is sampled respectively, according to sample, generate gray level co-occurrence matrixes; According to co-occurrence matrix, calculate the sample characteristics of all kinds of atural object, comprise energy, contrast, unfavourable balance square, entropy and correlativity;

(3) calculate average and the variance of all kinds of ground object sample feature, as the initial parameter value of parameter estimation, carry out gauss hybrid models parameter estimation, obtain the class weight α of all kinds of atural object of remote sensing images _k, covariance sigma _k, k=1,2 ... C;

(4) computed image nicety of grading index coefficient k appa, comprises following sub-step:

(4.1) categorised decision face h is set up _j:

h_{j} = \frac{μ_{j} + μ_{j + 1}}{2}, 1 \leq j < C,

Wherein, μ _j(j=1,2 ..C) is average, is the center of the corresponding Gaussian spatial of jth class atural object;

(4.2) P is calculated _kk, P _kadd:

If k=1, then

\{\begin{matrix} P_{kk} = P_{k} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) \end{matrix};

If 1 < k < C,

\{\begin{matrix} P_{kk} = P_{k} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 2}^{C} P_{j} [Φ (\frac{h_{k} - μ_{k}}{σ_{k}}) - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix};

If k=C,

\{\begin{matrix} P_{kk} = P_{k} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \\ P_{kadd} = Σ_{j = 1}^{C - 1} P_{j} [1 - Φ (\frac{h_{k - 1} - μ_{k}}{σ_{k}})] \end{matrix};

In formula, Φ is gauss of distribution function; P _kit is the class weight calculated by step (3); P _kkbelong to kth class under representing truth, be also correctly classified as the pixel ratio of kth class; P _kaddrepresent reality and do not belong to kth class, and be divided into the pixel ratio of kth class by mistake;

(4.3) nicety of grading index k appa is calculated:

kappa = 1 - \frac{Σ_{k = 1}^{C} P_{kadd}}{1 - Σ_{k = 1}^{C} [P_{k} (P_{kk} + P_{kadd})]} .

2. method according to claim 1, is characterized in that, carries out gauss hybrid models parameter estimation and comprise following sub-step in described step (3):

(2.1) every class average in present image is established to be respectively m ₁... m _c, then initial parameter value is made to be θ °={ a ₁..., a _c, m ₁... m _c, δ ₁... δ _c; Wherein a ₁..., a _cfor the class weight of every class atural object, this weight is determined according to image atural object prior distribution by user, as without priori, is just defaulted as and is evenly distributed, a _i=1/C, meets δ ₁... δ _cfor eigenwert variance, C is classification number; G (x, μ _j, δ _j) be the probability density function of Gaussian distribution:

g (x, μ_{j}, δ_{j}) = \frac{1}{{(2 π)}^{1 / 2} {| δ_{j} |}^{1 / 2}} \exp (- [\frac{1}{2} {(x - μ_{j})}^{T} {δ_{j}}^{- 1} (x - μ_{j})]),

In formula, μ _jfor average, δ _jfor variance;

(2.2) μ is made _j=m _j, j=1,2..C, iterative computation

β_{j} (x) = \frac{a_{j} g (x, μ_{j}, δ_{j})}{Σ_{k = 1}^{C} g (x, μ_{k}, δ_{k})}, j = 1,2, . . . C,

(formula 1)

a_{j} = \frac{1}{n} Σ_{i = 1}^{N} β_{j} (x_{i}), j = 1,2 . . C,

(formula 2)

δ_{j} = \frac{Σ_{i = 1}^{N} β_{j} (x_{i}) (x_{i} - μ_{j}) {(x_{i} - μ_{j})}^{T}}{Σ_{i = 1}^{N} β_{j} (x_{i})}, j = 1,2, . . . C,

(formula 3)

N is the number of pixels of all kinds of atural object in sample; X is the variable representing pixel, β _jx () represents that current pixel belongs to the posterior probability of jth class;

By (formula 1)-(formula 3) iterative computation, until reach predetermined iteration precision ζ, iteration stopping, being fixed average is m _jtime class weight α _jwith variance δ _jestimated result.

3. method according to claim 2, is characterized in that, in described step (2.2), iteration stopping condition is: || θ ^t+1-θ ^t|| < ζ,

Wherein, θ ^tbe the t time iteration result, t is iterations, and ‖ ‖ represents and asks Euclidean distance, and predetermined iteration precision ζ is 10 ^-5-10 ^-6, according to computational accuracy and speed of convergence choosing comprehensively.