CN101478693A - 一种对星载光学遥感图像压缩质量进行评价的方法 - Google Patents

Abstract

本发明提供了一种对星载光学遥感图像压缩质量进行评价的方法。方法包括以下步骤：获取模拟原始图像；获取经压缩后再解压的图像；分别分析模拟原始图像和压缩后图像的动态范围；分别分析模拟原始图像和压缩后图像的纹理细节；分析压缩对图像的应用影响；分析模拟原始图像和压缩后图像的相互关系；综合以上对模拟原始图像和压缩后图像进行的各项分析结果之后，对图像压缩质量进行评价。本发明针对卫星传感器的特性生成模拟图像，采用客观评价为主、主观评价为辅的综合方法全面、系统、有针对性地评价图像压缩质量，构建了星载光学遥感图像压缩质量评价的完整体系。

Description

一种对星载光学遥感图像压缩质量进行评价的方法

技术领域

本发明属于卫星遥感图像处理与分析领域，尤其涉及一种对星载光学遥感图像压缩质量进行评价的方法。

背景技术

为了有效降低数据传输带宽或码速率，缓解地面接收站的压力，数据压缩技术被广泛应用于遥感卫星领域，特别是几乎所有的高分辨率遥感卫星均采用了有损压缩技术。虽然遥感图像压缩能带来传输和存储效率的提高，但不能以过多损失图像的质量为代价，必须使压缩造成的图像质量损失控制在遥感图像应用可接受的范围内。因此，必须提出合理的遥感图像压缩质量评价方法。

目前，压缩领域的研究主要集中在压缩算法的编码性能上，对图像压缩质量的评价方法研究较少，特别地，针对遥感卫星数据压缩质量评价，研究更少。本发明针对星载光学遥感图像压缩质量评价进行。

现有技术方案1如下：

1)图像获取

国外中高分辨率卫星(如SPOT5、IKONOS等)的图像和航空影像(75cm)。

2)评价方法

(1)主观评价方法

主观评价方法就是让观察者根据一些事先规定的评价尺度或自己的经验，对测试图像按视觉效果提出质量判断，并给出质量分数，对所有观察者给出的分数进行平均，所得的结果即为图像的主观质量评价。

(2)客观评价方法

峰值信噪比PSNR

PSNR定义为与均方误差MSE相关的一个量，PSNR和MSE分别定义为

PSNR＝-10log₁₀MSE/(255×255)   (9)

$\mathrm{MSE}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[f(i,j)-g(i,j)]}^2$

其中M，N分别为图像每行和每列的元素数，f(i，j)为原始图像在(i，j)点的灰度值，g(i，j)为解压图像在(i，j)点的灰度值。

相似度

$\mathrm{XSD}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}f(i,j)g(i,j)}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}f{(i,j)}^2}\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}g{(i,j)}^2}}---\left(10\right)$

M，N、f(i，j)和g(i，j)的意义同上。

拉氏逼真度

$L\_\mathrm{Fidelity}=-\underset{i=2}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=2}{\overset{N-1}{\mathrm{\Σ}}}[f\′(i,j)g\′(i,j)]/\underset{i=2}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=2}{\overset{N-1}{\mathrm{\Σ}}}{f\′}^2(i,j)---\left(11\right)$

其中

f′(i，j)＝f(i+1，j)+f(i-1，j)+f(i，j-1)+f(i，j+1)-4f(i，j)

g′(i，j)＝g(i+1，j)+g(i-1，j)+g(i，j-1)+g(i，j+1)-4g(i，j)。

拉氏均方误差

其中f′(i，j)和g′(i，j)意义同上。

现有技术的技术方案二如下：

1)灰度均值

图像灰度均值可以从整体上反映图像的亮度，其计算公式为

$\mathrm{\μ}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}f(i,j)---\left(13\right)$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值；

2)灰度方差

灰度方差反映了图像灰度的层次信息，其计算公式为

${\mathrm{\σ}}^2=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[\mathrm{\μ}-f(i,j)]}^2---\left(14\right)$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值，μ为图像的灰度均值；

3)信噪比

信噪比是图像中的有用信息与噪声信号的比值。LSD_max为局域方差极大值。

$\mathrm{SNR}=\mathrm{\μ}/\sqrt{{\mathrm{LSD}}_{\max }}---\left(16\right)$

4)对比度

图像对比度是基于灰度共生矩阵的评价图像纹理的参数，它反映了图像中的目标与背景相比可辨认的清晰程度，计算公式为

$f=\underset{n=0}{\overset{L-1}{\mathrm{\Σ}}}{n}^2\left\{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{^j=0}{\overset{L-1}{\mathrm{\Σ}}}\stackrel{\mathrm{\Λ}}{p}(i,j)\right\}---\left(16\right)$

其中，|i-j|＝n。p(i，j)的定义见式(28)

5)角二阶矩

图像角二阶矩是灰度共生矩阵的二次统计量，具体为像素值平方的和，也称为能量，是图像灰度分布均匀性的度量。其计算公式为

$f=\underset{^i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}\stackrel{\mathrm{\Λ}}{p^2}(i,j)---\left(17\right)$

p(i，j)的定义见式(28)

6)信息熵

熵值是图像所具有信息量的度量，纹理的复杂度越高就意味着图像信息量越大，其熵值也越大，这正是图像质量所最需要表征的。其计算公式为

$H=-\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}p\left(x_i\right){\log }_{2}p\left(x_i\right)---\left(18\right)$

其中，p(x_i)为灰度等于x_i的概率密度。

7)边缘信号

边缘是图像关于形状特征和细节的重要信息。用45°，135°两个斜线方向上的归一化边缘算子E₁、E₂分别对图像进行卷积计算，相加后得到图像边缘e(x，y)，即

e(x，y)＝E₁(f(i，j))+E₂(f(i，j))＝E(f(i，j))   (19)

定义 ${\mathrm{\σ}}_e^2=\frac{1}{M\×N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{e}^2(i,j)$ 为边缘e(i，j)的平均能量。

现有技术的缺点

1)现有技术一主要存在以下缺点：

模拟图像的获取存在一定的问题。国外中高分辨率卫星几乎均采用了星上实时压缩技术，其图像均为解压图像，直方图可能不连续。而未经压缩的遥感图像，其直方图应是连续的。因此，用解压后的卫星图像，作为压缩图像质量评价的模拟原始图像，是不合理的。

模拟图像的特性和数量未做深入的研究。

主观评价方法对参评人员的身份和数目未做深入的研究。

图像的显示(放大倍率、增强方法)方式没有相应研究。

客观评价方法过于单调，不能构成对图像压缩质量的全面评价。

2)现有技术二主要存在以下缺点：

现有技术二中提到的评价方法是通用的图像质量评价方法，针对性不强。

这些评价方法只是针对于图像纹理细节方面的评价，而压缩对图像质量的影响，不仅限于对图像纹理细节的影响。

发明内容

本发明所要解决的技术问题是对星载光学遥感图像压缩质量进行全面系统的评价。

本发明创造性地提出了如何科学生成压缩质量评价的模拟图像的方法。将航片进行数字化扫描，所得图像直方图连续；并且结合卫星传感器特性，生成的模拟图像的动态范围、噪声特性、量化比特数等与真实的卫星尽量一致，因而可以作为压缩质量评价的模拟图像。

本发明对获取的模拟图像的特性和数量做了详尽描述。模拟图像具有代表性，充分反映遥感图像的地域特征、地物类型和成像时相等信息。

本发明之主观评价对参评人员的身份做了明确规定，即参评人员均为训练有素的遥感图像判读专家，并规定了参评专家的最少数目；明确指出了主观评价的尺度。

本发明之客观评价创造性地从直方图的动态范围、图像的纹理细节、应用影响和图像的相互关系四个方面对图像压缩质量进行全面、系统、有针对性地评价。

本发明的技术方案如下：

一种对星载光学遥感图像压缩质量进行评价的方法，其特征在于包括如下步骤：

获取模拟原始图像；

获取压缩后的图像；

分别分析模拟原始图像和经压缩后再解压的图像的动态范围，其中包括计算图像的低端、中端、高端灰度值以及图像的标准方差和灰度均值；

分别分析模拟原始图像和经压缩后再解压的图像的纹理细节，其中包括计算图像的分块标准方差、信噪特性、信息熵值、角二阶矩、对比度、边缘信号；

分析压缩对图像的应用影响，其中包括计算经压缩后再解压的图像的几何畸变；

分析模拟原始图像和经压缩后再解压的图像的相互关系，其中包括计算图像的差值图像、直方图相关程度、峰值信噪比PSNR、图像相关系数、综合评价因子；

综合以上对模拟原始图像和经压缩后再解压的图像进行的各项分析结果之后，对图像压缩质量进行客观评价。

所述模拟原始图像的获取可以通过航片扫描间接获取，或者用与星上传感器技术指标相同或相近的机载传感器的图像作为模拟原始图像；

并且所述模拟原始图像的数目不少于10幅。

所述图像的低端、中端、高端灰度值分别定义为累计直方图面积为为5％，50％和95％所对应的灰度值；

所述图像的灰度均值的计算公式为

$\mathrm{\μ}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}f(i,j)---\left(20\right)$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值；

所述图像的标准方差的计算公式为

$\mathrm{\σ}=\sqrt{\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[\mathrm{\μ}-f(i,j)]}^2}---\left(21\right).$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值，μ为图像的灰度均值；

所述图像的分块标准方差的计算方法如下：以n*n的窗口作为分块单位求取标准方差，其中n为整数，最后计算各分块图像标准方差的平均值作为图像标准方差；

所述峰值信噪比PSNR的计算公式如下

PSNR＝-10log₁₀MSE/(255×255)   (22)

其中MSE为均方误差，由下式计算

$\mathrm{MSE}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[f(i,j)-g(i,j)]}^2$

其中M，N分别为图像每行和每列的元素数，f(i，j)为原始图像在(i，j)点的灰度值，g(i，j)为解压图像在(i，j)点的灰度值；

所述图像相关系数的计算方法如下：

设图像匹配的目标窗口灰度矩阵为G＝(g_i，j)i＝1，2，...，m；j＝1，2...，n；其中m与n是G的行列数，取值为奇数。与G相应的灰度函数为g(x，y)，(x，y)∈D。搜索区灰度矩阵为G′＝(g′_i，j)i＝1，2，...，k；j＝1，2...，l；其中k与l是G′的行列数，取值为奇数；与G′相应的灰度函数为g′(x，y)，(x，y)∈D′；

采用图像相关系数作为图像匹配的匹配测度，相关系数是标准化的协方差函数，定义为

$\mathrm{\ρ}=S_{\mathrm{xy}}/\sqrt{S_{\mathrm{xx}}{S}_{\mathrm{yy}}}---\left(23\right)$

其中，S_xy是两幅影像图像窗口之间的协方差函数，S_xx、S_yy则分别为两幅影像上图像窗口内的方差；

由离散灰度数据对相关系数的估计为

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}-\stackrel{\‾}{g})(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i,j}-\stackrel{\‾}{g})}^2\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}^2}}---\left(24\right)$

其中 $\stackrel{\‾}{g}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j},$ ${\stackrel{\‾}{g}}_{r,c}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}$

对上式进行优化后得到

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}{g}_{i+r,j+c}^{\′}-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i,j}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2][\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i+r,j+c}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}\right)}^2]}}---\left(25\right)$

在以上式中，(c，r)为搜索区的待匹配窗口中心坐标。若ρ(c₀，r₀)>ρ(c，r)(c≠c₀，r≠r₀)，则c₀，r₀为搜索区图像相对于目标区图像位移的行、列参数。对于一维相关，有r＝0；

所述综合评价因子的计算方法如下

综合评价因子为峰值信噪比PSNR与相关系数ρ的乘积PSNR*ρ。

在所述客观评价结果的基础上再结合主观评价的结果以对图像压缩质量给出全面的评价。

所述主观评价尺度包括辐射精度、清晰度、纹理、几何特性、色调；

在评价过程中进一步采用相对尺度判据，用其评价一组图像对应原始图像的相对质量；主观评价结果表示为平均意见的分数，由以下公式表示：

$\stackrel{\‾}{c}=\frac{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{c}_{i,j}}{I*J}---\left(26\right)$

其中，I为判读人员的人数，应多于10人，J为图像的数量，c_i，j代表第i名判读人员对第j幅图像的打分，当采用10分制时，c_i，j取1，2，...10；

判读人员在评判图像质量时，可放大图像观看细节，还可采用图像增强方式考察信息的微小变化。

本发明技术方案带来了如下有益效果：2005年7月，航天504所和中国资源卫星应用中心在西安组织召开了“资源一号卫星02B星(CBERS-02B)HR相机和环境减灾1A/B星(HJ-1A/B卫星)CCD相机图像数据压缩质量评价”专项评审与研讨会，评价方案采用本发明技术方案。专家评审发现CCD相机8：1的压缩对图像的处理和应用影响很大。卫星研制部门及时更新方案，采用4：1的压缩算法。2006年2月通过新一轮评价，专家一致认为：新压缩方案不仅满足广大遥感用户的应用要求，同时也符合数据下传的要求。本发明促成了HJ-1A/B卫星CCD相机压缩方案的改进和完善，即从8：1的压缩比降低至4：1，为进一步保证图像质量起到积极作用。此发明技术能全面评价压缩图像质量的影响，有利于在分辨率和压缩比之间达成和谐的统一，这对合理设计星上实时压缩方案具有很强的指导意义。质量评价的质量因子可以广泛推广。截至目前，这一发明成果大量推广到卫星在轨测试的质量评价，MTFC图像补偿的质量评价等。

附图说明

图1是城市图像。

图2是农田图像。

图3分前后和分奇偶压缩的效果比较

具体实施方式

在本部分中将给出本发明的实施例，并予以详细描述，以便更好地说明本发明特点和功能，使得本发明更易于理解，而不是用来限定本发明的保护范围。

首先介绍一下灰度共生矩阵、基于相关系数的图像匹配、基于最小二乘法的图像匹配。

灰度共生矩阵

从灰度值为i的像素点出发，距离为(Dx，Dy)的另一个像素点出现灰度值为j，定义这2个灰度值在整个图像中发生的概率，或称为频度即为灰度共生矩阵，用数学表示则为

P(i，j，δ，θ)＝{(x，y)|f(x，y)＝i，f(x+Dx，y+Dy)＝j}   (27)

其中(x，y)为二维图像的行、列号，x＝0，1，…M-1，y＝0，1，…N-1；

i，j＝0，1，…L-1，当图像采用8比特量化时L＝256；

δ取值为Dx，Dy的欧氏距离，实际应用中一般取值为1，即相邻元素；

θ为方向，一般取值0°，45°，90°，135°等4个方向。

采用以上方法，就将(x，y)的空间坐标转化为(i，j)的“灰度对”的描述，这一概念可以理解为像素对或灰度级对的直方图，要求像素对的距离不变，像素的灰度差不变。

根据上述定义，所构成的灰度共生矩阵是一个集合，集合中的一个元素P(i，j，δ，θ)为第i行，第j列矩阵元素，表示所有在θ方向上，相邻间隔为δ的像素，一个为灰度i值，另一个为灰度j值得相邻点对数。

可以看出，灰度共生矩阵反映了图像灰度关于方向，相邻间隔，变化幅度的综合信息，可以作为分析图像基元和排列结构的信息。一幅卫星图像在确定了方向和距离后，也就确定了其共生矩阵。其重要特点是以主对角线为对称轴，两边对称。若该共生矩阵的主要数值集中在矩阵的对角线附近，就说明该方向上灰度变化的频数较低，不同地物反差不甚明显，纹理相对粗糙；若主对角线数值基本为零，数据主要集中在左下角和右上角，就说明此方向上灰度变化较为剧烈，纹理细致。

为了利用灰度共生矩阵所提供的图像灰度方向，间隔和变化幅度的信息，在共生矩阵的基础上抽取的纹理特征，称为二次统计量。为了使二次统计量的表达式清晰简单，可作正规化处理如下

$\stackrel{\mathrm{\Λ}}{p}(i,j)=P(i,j)/R---\left(28\right)$

R为正规化常数，含义为相邻点对的组合数，在θ＝0°，δ＝1，图像大小为N*N时，P(i，j)的元素之和为2N(N-1)，正规化处理后，其元素之和为1。

基于相关系数的图像匹配

若图像匹配的目标窗口灰度矩阵为G＝(gi，j)i＝1，2，...，m；j＝1，2...，n；m与n是G的行列数，一般情况下取奇数。与G相应的灰度函数为g(x，y)，(x，y)∈D。搜索区灰度矩阵为G′＝(g′i，j)i＝1，2，...，k；j＝1，2...，1；k与1是G′的行列数，一般情况下也取奇数。与G′相应的灰度函数为g′(x，y)，(x，y)∈D′。

图像匹配的匹配测度很多，相关系数是最常使用的一种匹配测度。相关系数是标准化的协方差函数，定义为

$\mathrm{\ρ}=S_{\mathrm{xy}}/\sqrt{S_{\mathrm{xx}}{S}_{\mathrm{yy}}}---\left(29\right)$

其中，S_xy是两幅影像图像窗口之间的协方差函数，S_xx、S_yy则分别为两幅影像上图像窗口内的方差。从式(1)不难知道，相关系数|ρ|≤1。

由离散灰度数据对相关系数的估计为

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}-\stackrel{\‾}{g})(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i,j}-\stackrel{\‾}{g})}^2\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}^2}}---\left(30\right)$

$\stackrel{\‾}{g}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}---\left(31\right)$

${\stackrel{\‾}{g}}_{r,c}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}---\left(32\right)$

考虑到计算工作量，相关系数的优化公式为

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}{g}_{i+r,j+c}^{\′}-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i,j}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2][\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i+r,j+c}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}\right)}^2]}}---\left(5\right)$

(c，r)为搜索区的待匹配窗口中心坐标。若ρ(c0，t0)>ρ(c，r)(c≠c0，r≠r0)，则c0，r0为搜索区图像相对于目标区图像位移的行、列参数。对于一维相关，有r＝0。

基于最小二乘法的图像匹配

两个二维图像之间的变形，不仅仅存在着相对移位，而且还存在着图形变形，为了简化处理，假设压缩只使图像在水平和垂直方向发生偏移而不存在旋转。基于最小二乘法的图像匹配需要迭代求解，直到结果满意为止。假设原始图像为g₁(x，y)，对应的随机噪声为n₁(x，y)；解压缩图像为g₂(x，y)，对应的随机噪声为n₂(x，y)；g₂(x，y)相对于g₁(x，y)发生的零次几何变形为a，b，则：

g₁(x，y)+n_i(x，y)＝g₂(x+a，y+b)+n₂(x，y)   (33)

为解求相对位移da，db，对公式(1)进行线性化，误差方程为：

v＝C₁da+C₂db-Δg   (34)

其中： $C_1=\stackrel{.}{g_x},$ $C_2=\stackrel{.}{g_y};$ 由于在数字图像匹配中，灰度均按照规则网格的离散阵列，且采样间隔为常数Δ，故偏导可用差分代替。

${\stackrel{.}{g}}_x={\stackrel{.}{g}}_I(I,J)=\frac{1}{2}[g_2(I+1,J)-g_2(I-1,J)]---\left(35\right)$

${\stackrel{.}{g}}_y={\stackrel{.}{g}}_J(I,J)=\frac{1}{2}[g_2(I,J+1)-g_2(I,J-1)]---\left(36\right)$

初值a＝0，b＝0，观测值Δg是相应像元的灰度差。

逐个像元(在目标区内)建立误差方程

v＝CX-L   (37)

其中X＝[da   db]，C＝[C₁ C₂]，L＝Δg

从而X＝(C^TC)^-1(C^TL)   (38)

在以上内容的基础上，在下面给出本发明的具体实施例

一、模拟图像的获取

1)图像特性

遥感图像的获取应有代表性，应充分反映遥感图像的地域特征和地物类型，成像时相等信息。地物类型应该涵盖的地物包含：城市，山地，沿海地带，沙漠，湖泊，植被等。成像时相应该涵盖赤道至北纬60度的范围。

由于需要评价压缩前后图像质量的变化，因此所选图像不能经过压缩。若难以获取未经压缩的原始图像，可以通过航片扫描间接获取。

例如，对航摄比例尺为1：25000的航片进行数字化扫描，扫描分辨率为100微米，就可得到空间分辨率为2.5m的图像。

2)图像数目

模拟图像数目不小于10幅，有条件时图像数量超过20幅。

二、遥感图像压缩质量评价方法

(1)主观评价方法

判断图像质量的一种方法是主观评价，在电视领域常用的评价尺度有两种，即品质尺度和妨碍尺度。前者由一些未受过训练的、对图像质量评价不内行的观测者判断，此时得到的图像质量代表平均观测者的一般感觉；而后者则由训练有素的专家来进行。这些人在图像处理方面是有经验的，能够在图像质量方面提出严格的判断。专家观测者往往具有注意细小程度图像质量下降的能力，而这些正是外行的观测者所缺乏的。由于遥感图像不仅用于目视解译，也用于计算机解译和各种目的用途，因此采用电视领域妨碍尺度的评价尺度，聘请遥感图像判读专家进行遥感解压图像的评价较为恰当。

为了保证统计上的可靠性，在主观质量测试中通常至少应有10名评判者。专家在评价图像质量时，需要放大图像观看细节(通常放大为原图的200％倍)；同时也需要采用图像增强方式考察信息的微小变化。

本发明的主观评价尺度包括下面几个方面：

辐射精度

辐射精度是反映图像辐射状态的指标，同一地区的不同图像，灰度分布范围越大，则表明图像信息越丰富；由于不同地物的电磁辐射(反射或发射)强度不同，在影像上构成了灰阶的等级，灰阶代表了由最暗到最亮之间不同亮度的层次级别。一般在目视判读时，影像上的中间层级越多，所能够呈现的画面效果也就越细腻。

清晰度

图像的清晰度表示图像内部边缘的清晰程度，若图像质量较好，将反映为图像的边缘清晰，边缘相邻像素的灰度差异较大，易于人眼的判识；若图像质量不好，将反映为图像的边缘模糊，边缘相邻像素的灰度差异较小，不易于人眼的判识。

纹理

也叫影像结构，是指与色调配合看上去平滑或粗糙的纹理的粗细程度，即图像上目标物表面的质感。草场及牧场看上去平滑，成材的老树林看上去很粗糙。海滩的纹理能反映沙粒结构的粗细，沙漠中的纹理可表现沙丘的形状以及主要风系的风向。评价一幅卫星图像质量优劣与否，首先应该从其包含的信息量出发。不同卫星获取的同一幅遥感图像上的同一种地物，若纹理粗糙，说明该图像包含信息少。若图像上充满了细纹理，说明包含了更多的信息。

几何特性

几何特性指目标物在影像上所呈现的几何图形，在遥感影像上能看到的是目标物的顶部或平面形状。地物在影像上的形状受空间分辨率、比例尺、投影性质等影响。如果地物产生变形，说明产生了几何畸变。

色调

色调指影像上黑白深浅的程度，是地物电磁辐射能量大小或地物波谱特征的综合反映。色调用灰阶(灰度)表示，同一地物在不同波段的图像上会有很大差别；同一波段的影像由于成像时间和季节的差异，地物光谱反射或辐射量不同。即使同一地区同一地物的色调也会不同。同一波段、同一时相的影像，可以观察两种传感器的光谱响应是否一致。

在评价过程中进一步采用相对尺度判据，用其评价一组图像对应原始图像的相对质量。主观评价结果表示为平均意见的分数，由以下公式表示：

$\stackrel{\‾}{c}=\frac{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{c}_{i,j}}{I*J}---\left(40\right)$

其中，c_i，j代表第i名判读人员对第j幅图像的打分，当采用10分制时，c_i，j取1，2，...10。

(2)客观评价方法

1)图像的动态范围

低端、中端、高端灰度值

图像的低端、中端、高端灰度值分别定义为累计直方图面积为5％，50％和95％所对应的灰度值。

标准方差

这一参数反映了图像灰度层次的丰富程度，从直方图角度来说，它反映了直方图的大致分布宽度。在图像比较分析中，图像的方差较大，说明图像灰度层次丰富，在目视效果中，地物更加易于识别和分类，图像质量较为理想。其计算公式为

$\mathrm{\σ}=\sqrt{\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[\mathrm{\μ}-f(i,j)]}^2}---\left(41\right)$

灰度均值

本发明技术方案中用到的灰度均值算法与现有技术方案二中提到的灰度均值相同。

2)图像的纹理细节

分块标准方差

此处定义的图像分块标准方差与图像动态范围一节中标准方差的定义完全一致，但在具体操作时有所不同。即采用分块思想，以n*n(n为3、5、7、9或其它整数)的窗口为单位求取标准方差，最后计算各分块图像标准方差的平均值作为图像标准方差。分块标准方差可以体现图像局部纹理灰度层次的丰富程度。

信噪特性

可以使用比特分割分析图像的信噪特性。所谓比特分割就是将图像数据分成不同的比特位，依此取出某一比特位上的值(0或1)形成二值图像。大量研究及实验表明，通常每像素为8bit的灰度图像扫描数据中，低位的两个比特基本属于噪声信号。

信息熵值

本发明技术方案中用到的信息熵值算法与现有技术方案二中提到的信息熵值算法相同。

角二阶矩

本发明技术方案中用到的角二阶矩算法与现有技术方案二中提到的角二阶矩算法相同。

对比度

本发明技术方案中用到的对比度算法与现有技术方案二中提到的对比度算法相同。

边缘信号

本发明技术方案中用到的边缘信号算法与现有技术方案二中提到的边缘信号算法相同。

3)应用影响

几何畸变

解压后的图像与原始图像可能发生几何位置的偏移，可以使用不同算法得到该偏移值。参考前面对基于相关系数的图像匹配的说明，使用相关系数法可以获取整像元级意义上的几何位置偏移；参考前面对基于最小二乘法的图像匹配的说明，使用最小二乘法可以获取子像元级意义上的几何位置偏移。为了保证匹配精度，可以在图像上选取若干控制点获得位置偏移的程度；也可以在图像上等间隔布设网格点，当网格点符合匹配要求时，可以参与位置偏移的计算。

4)图像的相互关系

关联的图像指标是通过解压后图像与原始图像进行比较获取的。从获取的途径看可以分为两类：按照像素逐一比较(大多为减法)法和基于相关系数的匹配法。

差值图像

差值图像代表解压后的图像与原始图像的差别，见下式所示：

$\mathrm{Dif}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}|g(i,j)-f(i,j)|---\left(42\right)$

差值图像的平均值

$\mathrm{Diff}\_\mathrm{ave}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\mathrm{\Δf}(i,j)---\left(43\right)$

差值图像的最大值

Diff_max＝max(Δf(i，j))   (44)

直方图相关程度

直方图相关程度，反映解压后的图像与原始图像的直方图接近的程度。图像直方图是用于表达图像灰度值分布情况的统计图表，其横坐标是灰度值，纵坐标是出现这一灰度值的数目，即为取不同灰度值的面积或像素在整幅图像上所占的比例。假设原始、解压缩后图像的灰度概率密度函数分别为p(k)和q(k)，其中k(k＝0，1，2…K-1，当采用8比特量化时，K＝256)为灰度级别。则直方图的相关系数ρ_H为：

${\mathrm{\ρ}}_H=\frac{\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}(p\left(k\right)-\stackrel{\‾}{p})(q\left(k\right)-\stackrel{\‾}{q})}{\sqrt{\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{(p\left(k\right)-\stackrel{\‾}{p})}^2\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{(q\left(k\right)-\stackrel{\‾}{q})}^2}}---\left(45\right)$

峰值信噪比PSNR

本发明技术方案中用到的峰值信噪比PSNR算法与现有技术方案一中提到的峰值信噪比PSNR相同。

图像相关

图像匹配在测绘遥感领域得到广泛的应用，为图像融合、超分辨率技术、自动几何校正及自动生成数字地面模型(DEM)等应用提供基础。不同的图像匹配算法采用不同的测度，目前较多使用相关函数测度、协方差函数测度、相关系数测度、差平方和测度及差绝对值和测度。相关系数是灰度线性变换的不变量(证明从略)，可以有效消除由于成像时间造成的图像之间亮度与对比度的不同影响，能够较好用于自动几何校正的匹配过程。图像的相关系数ρ是标准化的协方差函数，详细定义见附录2所示。

综合评价因子

综合评价因子PSNR*ρ为峰值信噪比与相关系数的乘积。不同图像的冗余度不同，即使采用相同的压缩方法和压缩比，压缩效果可能会差异较大。在压缩方法和压缩比相同的情况下，对于纹理较丰富的图像，峰值信噪比较低、相关系数较高；而对于纹理较单一的图像，峰值信噪比较高、相关系数较低。

大量实验证明，综合评价因子能较为稳定地反映压缩前后图像质量的变化。

举例说明：

图1(a)和图2(a)分别为城市和农田的遥感图像。与图2(a)相比，图1(a)纹理较丰富。将两者用SPIHT算法进行8：1压缩，解压图像分别为图1(b)和图2(b)。

经计算得出，城市图像的峰值信噪比PSNR为35.61dB，压缩前后图像的相关系数ρ为0.91，综合评价因子PSNR*ρ为32.55dB；农田图像的峰值信噪比PSNR为37.81dB，压缩前后图像的相关系数ρ为0.86，综合评价因子PSNR*ρ为32.38dB。城市图像的PSNR大于农田图像的PSNR，而ρ小于农田图像的ρ，但综合评价因子PSNR*ρ十分接近。说明PSNR*ρ是比PSNR和ρ更稳定的图像压缩质量评价因子。

5)卫星的具体数传方式对图像质量的影响

对于卫星下传图像，在卫星内部传输可能采用若干方案，其中：分奇偶压缩和分前后压缩两种方案被普遍采用。分奇偶压缩是指将每片CCD数据分奇偶两路并行输出、压缩再合并下传；分前后压缩是指将每片CCD数据分成前后两块数据，分别压缩再下传，地面系统完成解压。

在环境减灾1A/B星(HJ-1A/B卫星)CCD相机图像数据压缩质量评价实验工作中，得出了如下结论：分奇偶压缩易造成虚假纹理(块状纹理)；而分前后压缩对纹理的损失较大，并随着压缩比的提高而逐步增大。

举例说明：

图3(a)为原始图像，分别采用分前后和分奇偶两种方案对原始图像进行8:1压缩，解压图像分别如图3(b)和图3(c)所示。对比图3中的(a)、(b)和(c)三张图可以看出，分前后压缩对图像纹理的损失较大，而分奇偶压缩引入了虚假块状纹理。

本发明的保护范围并不局限于上述内容，熟悉本技术领域的技术人员在本发明揭露的技术范围内的基础上所做的方案的变形、变化或者替换，都应涵盖在本发明保护范围之内。

Claims (7)

1、一种对星载光学遥感图像压缩质量进行评价的方法，其特征在于包括如下步骤：

获取模拟原始图像；

获取经压缩后再解压的图像；

分别分析模拟原始图像和经压缩后再解压的图像的动态范围，其中包括计算图像的低端、中端、高端灰度值以及图像的标准方差和灰度均值；

分别分析模拟原始图像和经压缩后再解压的图像的纹理细节，其中包括计算图像的分块标准方差、信噪特性、信息熵值、角二阶矩、对比度、边缘信号；

分析压缩对图像的应用影响，其中包括计算图像的几何畸变；

分析模拟原始图像和经压缩后再解压的图像的相互关系，其中包括计算图像的差值图像、直方图相关程度、峰值信噪比、图像相关系数、综合评价因子；

综合以上对模拟原始图像和经压缩后再解压的图像进行的各项分析结果之后，对图像压缩质量进行客观评价。

2、如权利要求1所述的方法，其特征在于：

所述模拟原始图像的获取可以通过航片扫描间接获取，或者用与星上传感器技术指标相同或相近的机载传感器的图像作为模拟原始图像；

并且所述模拟原始图像的数目不少于10幅。

3、如权利要求1所述的方法，其特征在于：

所述图像的低端、中端、高端灰度值分别定义为累计直方图面积为5％，50％和95％所对应的灰度值；

所述图像的灰度均值的计算公式为

$\mathrm{\μ}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}f(i,j)---\left(1\right);$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值；

所述图像的标准方差的计算公式为

$\mathrm{\σ}=\sqrt{\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[\mathrm{\μ}-f(i,j)]}^2}---\left(2\right)$

其中，M，N分别为图像每行和每列的元素数，0≤i≤N-1，1≤j≤M-1，f(i，j)为原始图像在(i，j)点的灰度值，μ为图像的灰度均值。

4、如权利要求1所述的方法，其特征在于：

所述图像的分块标准方差的计算方法如下：以n*n的窗口作为分块单位求取标准方差，其中n为整数，最后计算各分块图像标准方差的平均值作为图像标准方差。

5、如权利要求1所述的方法，其特征在于：

所述峰值信噪比PSNR的计算公式如下

PSNR＝-10log₁₀MSE/(255×255)                (3)

其中MSE为均方误差，由下式计算

$\mathrm{MSE}=\frac{1}{M*N}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{[f(i,j)-g(i,j)]}^2---\left(4\right)$

其中M，N分别为图像每行和每列的元素数，f(i，j)为原始图像在(i，j)点的灰度值，g(i，j)为解压图像在(i，j)点的灰度值；

所述图像相关系数的计算方法如下：

设图像匹配的目标窗口灰度矩阵为G＝(g_i，j)i＝1，2，...，m；j＝1，2...，n；其中m与n是G的行列数，取值为奇数；与G相应的灰度函数为g(x，y)，(x，y)∈D；搜索区灰度矩阵为G′＝(g′_i，j)i＝1，2，...，k；j＝1，2...，l；其中k与l是G′的行列数，取值为奇数；与G′相应的灰度函数为g′(x，y)，(x，y)∈D′；

采用图像相关系数作为图像匹配的匹配测度，相关系数是标准化的协方差函数，定义为

$\mathrm{\ρ}=S_{\mathrm{xy}}/\sqrt{S_{\mathrm{xx}}{S}_{\mathrm{yy}}}---\left(5\right)$

其中，S_xy是两幅影像图像窗口之间的协方差函数，S_xx、S_yy则分别为两幅影像上图像窗口内的方差；

由离散灰度数据对相关系数的估计为

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}-\stackrel{\‾}{g})(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i,j}-\stackrel{\‾}{g})}^2\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(g_{i+r,j+c}^{\′}-{\stackrel{\‾}{g}}_{r,c}^{\′})}^2}}---\left(6\right)$

其中 $\stackrel{\‾}{g}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j},$ ${\stackrel{\‾}{g}}_{r,c}=\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}$

对上式进行优化后得到

$\mathrm{\ρ}(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}{g}_{i+r,j+c}^{\′}-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i,j}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2][\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g^2}_{i+r,j+c}-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}\right)}^2]}}---\left(7\right)$

在以上式中，(c，r)为搜索区的待匹配窗口中心坐标；若ρ(c₀，r₀)>ρ(c，r)(c≠c₀，r≠r₀)，则c₀，r₀为搜索区图像相对于目标区图像位移的行、列参数；对于一维相关，有r＝0；

所述综合评价因子的计算方法如下

综合评价因子为峰值信噪比PSNR与相关系数ρ的乘积PSNR*ρ。

6、如权利要求1-5中任意一项所述的方法，其特征在于：

在所述客观评价结果的基础上再结合主观评价的结果以对图像压缩质量给出全面的评价。

7、如权利要求6所述的方法，其特征在于：

所述主观评价的尺度包括辐射精度、清晰度、纹理、几何特性、色调；

在评价过程中进一步采用相对尺度判据，用其评价一组图像对应原始图像的相对质量；主观评价结果表示为平均意见的分数，由以下公式表示：

$\stackrel{\‾}{c}=\frac{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{c}_{i,j}}{I*J}---\left(8\right)$

其中，I为判读人员的人数，应多于10人，J为图像的数量，c_i，j代表第i名判读人员对第j幅图像的打分，当采用10分制时，c_i，j取1，2，...10；

判读人员在评判图像质量时，可放大图像观看细节，还可采用图像增强方式考察信息的微小变化。

Images