系数估计

2024-06-02

系数估计（精选八篇）

系数估计篇1

合成孔径雷达(SAR)图像的散射系数和地物表面的电磁特性紧密相关,可以据此理解散射系数和图像数字值(DN,即图像像素值)之间的关系[1,2]。对于城市和山区,由于发生角反射,因此对应的图像往往具有较高的DN值;对于河流和湖波,由于发生镜面发射,雷达接收不到大部分电磁波,因此对应的图像具有较低的DN值;对于森林区域,由于发生漫散射,雷达可以接收一部分电磁波,因此对应图像的DN值在上述两者之间。总之,DN值依赖照射角度以及地物目标反射电磁波的特性,由此可对SAR图像进行解译。

对于特定地物表面以及成像条件,理想状态下的DN值应该能够正确反映地物表面的真实粗糙程度。但是,成像处理过程中(如传感器平台和信号处理子系统等等)总会发生一些误差甚至错误,因此,实际的DN值一般不能完全精确地反映地物表面的散射系数。如何通过实际的DN值获得散射系数的真实值,这一问题的解决直接依赖于雷达校正方程[3]。可以证明,校正方程具有线性形式,其中的校正因子由雷达系统参数(如天线增益和发射功率等等)决定。当然可以通过内部校准和外部校准来估计这些参数,由此获得校正因子,这就是校正因子估计的一般方法[3]。

基于散射系数的Rayleigh模型,本文提出了一种估计校正因子的简单而有效的方法。只要分别对线性形式的校正方程取均值和方差,并使用经典的Rayleigh模型描述散射系数的统计特性,校正因子就可以由图像的DN值估计得出。基于散射系数Rayleigh模型的校正因子估计具有解析的表达式,易于计算,便于实现。Monte Carlo仿真验证了这种估计方法的有效性。

1 校正因子估计式的推导

对于SAR图像,校正方程可以写为如下的线性形式[3]:

σ $_{d B}^{0}$ =A+B·lgDN (1)

其中,σ $_{d B}^{0}$ 表示以dB为单位的散射系数,lg代表常用对数,DN代表图像的像素值,A和B是待估计的校正因子。使用对数换底公式,经过适当化简,式(1)可以写为如下形式:

lnσ0=C+D·lnDN (2)

其中,ln代表自然对数,σ0表示不以dB为单位的散射系数(即 $σ_{d B}^{0} = 10 \cdot \lg σ^{0}) ‚ C = \frac{A l n 10}{10} ‚ D = \frac{B}{10}$ 。分别对式(2)两边取均值和方差,可得:

E(lnσ0)=C+D·E(lnDN) (3)

Var(lnσ0)=D2·Var(lnDN) (4)

其中,E(·)和Var(·)分别表示均值和方差。由式(3)和式(4)可以获得C和D,由此,校正因子A和B就可以通过下面两式估计得到:

$A = \frac{10}{\ln 10} [E (\ln σ^{0}) - E (l n D Ν) \cdot \sqrt{\frac{V a r (\ln σ^{0})}{V a r (l n D Ν)}}] (5)$

$B = 10 \sqrt{\frac{V a r (\ln σ^{0})}{V a r (l n D Ν)}}$ (6)

其中,E(lnDN)和Var(lnDN)可以通过图像样本估计如下:

$\hat{E} (\ln D Ν) = \frac{1}{Μ} \sum_{i = 1}^{Μ} \ln D Ν_{i}$ (7)

$V \hat{a} r (\ln D Ν) = \frac{1}{Μ - 1} \sum_{i = 1}^{Μ} [\ln D Ν_{i} - \hat{E} (\ln D Ν)]^{2} (8)$

其中,DNi代表图像样本,M是样本个数。显然,由式(5)和式(6)可知,只要知道散射系数的统计分布并且获得它的对数的均值和方差,就可以方便地由观察到的DN值估计出校正因子。

2 基于Rayleigh模型的校正因子估计

散射系数最经典的统计模型是Rayleigh模型[1,2]。当地表分辨单元内存在大量散射体,同时又没有占主导地位的散射体,此时,Rayleigh模型可以精确描述雷达回波的统计特性。对于线性检波器,Rayleigh模型的概率密度函数如下[2,4,5]:

$f (x) = \frac{x}{2 γ} e^{- \frac{x^{2}}{4 γ}}$ (9)

其中,γ(γ>0)是尺度参数。图1(a)给出了线性检波器下Rayleigh模型的概率密度函数。显然,此时的Rayleigh模型具有明显的单峰特性。假定X为服从线性检波器下Rayleigh模型的随机变量,则Rayleigh模型的对数X′(X′=lnX)的概率密度函数可以写为如下形式:

$f_{X^{'}} (x) = \frac{e^{2 x}}{2 γ} e^{- \frac{e^{2 x}}{4 γ}}$ (10)

本文使用Rayleigh模型描述散射系数,由此,可以获得散射系数Rayleigh模型的对数的均值和方差如下: $E (X^{'}) = E (\ln σ^{0}) = - \frac{C_{e}}{2} + \ln 2 + \frac{\ln γ}{2}$ (11)

$V a r (X^{'}) = V a r (\ln σ^{0}) = \frac{π^{2}}{24}$ (12)

其中,Ce是欧拉常数(Ce≈0.5772)。把式(11)和式(12)代入到式(5)和式(6)中,经过适当化简,可得线性检波器下校正因子的估计式如下:

$\hat{A} = \frac{10}{\ln 10} [- \frac{C_{e}}{2} + \ln 2 + \frac{\ln γ}{2} - \frac{π \hat{E} (l n D Ν)}{\sqrt{24 V \hat{a} r (l n D Ν)}}] (13)$

$\hat{B} = \frac{10 π}{\sqrt{24 V \hat{a} r (l n D Ν)}}$ (14)

其中, $\hat{E} (\ln D Ν)$ 和V $\hat{a}$ r(lnDN)可以通过式(7)和式(8)由图像样本估计得到。

对于平方率检波器,Rayleigh模型的概率密度函数可以写为如下形式[2,4]:

$f (x) = \frac{1}{4 γ} e^{- \frac{x}{4 γ}}$ (15)

显然,此时的Rayleigh模型实际上就是传统的指数分布。图1(b)给出了平方率检波器下Rayleigh模型的概率密度函数,可见,此时的Rayleigh模型已经不再具有单峰特性。本文同样使用Rayleigh模型描述散射系数,类似地,可以获得散射系数Rayleigh模型的对数的均值和方差如下:

E(X′)=E(lnσ0)=-Ce+2ln2+lnγ (16)

$V a r (X^{'}) = V a r (\ln σ^{0}) = \frac{π^{2}}{6}$ (17)

其中,X′=lnX,X是服从平方率检波器下Rayleigh模型的随机变量。把式(16)和式(17)代入到式(5)和式(6)中,经过适当化简,可得平方率检波器下校正因子的估计式如下:

$\hat{A} = \frac{10}{\ln 10} [- C_{e} + 2 \ln 2 + \ln γ - \frac{π \hat{E} (l n D Ν)}{\sqrt{6 V \hat{a} r (l n D Ν)}}]$ (18)

$\hat{B} = \frac{10 π}{\sqrt{6 V \hat{a} r (l n D Ν)}}$ (19)

3 Monte Carlo仿真

Monte Carlo仿真的一个重要问题就是如何产生服从特定分布的样本。本文使用Rayleigh模型描述散射系数,并通过文献[6]所提方法产生服从Rayleigh模型的样本,然后通过式(1)获得DN的样本值。对于一组DN样本值,使用本文所提方法估计出校正因子,这就完成了一次仿真实验。由于样本产生的随机性,往往需要独立运行多次仿真实验,并把多次仿真结果的均值和标准差作为最终的估计结果,这就是Monte Carlo仿真的一般步骤。总之,基于散射系数Rayleigh模型的校正因子估计的软件编程实现方法可以归纳为以下5步:(1) 产生服从Rayleigh模型的样本;(2) 生成DN样本值;(3) 根据式(7)和式(8)计算DN的对数的样本均值和样本方差;(4) 根据式(13)和式(14)(对于线性检波器)或者式(18)和式(19)(对于平方率检波器),获得校正因子的一次估计值;(5) 重复Monte Carlo仿真,计算多次估计值的平均值和标准差,获得校正因子的最终估计结果。采用MATLAB语言格式,基于散射系数Rayleigh模型的校正因子估计的核心编程可以描述如下(N为Monte Carlo仿真次数)。

使用MATLAB 6.5软件平台,在Windows XP操作系统和Intel 1.86GHz主频以及256MB内存的硬件环境下进行校正因子估计实验。为了精确评价所提方法的估计性能,本文采用仿真数据进行实验,即先仿真产生服从Rayleigh模型的样本,然后再生成DN的样本值。为了克服样本产生的随机性这一缺陷,本文采用Monte Carlo仿真方法,即使用多次仿真结果的平均值作为最终的估计值。校正因子估计首先要通过式(7)和式(8)计算DN样本对数的均值和方差的估计值,因此,样本数量是影响校正因子估计性能的重要因素。图2给出了线性检波器下校正因子估计的性能随样本数量的变化关系。可见,校正因子的估计值均随着样本数量的增加而逐渐趋向真实值,其估计标准差均随着样本数量的增加而逐渐降低,即较多的样本数量可以获得较高的估计性能。为了克服样本产生的随机性,Monte Carlo仿真需要独立运行多次,并把多次仿真结果的均值作为最终的参数估计值。因此,运行次数是影响校正因子估计精度的又一因素。图3给出了线性检波器下校正因子的估计精度随运行次数的变化关系。可见,校正因子的估计值总体上均随着运行次数的增加而趋向真实值,即较多的运行次数可以获得较高的估计精度。在Rayleigh模型参数和校正因子分别取不同值时,对应于线性检波器和平方率检波器,校正因子估计的Monte Carlo仿真结果分别如表1、表2和表3所示。其中,括号中的数值为多次仿真结果的标准差。可见,不论Rayleigh模型参数和校正因子真值取为何值,也不论检波器属于何种类型,本文提出的校正因子估计均能获得较高的估计精度。

(A=10,B=10,γ=1,对于不同的样本数量,独立运行次数均为100)

(A=10,B=10,γ=1,对于不同的运行次数,样本数量均为10000)

仿真结果见表1-表3。

4 结语

由于实际成像过程中误差的存在,图像的DN值一般不能直接反映地物表面的散射系数,而需要使用校正方程对图像进行校正,其关键在于如何高效估计校正方程中的校正因子。基于散射系数的Rayleigh模型,本文提出了一种估计SAR图像校正因子的简单方法。只要知道散射系数的统计分布,就可以依据DN观察值方便地估计出校正因子。本文使用经典的Rayleigh模型描述散射系数的统计特性,由此获得了解析形式的校正因子估计式。Monte Carlo仿真表明,基于散射系数Rayleigh模型的校正因子估计可以获得较高的估计精度,因此它是校正因子估计的有效方法。使用本文所提方法,可以高效估计出校正方程的校正因子,从而通过图像的DN值精确获得地物表面的散射系数,最终获得不同地物目标的电磁特性和图像灰度值的定量关系,为图像的精确解译打下坚实的基础,这就是本文所提方法的应用前景。

摘要：合成孔径雷达(SAR)图像的校正因子一般通过雷达系统参数的估计而获得。基于散射系数的Rayleigh模型,提出一种估计校正因子的简单而有效的方法。分别对线性形式的校正方程取均值和方差,并使用Rayleigh模型描述散射系数,针对线性检波器和平方率检波器,获得校正因子估计的解析表达式。Monte Carlo仿真表明基于散射系数Rayleigh模型的校正因子估计具有较高的估计精度。

关键词：合成孔径雷达图像,散射系数,Rayleigh模型,校正因子,Monte Carlo仿真

参考文献

[1]Ulaby F T,Dobson M C.Handbook of radar scattering statistics for ter-rain[M].Dedham:Artech House,1989.

[2] Skolnik M I.Introduction to radar systems[M].3rd ed.New York: McGraw-Hill,2001.

[3] Curlander J C,Mcdonough R N.Synthetic aperture radar:systems and signal processing[M].New York:John Wiley & Sons,1991.

[4] Oliver C,Quegan S.Understanding synthetic aperture radar images[M].Boston:Artech House,1998.

[5] Kuruoglu E E,Zerubia J.Modeling SAR images with a generalization of the Rayleigh distribution[J].IEEE Transactions on Image Processing,2004,13(4):527-533.

系数估计篇2

应用支持向量回归估计预测陀螺误差系数

针对目前小样本容量陀螺误差系数预测精度不高的.问题,本文将支持向量回归估计引入到陀螺误差系数的预测研究中.通过对某型陀螺某项误差系数的预测,并且对比分析该方法与目前通用的AR模型预测方法的预测效果,结果表明本文采用的支持向量回归估计具有更高的预测精度.

作者：焦巍王宏力刘光斌 JIAO Wei WANG Hong-li LIU Guang-bin 作者单位：第二炮兵工程学院,西安,710025 刊名：电光与控制 ISTIC PKU英文刊名：ELECTRONICS OPTICS & CONTROL 年，卷(期)： 13(5) 分类号：V241.5 关键词：支持向量机支持向量回归估计误差系数预测 AR模型

系数估计篇3

1. 引言

考虑线性回归模型: y = Xβ + e,E( e) = 0,Cov( e) =σ2In. ( 1. 1)

其中,为n×1随机观测向量,X为n×p的设计矩阵且已中心化和标准化,rank( X) = p,β为p×1的未知参数向量,e为n×1随机误差向量,In为n阶单位矩阵. 存在p×p正交矩阵Φ使得X'X = ΦΛΦ',其中Λ = diag( λ1,λ2,…, λp) ,Φ = ( φ1,φ2,…,φp) ,这里λ1≥λ2≥…≥λp为X'X的特征值,φ1,φ2,…,φp为对应的标准正交化特征向量. 令: Z = XΦ,α = Φ'β,则得到线性回归模型( 1. 1) 的典则形式:

y = Zα + e,E( e) = 0,Cov( e) = σ2In.( 1. 2)

α的LS估计为: 从而原始参数β的LS估计为当设计矩阵X的列向量之间出现复共线关系时,称模型( 1. 1) 为病态线性回归模型. 由文献可知,当且仅当X'X存在很小特征值时模型出现病态. 不妨假设非常大. 因而在均方误差意义下LS估计不再是一个好的估计. 为了解决这个问题,统计学家们做了大量工作. 目前应用最为广泛的有两种方法: 一是Hoerl和Kennard于1970年提出的岭估计. 对模型( 1. 1) ,回归系数β的岭估计定义为: 二是主成分估计. 岭估计是以牺牲无偏性换取方差部分的大幅度减小,达到最终降低其均方误差的目的. 但从上述分析可知,真正使得LS估计变坏的原因在于λr +1,…,λp很小, 因而增大λr +1,…,λp是有必要的. 对于主成分估计,虽然后面p - r个主成分对因变量影响较小,但毕竟是影响y的一些因素,若轻易剔除,显然有失真之弊. 为了弥补这些不足,本文提出主成分—岭估计的设想.

2. 主成分 — 岭估计的定义

( 1 - c) y = Z 2 α 2 +e/ 2. ( 2. 2)

其中E( e) = 0,Cov( e) = σ2In,c的取值可根据实际需要预先取定. 易见,以上几个模型本质是相同的.

定义: 在模型( 2. 2) 中,回归系数α =(α1,α 2 )T的主成分—岭估计估计定义为:

3. 主成分 — 岭估计的基本性质

引理1 α1,α2的估计α*1,α*2( k) 具有下列性质: ( 1) c'α* 1是c'α1的最佳线性无偏估计; ( 2) α*2( k) 是α2的一个有偏估计;

引理2 β*具有以下基本性质: ( 1) β*是最小二乘估计向原点的一种压缩,且存在k > 0,使得β*是岭估计向原点的一种压缩; ( 2) β*比岭估计具有更小的偏差.

系数估计篇4

考虑推广增长曲线模型中回归系数的线性估计问题和在线性估计类中的可容许性,利用增长曲线模型结果,得到了线性估计在线性估计类中可容许的两个充分条件.

作者：焦万堂吴志德李俊海 JIAO Wan-tang WU Zhi-de LI Jun-hai 作者单位：焦万堂,李俊海,JIAO Wan-tang,LI Jun-hai(河南工业大学理学院,郑州,450001)

吴志德,WU Zhi-de(郑州大学数学系,郑州,450001)

系数估计篇5

对于一个城市而言,其科技力量的强弱、科技进步能力的大小,将成为决定城市可持续发展的主要因素。资本产出弹性系数和劳动力产出弹性系数是定量测算科技进步对经济增长贡献率的重要环节。

1 资本与劳动的产出弹性系数的概念

资本的产出弹性系数(用α表示),指在其他条件不变的情况下,资金增加1%时,产出增加(α%),一般为0.3-0.35。

劳动的产出弹性系数(用β表示),指在其他条件不变的情况下,劳动增加1%时,产出增加β%,一般为0.7-0.65。

科技进步贡献率测算模型采用国家计委、国家统计局推荐的增长速度方程:

Ea=(a/y)×100% (1)

a=y-αK+βL (2)

Ea为科技进步对产值增长速度的贡献,即在产值增长速度中的科技进步因素所占的比重;a为科技进步的年平均增长速度;(1)式中,y为产出的年平均增长速度;资金是固定资产和流动资金之和;K为资金的年平均增长速度;L为劳动者的年平均增长速度;α为资本的产出弹性系数;β为劳动的产出弹性系数。(2)式简化为;技术进步速度=国民收入的增加率-资本对国民收入的贡献率×资本的增加率+劳动对国民收放的贡献率×劳动力的增加率。

2 资本与劳动产出弹性系数的估计

2.1 份额估计法

国民收入的使用分为:积累和消费。粗略地讲:积累主要用于资金投入;消费主要用于劳动投入。因此,就全国或省、市、区大范围而言,取积累率为α消费率为β具体到一个企业,固定资产折旧加流动资金占总成本的比率为α;工资加福利开支占总成本的比率为β。即:

α=积累率=(固定资产折旧+流动资金)/总成本 (3)

β=消费率=(工资总额+福利开支)/总成本 (4)

以上主要针对工业而言,就我国而言,可以取α=0.25,β=0.75。

尽管这种方法很粗糙,但它简便易行。也可以资本形成率和最终消费来代替(以支出法计GNP时的统计指标,在一般年鉴中有)。

西方有些国家根据国民收入Y中,工资V和利润M所占的份额,估计α及β的值。即:

α=利税总额/国民收入 (5)

β=工资总额/国民收入 (6)

2.2 回归法

运用数学模型测算技术进步因子贡献的研究需要引入生产函数。不同经济学流派对生产函数的理解不同,因此生产函数的表达形式也各不相同。如果用Q表示产量,用X1,X2……Xn表示不同要素投入量,则可以得到生产函数的一个一般表达式:

Q=F(X1,X2…Xn)

这个一般形式,说明投入量和产出之间存在的映射关系。目前理论界认为比较成熟的是柯布—道格拉斯(Cobb-Douglas)生产函数为对象,以线性回归分析来对资本与劳动的产出弹性加以定量研究, C—D生产函数是美国经济学教授道格拉斯((P-H-Douglas)与数学家柯布(C-W-Cobb)合作,研究分析了美国制造业1899-1922年的历史资料后得出的,1942年首届诺贝尔经济学奖获得者丁伯根(J-Tinbergen )对此方法作了改进,后称为C-D生产函数[1]。C-D函数的表达式为:

Q=AKαLβ (7)

式中:

α+β=1 0<α<1

对上式两边分别求自然对数,有:

lnQ=lnA+αlnK+βlnL (8)

其中α、β即为资本与劳动的产出弹性系数即边际报酬,Q为国内生产总植(GNP),K为资本量,L为劳动投入量。(8)式中Q、K、L为历年的采集得到的数据,而A为常数。这样,α、β的估计就符合一个线性回归模型[2]。

线性回归方法与份额分析方法是估计参数α与β的常用方法。线性回归方法需要较多年份的历史数据,因此,本文选用13年以上的数据资料进行分析。王新雷、张培基[3]等在《电力工业科技进步评价方法研究》中,参考了国内外大量经济学资料和国内其他部委测算方法中的α和β,最后,确定整个电力工业的α=0.3,β= 0.70。

3 资本与劳动产出弹性系数估计的软件实现

3.1 语言选择

本文选择了Visual Basic6.0面向对象的程序设计语言,因为Visual Basic6.0可以方便的访问数据库,同时在数据计算方面也有着强大的功能。数据库选择了Access数据库,因为科技贡献率的测算的数据量不大,一年只有一条记录。

3.2 数据库的建立

采用线性回归分析的方法来确定α、β的值时,技术进步对常数项,数据资料需求同城市劳动力结构及历年的国内生产总值(GNP)统计表,如表1所示。

因此,在数据库sciadvance.mdb里建立一个数据表advan。表结构如图1所示。其中gnp、money、labor三个字段设置为单精度数字。

3.3 界面设计

用于实现“线性回归分析的方法来确定α、β的值”的程序,共有三个窗体,包括数据输入窗体、数据选择窗体、数据计算窗体,是整个“区域科技发展状况评价计量”软件的一部分。本文重点介绍两个窗体。数据录入窗体如图2所示;计算α、β值的窗体,如图3(b)所示,选择计算用数据的窗体如图3(a)所示,为保证计算出的α、β的客观性,数据的选择要求13年以上[4]。

3.4 算法实现

本文省略数据选择窗体的程序设计,重点介绍计算生成α、β值的方法。

由式(8)可知,式中的Q对应于数据库中“gnp”字段,K对应于数据库中的“money”字段,L对应于数据库中的“labor”字段,例用线性回归法确定α、β和lnA的值。

设lnQ=Y, lnK=X1,lnL=X2,lnA=b0,α=b1,β=b2则,式(8)可以写作:

Y=b0 +b1X1 +b2X2

(Y1,X11,X21)

(Y2,X12,X22)

……

(Yn,X1n,X2n)

是一个样本,其中n≥13。采用最小二乘法[5]估计b0、b1、 b2的值。引入偏差平方和:

$G (b_{0}, b_{1}, b_{2}) = \sum_{k = 1}^{n} (y_{k} - b_{0} - b_{1} X_{1 k} - b_{2} X_{2 k})^{2}$

最小二乘估计就是求B0,B1,B2,使得:

G(B0,B1,B2)=minG(b0,b1,b2)

因为G(b0,b1,b2)是b0,b1,b2的非负二次型,故其最小值一定存在。分别求G对于b0,b1,b2偏导数,并令其为0。

$\begin{array}{l} \frac{\partial G}{\partial b_{0}} = \sum_{k = 1}^{n} (- 2) (y_{k} — b_{0} — b_{1} X_{1 k} — b_{2} X_{2 k}) = 0 \\ \frac{\partial G}{\partial b_{1}} = \sum_{k = 1}^{n} (- 2 X_{1 k}) (y_{k} — b_{0} — b_{1} X_{1 k} — b_{2} X_{2 k}) = 0 \\ \frac{\partial G}{\partial b_{2}} = \sum_{k = 1}^{n} (- 2 X_{2 k}) (y_{k} — b_{0} — b_{1} X_{1 k} — b_{2} X_{2 k}) = 0 \end{array}} (9)$

上述方程组称为正规方程组,其等价形式为:

$\begin{array}{l} \sum_{k = 1}^{n} Y_{k} = n b_{0} + b_{1} \sum_{k = 1}^{n} X_{1 k} + b_{2} \sum_{k = 1}^{n} X_{2 k} \\ \sum_{k = 1}^{n} Y_{k} X_{1 k} = b_{0} \sum_{k = 1}^{n} X_{1 k} + b_{1} \sum_{k = 1}^{n} X_{1 k}^{2} + b_{2} \sum_{k = 1}^{n} X_{1 k} X_{2 k} \\ \sum_{k = 1}^{n} Y_{k} X_{2 k} = b_{0} \sum_{k = 1}^{n} X_{2 k} + b_{1} \sum_{k = 1}^{n} X_{1 k} X_{2 k} + b_{2} \sum_{k = 1}^{n} X_{2 k}^{2} \end{array}} (10)$

将数据代入(10)式,可以得到本模型的正规方程组,再对这一方程组求解,即可得到b0、b1、b2即lnA、α、β的值。

程序流程如图4所示。

3.5 软件的有效性及正确性检验

以某市的数据作为样本进行实证研究,运算结果如图5所示。由图中可以看出,α=0.8807,β=0.3214,α+β>1,说明该市的生产函数属于规模报酬递增型,且资本的产出弹性系数较之于劳动的产出弹性系数大,说明该市由于增加了产业结构调整,促进了资本密集型产业的发展。用SPSS对本组数据进行线性回归分析,得到的结论与程序得到的结论相同。从而证明计算结果无误。

摘要：运用Visual Basic与Access数据库相结合的方法,使用线性回归法估计资本与劳动的产出弹性系数。程序实现了将数据库中的数据(GNP、资本量、劳动力)使用线性回归法计算资本与劳动的产出弹性系数的功能,并将计算结果输出到一个图片框中。

关键词：Visual Basic,Access,线性回归,资本产出弹性系数,劳动产出弹性系数

参考文献

[1]科技进步对经济增长贡献的测算方法研究[J].http://www.hnst.gov.cn/tongji/show3.asp?InfoID=5.

[2]陈魁.应用概率统计[M].清华大学出版社2006.

[3]王新雷,张培基,毛晋,李彦斌.电力工业科技进步评价方法研究.科技管理研究[J].1998,2.

[4]求是科技.Visual Basic 6.0程序设计与开发技术大全[M].人民邮电出版社,2004.

系数估计篇6

在多尺度变换域, 自然图像子带系数的统计建模是许多图像处理和低层计算机视觉应用领域的基础, 例如图像去噪、图像压缩编码和纹理图像检索等[1,2,3]。其中, 广泛应用的模型是变换域子带系数的广义高斯分布。研究表明, 应用广义高斯分布的建模是对自然图像变换子带系数最逼近和最成功的建模方式[3,4]。

Contourlet变换由Do和Vetterli在2002年提出, 是在离散域用滤波器组来定义和实现的[5]。它将图像的多尺度和多方向表示灵活而又有机地结合起来, 因而能准确、最优地刻画图像。目前, 基于Contourlet变换的理论和应用是研究的热点。由于Contourlet变换能够最稀疏地刻画图像的特征, 而其优良的统计特征可以代表图像的内容, 同时, 自然图像多重子带的边缘分布能充分地反映其特征。因此, 在Contourlet变换域进行广义高斯分布统计建模有着重要的意义。

子带系数的广义高斯统计建模已经在小波域得到了广泛的采用[3,4,6]。这种模型在实际使用中存在的主要局限是它的形状参数的准确估计问题。广义高斯模型形状参数估计早期的研究工作由Varanasi 和Aazhang进行了总结[7]。典型的估计方法有三种, 即矩估计法、熵匹配法和最大似然 (ML) 估计法。

本文提出一种混合的形状参数估计方法:当形状参数值较小时, 对于小样本采用熵匹配方法估计;而对于大样本利用最大似然方法估计;当形状参数值较大时, 采用矩方法估计。

1相关理论

1.1Contourlet变换

Contourlet变换通过滤波器组来实现多尺度变换和多方向变换[5]。首先, 多尺度变换是由拉普拉斯金字塔 (LP) 实现。然后, 在LP不同尺度下的各个差值子带上由方向滤波器组 (DFB) 进行不同数目的方向变换。

Contourlet变换通过多尺度表示和多方向表示的有机结合来有效地刻画图像。多尺度表示的目的是使图像的能量集中于奇异点, 而在此基础上的多方向表示是使能量集中于奇异线段 (或边缘轮廓线段) , 这样就实现了图像的稀疏表示。具体地, Contourlet变换由金字塔方向滤波器组结构实现, 如图1 (a) 所示。图中只画出了分解过程, 重建过程与之相反。多尺度分解由LP完成。它首先将输入图像分解为一个低通的概貌子带和一个带通子带 (或差值子带) , 然后, 多尺度分解过程可在概貌子带上迭代进行。方向分解是利用DFB只在LP的各带通子带上进行, 这是由于方向滤波器组本身不适合处理图像的低频部分。图1 (b) 给出了一个Contourlet分解的例子, 图中, 对Lena图像进行了2层LP分解。其中, 保持概貌子带不变, 而对两层带通子带由高频到低频相应地分别分解了8个和4个方向子带。

1.2广义高斯分布模型

随机变量X的广义高斯分布概率密度函数定义为[6,8]

p (x; $μ, σ^{2}, β) = \frac{β η (σ, β)}{2 Γ (1 / β)} \exp {- [η (σ, β) | x - μ |]}^{β}$ (1)

式中, μ、σ2和β分别表示分布的均值、方差和形状参数。Γ (·) 为Gamma函数, 其表达式为:

$Γ (x) = \int_{0}^{\infty} t^{x - 1} e^{- t} d t x > 0$

此外, η (σ, β) 的表达式为:

$η (σ, β) = \frac{1}{σ} [\frac{Γ (3 / β)}{Γ (1 / β)}]^{1 / 2}$ (2)

式 (1) 中, 当β=1时对应于拉普拉斯分布;β=2时对应于高斯分布;而当β→∞时对应于均匀分布。参数β的大小决定广义高斯分布概率密度曲线的形状。β的值越小, 分布曲线越尖锐;反之, 分布曲线具有较大的拖尾。图2示出了随不同的β值相应的广义高斯分布概率密度曲线, 图中, 分布的均值为零, σ=6。

2Contourlet变换域分布参数的估计

在多尺度变换域, 给定子带图像数据, 广义高斯分布的三个参数, 即均值μ、方差σ2和形状参数β需要估计。通常, 对于子带样本值xi, i=1, 2, …, N;当子带样本总数N足够大时, 均值μ、方差σ2的估计可以采用下式进行:

$\hat{μ}_{X} = \frac{1}{Ν} \sum_{i = 1}^{Ν} x_{i} \hat{σ}_{X}^{2} = \frac{1}{Ν} \sum_{i = 1}^{Ν} (x_{i} - \hat{μ}_{X})^{2}$ (3)

因此, 下面的研究主要集中在形状参数β的估计上。

2.1矩估计方法

在多尺度变换域, 形状参数β的矩估计方法由Mallat提出[4]。对于广义高斯分布的随机变量X, 有:

$\begin{array}{l} E [| x - μ |^{2 m}] = \int_{- \infty}^{\infty} | x - μ |^{2 m} \frac{β η (σ, β)}{2 Γ (1 / β)} e^{- [η (σ, β) | x - μ |]^{β}} d x \\ = σ^{2 m} \frac{Γ ((2 m + 1) / β) Γ^{m - 1} (1 / β)}{Γ^{m} (3 / β)} β > 0 m > 0 (4) \end{array}$

根据式 (4) , 所采用的估计形状参数β的方法描述如下。

取: $m = \frac{1}{2} E [| x - μ |] = σ [\frac{Γ^{2} (2 / β)}{Γ (1 / β) Γ (3 / β)}]^{1 / 2}$

定义比值函数:

$F_{1} (β) = \frac{E [| x - μ |]}{σ} = [\frac{Γ^{2} (2 / β)}{Γ (1 / β) Γ (3 / β)}]^{1 / 2}$ (5)

通常, $\hat{E} [| x - μ |] = \frac{1}{Ν} \sum_{i = 1}^{Ν} | x_{i} - \hat{μ}_{x} |$ 和 $\hat{σ}^{2} = \hat{σ}_{x}^{2}$ (参见式 (3) ) 。则形状参数β根据下式利用查询表进行估计

$\hat{β} = \hat{β}_{1} = F_{1}^{- 1} (\frac{\hat{E} [| x - μ |]}{\hat{σ}})$ (6)

2.2最大似然估计方法

给定N个样本X= (X1, X2, …, XN) , 假定其服从广义高斯分布, 均值为μ, 并且样本分量是独立同分布的。则联合概率分布的对数似然函数定义为:

L (x; $μ, σ^{2}, β) = \sum_{i = 1}^{Ν} \log p (x_{i}$ ;μ, σ2, β) (7)

可以采用最大似然的方法估计σ2和β, 即利用式 (7) 分别对σ和β进行求导, 并令导数为零, 则有[3]:

$\hat{σ} = (\frac{Γ (3 / \hat{β})}{Γ (1 / \hat{β})})^{1 / 2} (\frac{\hat{β}}{Ν} \sum_{i = 1}^{Ν} | x_{i} - \hat{μ} |^{\hat{β}})^{1 / \hat{β}}$ (8)

$\begin{array}{l} g \hat{β} = 1 + \frac{Ψ (1 / \hat{β})}{\hat{β}} - \frac{\sum_{i = 1}^{Ν} | x_{i} - \hat{μ} |^{\hat{β}} \log | x_{i} - \hat{μ} |}{\sum_{i = 1}^{Ν} | x_{i} - \hat{μ} |^{\hat{β}}} + \\ \frac{\log (\frac{\hat{β}}{Ν} \sum_{i = 1}^{Ν} | x_{i} - \hat{μ} |^{\hat{β}})}{\hat{β}} = 0 (9) \end{array}$

式中, digamma函数:

$Ψ (x) = \frac{d l o g (Γ (x))}{d x} = - γ - \frac{1}{x} + \sum_{k = 1}^{\infty} (\frac{1}{k} - \frac{1}{k + x}) ‚ γ = 0.577 \dots$ 为尤拉常数。由式 (8) 和式 (9) 可见, 当 $\hat{μ} = \hat{μ}_{x}$ 由式 (3) 确定后, $\hat{β}$ 可由式 (9) 确定, 从而 $\hat{σ}^{2}$ 可由式 (8) 确定。

由于式 (9) 为一超越方程, 因此 $\hat{β}$ 的求解可采用Newton-Raphson迭代算法实现。具体的算法如下:

β1=β0

$\begin{array}{l} W h i l e (g (β_{n} + ε) > 0)) \\ β_{n + 1} = β_{n} - \frac{g (β_{n})}{g^{'} (β_{n})} (10) \end{array}$

End while

$\hat{β} = β_{n + 1}$

算法中, β0为β的初始值, 它根据E[|x-μ|]/σ, 其值由式 (5) 确定。ε为控制精度。

2.3熵匹配方法

熵匹配方法特别适合于在小样本情况下对形状参数的估计[6]。它是利用微分熵和离散熵的关系而进行估计的。对于概率密度为p (x) 的随机变量X, 其微分熵h (x) 定义为[9]:

$h (x) = - \int_{- \infty}^{\infty} p (x) \log_{2} p (x) d x$ (11)

而对于广义高斯分布, 微分熵可以表达为[6]:

$h_{G G} (σ, β) = - \int_{- \infty}^{\infty} p (x$ ;μ, σ2, β) log2p (x;μ, σ2, β) dx

$= - \log_{2} [\frac{β η (σ, β)}{2 Γ (1 / β)}] + \frac{1}{β \ln 2}$ (12)

令离散熵为 $Η (Δ) = - \sum_{i} p_{i} \log_{2} p_{i}$ , 其中, Δ为均匀量化的步长。对于给定样本数据, 首先按步长Δ均匀量化, 则pi为量化后的整数数据分布的概率密度。根据文献[9]中的定理9.3.1, 有:

hGG (σ, β) ≅H (Δ) +log2Δ (13)

将式 (2) 和式 (12) 代入式 (13) , 可得:

$F_{Η} (β) = Η (Δ) - \log_{2} \frac{σ}{Δ} = - \log_{2} [\frac{β Γ^{1 / 2} (3 / β)}{2 Γ^{(3 / 2} (1 / β)}] + \frac{1}{β \ln 2}$ (14)

令 $\hat{σ} = \hat{σ}_{X}$ (参见式 (3) ) 。H (Δ) 和Δ可根据样本数据的直方图求得。这样, 形状参数β根据下式利用查询表进行估计:

$\hat{β} = F_{Η}^{- 1} (Η (Δ) - \log_{2} \frac{\hat{σ}}{Δ})$ (15)

应当注意:对于式 (14) 中的FH (β) , 当β=2时, 它达到全局最大值, 为 $\log_{2} \sqrt{2 π e}$ 。由此可见, 熵匹配方法不适合估计形状参数β为2周围的值。

2.4基于Contourlet子带系数的形状参数估计方法

本文采用上述三种方法进行Contourlet变换子带系数的混合建模。根据Contourlet变换子带的大小N, 样本数据可分为小样本 (N～O (103) ) 和大样本N～O (105) 两种[6]。通常, 熵匹配估计方法适合于小样本情况。相比于矩估计方法, 它在形状参数β≤1时的估计效果最好。图3给出了随β的变化矩估计函数F1 (β) 和熵匹配函数FH (β) 的关系曲线, 当β≤1时, FH (β) 的斜率大于F1 (β) 。因此, 采用熵匹配方法对形状参数β的估计误差小于矩估计方法。最大似然估计方法也适合于β≤1的情况[7]。但是, 为了保证具有良好的逼近性能, 它通常用于大样本时更有效。此外, 当β>1时, 由于矩估计函数曲线具有单增特性, 因此可用F1 (β) 对较大的β值进行估计。

基于上述的讨论, 对于Contourlet 变换的每一个方向子带和概貌子带, 首先用矩估计方法进行形状参数的估计。若其值大于1, 则估计的结果保留;否则, 对于小样本, 采用熵匹配方法估计, 而对于大样本, 利用最大似然方法估计。

3实验结果

与基于小波变换的广义高斯分布建模方法相似, 本实验的建模方法也是将分布密度曲线和Contourlet子带系数的直方图相拟合[3]。在Contourlet变换实验中, 其LP分解选择 “9-7”双正交滤波器, 而DFB分解采用具有梯型结构的双正交、支撑域为 (23, 23) 和 (45, 45) 的PKVA滤波器[5]。对于图像Barbara (大小为512×512) , Contourlet变换进行三层 LP分解。保持概貌子带不变, 对于带通子带, DFB分解的方向子带数由细尺度到粗尺度分别为 (8, 4, 4) 。建模实验的结果如图4所示, 图中只列出了对三个典型的方向子带和概貌子带的建模结果。由图可见, 按本文提出的算法所估计的广义高斯分布概率密度曲线与实际的方向子带和概貌子带样本数据的分布得到了很好的拟合。

4结论

在本文中, 基于Contourlet变换子带系数的广义高斯分布建模, 提出了一种新的混合形状参数估计方法。对于形状参数β∈ (0, 1], 在小样本情况下采用熵匹配方法估计, 在大样本情况下利用最大似然方法估计。而当β大于1时, 矩估计方法被利用。实验结果表明, 本文所提出的方法可以有效地和准确地对Contourlet变换子带系数进行建模。

摘要：提出采用广义高斯概率密度建模的Contourlet变换系数的形状参数混合估计方法。当形状参数值较小时, 对于小样本采用熵匹配方法估计, 而对于大样本利用最大似然方法估计。当形状参数值较大时, 采用矩方法估计。实验结果表明, 所提出的方案可以有效准确地对Contourlet变换子带进行建模。

关键词：Contourlet变换,系数建模,广义高斯分布

参考文献

[1]Kivanc Mihcak M, Kozintsev I, Ramchandran K, et, al.Low-complexi-ty image denoising based on statistical modeling of wavelet coefficients[J].IEEE Signal processing Letters, 1999, 6 (12) :300-303.

[2]LoPresto S M, Ramchandran K, Orchard M T.Image coding based on mixture modeling of wavelet coefficients and a fast estimation-quantiza-tion framework[C]//Proc.IEEE Data Compression Conference, Snow-bird, UT, 1997:221-230.

[3] DO M N, Vetterli M.Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[J].IEEE Transactions on Image processing, 2002, 11 (2) :146-158.

[4] Mallat S G.A theory for multiresolution signal decomposition: the wavelet representation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence , 1989, 11 (7) :674-693.

[5] DO M N, Vetterli M.The contourlet transform: an efficient directional multiresolution image representation[J].IEEE Transactions on Image Processing, 2005, 14 (12) :2091-2106.

[6]Aiazzi B, Alparone L, Baronti S.Estimation based on entropy matching for generalized Gaussian pdf modeling[J].IEEE Signal Processing Let-ters, 1999, 6 (6) :138-140.

[7]Varanasi MK, Aazhang B.Parametric generalized Gaussian density es-timation[J].J.Acoust.Soc.Amer.1989, 86 (4) :1404-1415.

[8] Sharifi K, Leon-Garcia A.Estimation of shape parameter for generalized Gaussian distributions in subband decompositions of video[J].IEEE Transactions on Circuits and Systems for Video Technology, 1995, 5 (1) :52-56.

系数估计篇7

关键词：粒子群优化算法,洛伦兹曲线,基尼系数

居民收入分配关系到广大人民群众的生活水平, 基尼系数(Gini coefficient)是衡量一个经济体收入分配状况重要经济指标,洛伦兹(Lorenz curve)曲线是用来刻画收入分配不平等程度

经济指标。图1中所示的直线y =x为绝对平等线,曲线即为洛伦兹曲线。

洛伦兹曲线与基尼系数之间的关系式为:

其中SA为绝对平均线与洛伦兹曲线间的几何面积,SB为洛伦兹曲线下方的曲三角形面积,对SB的计算大致可分为离散法和连续法两种[1],其中连续法是用一个函数L(p)来拟合洛伦兹曲线,则。由此基尼系数可定义为:

对于洛伦兹曲线方程,国内外学者如Bas- mann[2]、Sarabia[3]、Chotikapanich[4]已提出了各种形式的函数来拟合洛伦兹曲线,已有的这些模型存在诸多问题。例如:拟合效果欠佳,有些模型甚至不满足洛伦兹曲线的条件。本文采用张奎、王原提出的Sarabia洛伦兹曲线的推广模型[5],数值实验说明此模型拟合效果远优于前人提出的模型。而对模型中参数的估计大都采用传统的方法,如:密度函数的Kernel估计法、无约束非线性最小二乘估计法以及回归分析等。现代智能优化算法在经济、工程等各个领域中有着广泛的应用。隶属于智能优化算法的粒子群优化算法对复杂非线性问题具有较强的寻优能力以及简单通用等显著特点,此算法概念清晰,容易实现,同时又有深刻的智能背景。鉴于此,本文对模型中的参数尝试采用粒子群优化算法进行估计,进而分别确立城镇与农村的洛伦兹曲线。由该算法估计参数确定的洛伦兹曲线计算2002-2011年我国城镇和农村的基尼系数,再利用“分组加权法”[6]得出总体的基尼系数。

1洛伦兹曲线方程与基尼系数

1.1洛伦兹曲线方程

L(p,τ)是定义在[0,1]区间上、取值于[0,1] 区间上的函数,若满足以下特征

则L(p,τ)被称为洛伦兹曲线方程[7],τ为洛伦兹曲线方程中的参数。文献[8]中已提出很多种洛伦兹曲线模型,大都是一次函数、二次函数、幂函数和指数函数。本文采用张奎,王原提出的对Sarabia洛伦兹曲线的推广模型:

其中α≥0,γ≥0,β∈(0,1],η≥1,α+γ≥1均为参数。p为收入群体的累计人口比例,0≤p≤1; L为该群体拥有的总收入比例,0≤L ≤1。

1.2基尼系数

首先根据2002-2011年我国农村和城市分组数据拟合出历年农村和城市的洛伦兹曲线L1(p)和L2(p),由此计算出各年的农村和城镇的基尼系数G1和G2,再利用Sundrum提出的“分组加权法”,计算总体基尼系数。总体基尼系数计算公式为:

式(4)中m为全体居民的人均收入,m1为农村居民人均纯收入,m2为城镇居民人均总收入,P1和P2分别为当年农村和城市的人口比例。“利用此式计算总体的居民收入基尼系数时的假设条件是‘城乡居民收入分布不重叠’,但我国全国和各省级区域城镇和农村收入分布大多存在重叠现象”[9],收入分组不重叠这一难题并没有得到解决,因此本文仍用式(4)计算总体的居民收入基尼系数。

2粒子群优化算法(PSO)及算法步骤

2.1 PSO算法原理

粒子群优化(Particle Swarm Optimization)算法[10]简称PSO,算法中以每个粒子的位置表示待优化参数的解,初始时随机生成一个具有n个粒子的种群,每个粒子的维数均为d 。由待优化问题的目标函数确定的适应度值评价每个粒子的性能,每次迭代都保留适应度值较优的粒子的位置。pi(t)是粒子i本身到目前为止寻到的最优解,称为个体最优解。 pg(t)是整个粒子群迄今搜索的最优解,称为全局最优解。第t+1次迭代计算,粒子i根据下列规则来更新自己的速度和位置:

其中w为惯性权重,c1、c2为学习因子,r1、r2为0~1之间的随机数,r为约束因子。

2.2 PSO算法基本步骤

步骤1:初始化粒子种群、个体最优解以及全局最优解;

步骤2:根据式(5)和式(6),更新每个粒子的位置和速度;

步骤3:计算每个粒子的适应度值。本文取适应度函数为拟合值和实际值的均方误差(MSE,mean squared error)即:

步骤4:将每个粒子的适应度值与个体最优值进行比较,如果较优则更新当前的个体最优解;

步骤5:将每个粒子适应值与全局最优值进行比较,如果较优则更新当前的全局最优解;

步骤6:如果达到预先设定的最大迭代次数或者目标函数的精度要求,则停止计算;否则返回步骤2。

3实例分析

表1和表2是根据历年国家统计局公布的“人民生活”相关数据整理分析所得。应用PSO算法对(pi,Li)这组数据进行估计,使拟合误差达到最小时所确定的L(p)方程中的参数即为所求的最优解。 PSO算法在解决不同问题时参数设定均不同,针对本文问题经过多组参数的调试,结果表明将参数设定[11]为:c1=1.7592,c2=1.7592,w =0.725,r =0.5,种群规模50,最大迭代次数100为较佳。选取此组参数可使算法快速收敛、粒子震荡轨迹的幅度随时间不断减小,并且能够快速搜索到全局最优解, 避免了过度探索。以下以我国2002年农村居民人均纯收入和城镇居民人均总收入的分组数据为例对洛伦兹曲线中的参数进行估计,并计算其对应的基尼系数。

3.1估计参数

表3和表4中的第4、5行数据都是通过Leven -berg-Marquardt无约束最小二乘算法和PSO算法分别估计出洛伦兹曲线方程中的参数,再根据此方程计算得到L(p)的估计值。表5、表6是两算法对L(p)估计值的绝对误差比较。表5中第2、4组数据,表6中第1、2、4、6组数据的PSO比LM算法的绝对误差小。显然,对于2002年数据洛伦兹曲线方程中的参数估计,PSO比LM算法较优。利用PSO算法的估计参数所确定的2002年农村拟合洛伦兹曲线方程为:

城镇的拟合洛伦兹曲线方程为:

再根据式(1),计算农村和城镇的基尼系数分别为G1=0.3243、G2=0.4091。由已求得农村和城镇的洛伦兹曲线和基尼系数,根据式(4)和表8中2002年相关数据计算我国当年总体的基尼系数为G =0.4606。

3.2计算基尼系数

类似于计算2002年基尼系数的过程,估计2003 -2011年我国农村和城镇的洛伦兹曲线中的参数及基尼系数的计算结果如表7和表8。

应用上述方法计算的我国2003-2011年居民的基尼系数如表9中最后一列所示,我国在2000年公布的基尼系数为0.412,根据本文的估计2009年我国总体基尼系数是0.5003与0.5相差0.003,2011年的基尼系数为0.4891,距0.5也非常接近。由此可知本文所估计的基尼系数的变化趋势与大多数学者的预测相符。

4总结

系数估计篇8

说话人识别分为两个主要过程, 即特征提取和模式识别。其中特征提取是说话人识别的关键问题。目前, 在说话人识别系统中最常用的特征参数主要有Mel频率倒谱系数 (Mel-frequency ceptrum coefficient, MFCC) 、线性预测倒谱系数 (linear predictive cepstrum coefficient, LPCC) 、基因周期、感知对数面积比系数 (perceptual log area ratio, PLAR) [1]等。MFCC特征参数更符合人耳的听觉特性, 比其他特征参数具有更优的识别特性, 所以是目前使用最广的特征参数。模式识别中常用的方法有支持向量机 (support vector machine, SVM) 、矢量量化方法、隐马尔可夫模型方法、动态时间规整方法和神经网络等。

1 MFCC参数的提取

1.1 MFCC参数

因为人耳对外来的信号有不同的调节作用, 所以对于不同频率的语音, 人耳有不同的感知能力。在1 000 Hz以下, 感知能力与频率成线性关系, 而在1 000 Hz以上, 感知能力则与频率成对数关系[2]。所以, 人们提出更适合人耳听觉系统的Mel频率。Mel频率的意义的是:1 Mel为1 000 Hz的音调感知程度的1/1 000[2]。实际频率和Mel频率之间的转换关系式为:

式 (1) 中f是实际频率, 单位是Hz。

MFCC参数就是在Mel频率域中提取出来的。先是将语音的频谱转化基于Mel标尺频率的非线性频谱, 然后再将所得到的频谱转化到倒谱域上。这样将信号的频谱转化到可以被人耳感知频域中, 可以达到更好的效果。

1.2 MFCC参数的提取过程

MFCC参数的提取过程如图2所示。

MFCC特征参数的提取及计算过程如下:

(1) 预加重。对语音信号进行预加重是为了去除口唇辐射的影响, 增加语音的高频分辨率。本文是通过将语音信号通过一阶FIR高通数字滤波器来实现预加重的。数字滤波器的传递函数是

式 (2) 中α为预加重系数, 0.9<α<1.0。

(2) 加窗分帧。语音信号是随时间变化的信号, 由于发声器官的惯性运动, 可认为在一小段的时间里语音信号近似不变, 即语音具有短时平稳性[3]。所以就可以对语音信号进行加窗分帧。在标准MFCC提取的过程中采用的是汉明窗。汉明窗的窗函数如下:

(3) DFT变换。经过FFT变换后求的语音信号的频谱。转换公式为:

式 (4) 中x (n) 是输入的语音信号, N为傅里叶变换的点数。

(4) Mel频率滤波。将上述得到的语音频谱通过Mel滤波器组。Mel滤波器是一组三角滤波器。在该滤波器组中每个滤波器的中心频率在Mel频率标尺上的跨度相同。

滤波器的传递函数为:

式 (15) 中m是滤波器的个数, f (m) 是滤波器的中心频率。

(5) 离散余弦变换DCT。DCT的主要作用是对上述滤波器的输出做变换, 最终得到MFCC特征参数c (n) 。

式 (6) 中M为三角滤波器的个数。

2 MFCC参数的改进算法

说话人识别的特征参数受到各种因素的制约。如说话人身处环境的噪音干扰。对于同一说话人, 特征参数还受到说话人的情绪、年龄, 说话时间等因素的影响。这些因素导致说话人在测试阶段和训练阶段发出的音有所差异。为了提高说话人MFCC参数特征的噪声鲁棒性和时间鲁棒性本文提出一种改进的MFCC算法。该算法对传统的MFCC算法做了两方面的改进。

2.1 窗函数的改进

说话人系统中通常采用Hamming窗的离散傅里叶变换来实现语音信号的频谱估计。Hamming窗的DFT在减少了频谱泄露但增加了频谱估计的方差值, 系统的稳定性也就随之降低。本文选用一种多窗频谱估计[4—6] (multitaper spectrum estimate) 的方法来代替原始的Hamming窗求频谱的方法。

多窗频谱估计 (简称Multitapering) 是传统的加汉明窗的DFT的一种扩展, 它的本质是采用具有不同权值的多个窗函数, 用它们频域的平均值来获得信号的频谱估计。通过多窗频谱估计方法获得的信号频谱, 减少了噪声对语音频谱的影响, 而且采用该方法对信号进行频谱估计的方差值和更小。所以, 通过该方法获得的信号频谱, 能够更好的反应出声道结构, 得到更加稳定的特征参数, 简称MTMFCC。多窗频谱估计的实现框图如图3所示[7]。

本文选取的窗函数是SWCE (sinusoidal window cepstral estimator) 。SWCE具有很高的频谱估计性能, 是一种具有最优均方差的倒谱估计方法。该方法采用多个近似的正弦窗函数进行加权来实现信号的频谱估计[8]。设x (n) 为经过预处理的语音信号, 它的功率谱为Sx (f) 。则信号具有N个采样点。则其的频谱可由式 (7) 求得:

式 (7) 中hi (n) 为第i个窗函数, αi为对应的窗函数的加权值, i=0, …, K-1。

式 (8) 中, i=1, 2, …, N;N为采样点数。

式 (9) 中K为窗函数的个数。

因此, 将式 (8) 、式 (9) 和式 (10) 代入式 (7) 中, 就可以得到SWCE对信号的频谱估计。

2.2 谱包络提取的改进

MFCC参数描述的主要是可以表征说话人声道特性的谱包络, 它是在信号的短时频谱上提取的。对浊音信号而言, 其幅度谱可以看作基音频率的各次谐波对谱包络的采样, 包含了谱包络和基音频率两方面信息[9]。而基因频率都会影响其对声道特性的准确描述, 进而影响说话人识别系统的识别率。本文提出一种对MFCC提取过程中谱包络提取的改进, 即在信号通过Mel滤波器前, 先经过一个滑动平均滤波器。利用滑动平均滤波器先对信号的幅度谱进行平滑, 得到信号谱包络的近似表示。这样就可以有效抑制基音频率对MFCC特征参数的影响。把这种特征参数称为SMFCC (smoothing MFCC) 参数[10]。

2.3 新特征的提取过程

在以上的两种改进算法中, 多窗频谱估计方法能够提高信号的噪声鲁棒性, 而谱包络提取法能够提高信号的时间鲁棒性。所以将两种方法进行结合, 就可以同时提高信号噪声和时间鲁棒性。

新特征的提取是先利用SWCE窗进行频谱估计, 得到语音信号的频谱后, 再利用滑动平均滤波器对频谱进行平滑。然后再通过Mel频率滤波器组, 这样得到更加稳定的特征参数。把新的特征参数称为MTSMFCC (multitapering smoothing MFCC) 。其提取框图如图4所示。

MTSMFCC参数的具体提取过程如下:

1) 语音信号进行预加重处理。

2) 经过预加重处理的语音信号, 利用多窗SWCE进行频谱估计, 得到信号的频谱X (k) 。通过实验证明窗的个数取4~8个较好。在本文中窗个数取K=8。保存X (k) 的数据, 设U (k) =X (k) 。

3) 通过滑动平均滤波器w (n) 对U (k) 进行平滑, 得到信号的平滑包络谱Y (k) 。其中w (n) 如式 (11) 所示。

式 (11) 中DFT取256个点, 在本实验中经过多次试验N取4时效果最好。

4) 将保存的频谱X (k) 与平滑谱Y (k) 进行比较, 得到Q (k) 。令Q (k) =max[X (k) , Y (k) ]。

5) 重复步骤 (2) 、 (3) 、 (4) , 经过多次试验本实验在迭代次数取9的时候效果最好。

6) 将得到的谱包络通过Mel频率滤波器组处理, 将信号的频谱变换到Mel频域。滤波器的个数取24个。

7) 做离散余弦变换得到MTSMFCC参数。

3 实验结果与分析

本实验采用的数据库是通过16个人在不同SNR (0, 5, 10, 15, 无噪音) 下每人发音三遍得到的。将其中两遍的发音作为训练数据, 另外一遍发音作为训练数据。样本的词汇量是50词。对应的训练样本的数量为1 600, 测试样本的数量为800。语音信号采样率频率为11.025 k Hz, 语音每帧长为256点, 帧移128点。

本文是在Windows系统平台下, 利用Microsoft Visual Studio 2008进行编程来提取特征。利用libsvm工具箱进行识别。实验结果见表1。

从实验结果可以看出:

(1) 在不同信噪比下, 改进的新特征参数MTSMFCC比原始MFCC、SMFCC以及MTMFCC的识别率都要高。说明将两种改进结合同时提高了语音信号的噪声鲁棒性和时间鲁棒性。

(2) 对于同一种特征参数, 随着信噪比的增加, 识别率随着增加。

(3) 在无噪声的情况下, 改进的特征参数反而比原始的特征参数识别率低。这是因为原始的MF-CC在理想状态下识别率很高。而本文的改进算法主要针对有噪声的情况。

4 总结

本文针对Hamming窗的DFT虽然降低了频谱泄漏, 但却增加了频谱估计的方差值这一缺点提出的多窗频谱估计。并将得到的频谱通过滑动平均滤波器, 以消除基音频率的影响。在噪声环境下提高了识别率。但是此改进方法也存在缺点, 在理想环境下识别率却降低了。这就需要对这一算法做进一步的研究和改进, 力求在各种情况下识别率都有很大的提高。

摘要：语音的特征提取是说话人识别系统中的关键问题。在传统的Mel频率倒谱系数 (MFCC) 参数的基础上, 提出一种改进的MFCC特征提取算法。该算法着眼于语音的前端处理, 在预处理阶段, 利用SWCE窗函数, 对信号进行多窗频谱估计。并对得到的频谱进行平滑处理, 得到信号的谱包络。然后对信号的谱包络进行计算, 得到改进的MFCC参数。实验表明, 在不同噪声环境下, 与传统的MFCC算法相比, 改进的算法识别率提高四个百分点以上。

关键词：Mel频率倒谱系数,多窗频谱估计,滑动平均滤波,谱包络,说话人识别

参考文献

[1] 尹聪, 白静, 龚宬, 等.基于PLAR的说话人确认系统的噪音鲁棒性研究.清华大学学报 (自然科学版) , 2013;53 (6) :791 —795Yin Cong, Bai Jing, Gong Cheng, et al.Research on noise robust of speaker verification based on PLAR.Journal of Tsinghua University (Natural Science Edition) , 2013;53 (6) :791—795

[2] 张晶, 范明, 冯文全, 等.基于MFCC参数的说话人特征提取算法的改进.电声技术, 2009;33 (9) :61—64, 69Zhang Jing, Fan Ming, Feng Wenquan, et al.Improved speaker feature parameter extraction algorithm based on MFCC.Speech Technology, 2009;33 (09) :61—64, 69

[3] 张雪英.数字语音处理.北京:电子工业出版社, 2010:22—27Zhang Xueying.Digital speech processing.Beijing:Publishing House of Electronics Industry, 2010:22—27

[4] Sandberg J, Hansson-Sandsten M, Kinnunen T, et al.Multitaper estimation of frequency-warped cepstra with application to speaker verification.IEEE Signal Processing Letters, 2010;17 (4) :343—346

[5] Thomson D J.Spectrum estimation and harmonic analysis.Proceedings of the IEEE, 1982;70 (9) :1055—1096

[6] Riedel K S, Sidorenko A.Minimum bias multipletaper spectral estimation.IEEE Trans on Signal Processing, 1995;43 (1) :188—195

[7] Kinnunen T, Saeidi R, Sandberg J, et al.What else is new than the Hamming window?robust MFCCs for speaker recognition via multitapering.Processing Interspeech, 2010:2734—2737

[8] Hansson-Sandsten M, Sandberg J.Optimal cepstrum estimation using multiple windows.Processing ICASSP, 2009:3077—3080

本文来自 360文秘网(www.360wenmi.com)，转载请保留网址和出处

【系数估计】相关文章：

线性回归模型中回归系数的平衡LS估计06-24

空间估计05-10

函数估计06-04

估计预测06-26

参数估计软件05-24

交通状态估计05-29

密度核估计05-30

《大数的估计》教案05-25

样本量估计05-16