时间：2023-03-10 18:15:41
［关键词］零息债券 最小二乘法 三次样条
一、引言
近几年来，经济和金融的领域中，短期利率的模型被广泛地应用。这是因为短期利率模型能够很好的刻画瞬时短期利率的变动规律。正是这个原因，目前有很多的短期利率模型的被提出。同时这类的模型很好地被应用到一些利率的衍生品的定价，如：债券、利率互换、利率的远期合约等利率衍生品。
Vasicek（1977）建立Gauss扩散的短期利率模型。但这个模型有一个不如人意，瞬时的利率有可能是负值。Coxetal.（1985）（CIR）建立根均值的短期利率模型。相对Vasicek，在参数满足Feller条件下，CIR模型的瞬时利率是正值。但这两个模型都是考虑参数是常数。参数是常数可能不适合市场的数据。众所周知，市场的利率是有期限结构的性质。而常数的参数，很难刻画出这种性质。Ho和Lee(1986)提出部分参数是时间的函数，但是在Ho和Lee模型中，波动率仍然是一个常数，也就是在任何时刻波动率是常数。这种想象也有点违背市场的数据。因此，Hull和White（1990）延拓了Vasicek模型和CIR模型，
在Hull和White模型中，所有的参数都是时间的函数，并在他们的文章，利用这两个模型对于利率互换的产品定价分析。进一步，Rogers（1995）在理论上分析Hull和White模型。
而在这篇文章中，我们将根据市场的数据来估计Hull和White（1990）延拓的Vasicek模型。考虑这个模型是因为这个模型有很好的分析性质，特别是应用到定价其他的利率衍生品。
二、模型
首先，假设短期利率r在风险中性测度下满足如下的随机微分方程：
（1）
其中， 均为时间函数，为标准的Brown运动，设如果是常数时候，（1）是Vasicek模型。
基于（1），在时间 时刻，到期日 支付为1单位的零息债券价格有如下表达式，
（2）
其中，E（••）表示条件数学期望，表示为在t时刻的观察值。如果是一个常数时，基于Vasicek模型（2）有解析表达式，这些结果可以参考Hull（2009）以及姜礼尚等（2008）。考虑到参数是时间时候，Rogers(1995)也证明（2）有显示的表达式。
根据Feynman-Kac公式，在鞅的测度下，（1）适合下面的微偏分方程，
（3）
众所周知（3）有一个仿射结构解，
(4)
其中，。也就是在到期日零息债券价格为一个单位。即使，和有解析的表达式，但他们是多重积分计算问题（参考Rogers(1995)）。因此，我们把这类问题归结为求方程的问题。
基于（4），满足下面的常微分方程
（Model 1）
其中，其终端条件 (Model 1)能够被发现在Hull和White（1990）。
相应Model1的离散格式，设其中为了简化说明，设Model 1的离散格式为：
（5）
其中，（5）终端条件，
基于（4），债券价格为
（6）
其中，。因此，对于不同的到期日，（6）能算出其相应的价格。
我们的问题是在给定市场零息债券价格下，估计Model 1中的参数，首先，假设在t时刻观察到零息债券价格。理论上，当 T足够大时候，我们能估计未来任何时刻的参数值，事实上，市场上的零息价格期限是有限，因此我们将讨论基于离散的市场数据来估计依赖时间函数的参数。
其次，在我们的计算方法中，需要计算理论的价格和市场价格的残值，这里我们将使用常微分方程计算。即使债券价格有显示表达式，但它是一个三重积分计算问题，这就可能产生离散的误差。使用偏微分方程计算时有显著的误差（看图1）。
图1基于Vasicek模型通过常微方程和偏微分计算零息债券价格数值解与精确解的比较，这里假设参数是常数，基于Vasicek模型，当r比较小时，常微分方程的数值解比较好，这归因于计算其偏微分方程时，其在r大于零的值依赖于r小于零的值。另一原因在于偏微分数值结果依赖对空间变量r的离散误差。而常微分方程数值解仅依赖t，r是一个外在的变量。
三、参数估计计算方法
为了简化问题我们假设当前时刻，进一步我们假设有N个市场价格这就意味着到期日相同的债券有不同价格，这个根源于初始利率的不同，这个假设隐含着，市场上债券价格通过同一个模型而得到的。
进一步我们假设，市场价格满足如下方程
（7）
其中，被称为理论价格，表示理论价格和市场价格之间的误差，直接求解（7）可能是不适定，这个归因于理论价格可能是非线性的形式解而且市场价格报价存在误差。这就导致来自（7）的解可能是不适定的，因此，我们可以把问题转化为求下列极小值的问题
（8）
其中，表示所在领域，求（8）极小解，，可通过非线性最小二乘法，如果直接求解（8）时，对于至少有一个值是不能被计算，这个依赖于求解理论价格时离散格式。因此我们将通过三次样条差值来计算端点的问题。
四、实证结果
我们使用美国国家债券每天交易市场数据，这篇文章所有使用的数据都来源于www.ustreas.gov.由于其报价是通过收益率曲线来实现，因此需要通过现值表达式转化为债券的价格。设，当前时刻t，在T时刻其相应在收益率曲线上收益率为 ，那么其相应的零息债券价格为
（9）
如：t=0，y=0.015，T=5，其债券价格为0.92774.
我们考虑2009年整年的到期为5年债券每天数据 ，同时我们限制所有参数取值落在
1.另一类方法是正则化方法，这种方法不仅仅可以解决端点的问题，而且处理参数光滑性的问题，我们将在今后研究中将使用这类的方法估计参数。
2.选择美国国家债券数据由于这类债券几乎是无风险债券，而对于一般债券数据由于隐含违约风险因素，可能不满足我们所考虑的模型。
3.该收益率曲线刻画是实际收益率而不是名义收益率。
区间[0,1]。对于初始利率将适应美国国库券到期为一周的每天收益率为初始利率。
表1列出基于Model 1参数估计结果，这里我们假设离散的步长 在表1中最后一列的数据通过三次样条插值得到。通过观察
表1：这个表格描述Model 1的参数估计数值结果，最后一行是（8）的残值，所有的数据都乘以100。
数据，对于所有的数据变化大约或不超过0.1%，这就表明参数有足够的光滑性。
为了测试算法的稳定性问题，我们在收益率和初始利率上添加扰动项，即：
（10）
（11）
其中是待定常数，是均值为0方差为1的正太分布的随机变量。
在表2中，我们列出了RMES的值，通过观察表2中数据，对于初始利率相对于初始的收益率数据更加稳定，根均值误差几乎没有什么变化，虽然对于初始的收益率变动其根均值有微小的变化，但这样的变化整体来说是可以接受的，如我们预期的一样，我们的计算方法是稳定。
五、结论
我们主要考虑基于零息债券价格估计Model 1的参数，其中参数是依赖时间函数。我们所使用方法是众所周知的最小二乘法。端点的处理通过三次样条插值计算而得到。通过数值结果可知，基于债券市场数据，参数估计值是稳定的，同时参数的估计值的光滑性比较好。
参考文献：
[1] O. A. Vasicek, An Equilibrium Characterization of the Term Structure,Journal of Financial Economics, 5 (1977): 177-188.
[2] J. C. Cox, J. E. Ingersoll, and S. A. Ross, A Theory of the Term Structureof Interest Rate, Econometrica, 53 (1985): 385-407.
[3] T. S. Y. Ho and S. B. Lee, Term Structure Movements and Pricing InterestRate Contigent Claims, Journal of Finance, 41, (1986): 1011-1029.7
[4] J. Hull and A. White, Pricing Interest-Rate-Derivative Securities, The
Review of Financial Studies, Vol. 3(4) (1990), 573-592.
[5] L. G. G., Rogers, Which model for term-structure of interest rates shouldon use? Mathematical Finace, 65, 93-116, 1995
[6] J. C. Hull, Options, Futures and Other Derivatives (7th ed.), PrenticeEducation, Inc., 2009.
[7]姜礼尚，徐承龙，任学敏，李少华等，金融衍生产品定价的数学模型于案例分析，高等教育出版社，2008
关键词：规划模型；资产负债管理；随机情景生成；参数估计
中图分类号：F810文献标识码：A文章编号：1001-6260（2008）02-0093-06
一、引言
随着计算能力的显著提升和算法研究的巨大进步，随机规划正在成为一个强有力的工具，在金融机构和个人长期资产负债管理中发挥着越来越重要的作用，并已取得了巨大的经济效益。一个比较著名的案例是Mulvey等（2000）为Towers PerrinTillinghast公司开发的一套随机资产负债管理系统，自1991年以来已在欧洲、北美、亚洲等地区的19个国家（地区）里为数千家养老金公司及保险公司提供决策咨询服务。US WEST养老基金因此节约了4.5亿～10亿美元的机会成本。
于立勇（2004）认为，与其它金融资产负债管理（ALM）模型相比，随机规划模型的主要优点是：问题刻画方便，可以把来自资产负债中的多种风险源整合在一个框架中进行考量；具有长期视野、可适应不同水平的风险规避条件且能把交易费用、市场的不完备性、税收、交易费用和管理规则等因素纳入考虑范围，具有较大的灵活性。此外，它易于求出数值解，便于给出可操作性投资策略建议。
ALM的随机规划模型是Mulvey等人于1998年提出的，它一般都是围绕着一个称为情景生成器的随机预测系统及一个资产负债决策优化模型进行设计，模型各模块之间的关系如图1所示。随机预测系统用来产生大量具有代表性的情景元素以模拟未来的不确定性，每个情景描绘了一个多阶段规划期间内模型经济变量的演变路径。把这些情景输入资产负债决策优化模型并求出模型的全局最优解，以此给资产负债管理提供决策建议。
由于所有的决策建议都是基于生成的情景加以优化得到的，因此情景质量的高低自然决定决策建议的质量。所生成的情景在多大程度上体现了未来的不确定性是个非常重要的问题。目前，在学术界和实务界主要有以下几种情景生成的方法：历史数据重构法、Russell的向量自回归模型、ORTEC的带有均衡条件的向量自回归模型VaR法及随机微分方程法。前几种方法主要借助历史数据，研究各变量的时间序列特点生成情景。而随机微分方程法则是系统考量经济因素之间内在逻辑关系，在一个统一的框架下，用一系列的随机微分方程刻画各变量发展演变特征。相对而言，这种方法产生的情景自然比仅仅借助历史数据产生的情景更准确，更有代表性。
Towers Perrin公司所用的随机预测系统CAP：Link堪称这方面的典范，它由一套包含关键经济变量的随机微分方程构成，这些经济变量包括价格、工资、通货膨胀率、不同久期的利率、股票的红利收益和红利增长率等。用微分方程刻画各变量演变进程，并产生有代表性的情景元素，且这些情景包含的经济变量在全球多个国家内同步确定。为了保证情景生成的质量，这些随机微分方程的参数估计要尽可能地准确。
二、参数估计的方法
目前，随机微分方程常用的参数估计方法主要有：极大似然估计法（maximum likelihood，ML）、广义矩方法（generalized method of moments，GMM）和模拟矩估计法（simulated moment estimation，SME）。三种方法中，模拟矩估计法不需要参数向量与模型变量间具备明显的关系，而只需将它们用模拟值代替即可，故其能广泛地应用到资产定价模型参数估计中。
（一）极大似然估计法（ML）
极大似然估计法最初于1912年由英国统计学家R.A.费希尔提出，它利用样本分布密度构造似然函数来求出参数的最大似然估计。下面以一个短期利率的随机微分方程的参数估计为例，对极大似然估计法的应用进行简要的介绍。利率方程如下：
广义矩方法(generalized method of moments,GMM)是关于参数估计的又一种方法。GMM的一般表述是由Hansen (1982)提出的。GMM最大的优点是仅需要一些矩条件而不是整个密度。很多的估计量都可以视为GMM 的特例，如普通最小二乘估计量、工具变量法估计量、两阶段最小二乘估计量、非线性联立方程系统的估计量以及动态理性预期模型的估计量等，在很多情况下极大似然估计量也可看作是 GMM 的一个特例。许多计量经济学模型不是通过完全的分布假设而是通过矩条件来设定，例如带有不可观测的个体影响的动态平面数据模型和含有理性预期的微观经济模型，这些模型通常是使用GMM方法来估计的。
一般地，GMM估计方法就是极小化下式：
ML、GMM及SME都是利用极小化误差来估计参数值的，虽然广泛应用于经济模型的参数估计中，但他们本身还存在着一些不足之处。首先，它们对历史数据的依赖过于严重。因为经济形势往往瞬息万变，产生的情景一般要能代表未来的一种趋势。当形势变化较大时，ML估计便会产生较大误差。此外，ML估计需要有最大似然函数，当模型比较复杂时，这点往往难以满足。GMM模型和SME方法需要把模拟的矩代替模型自身矩，但如果矩本身是不够稳定的，估计必然会产生问题。
一般来讲，优良的参数估计模型应具备如下特点：满足误差最小原则、保证出自模型的样本概率最大，以及产生与真实描述统计一致的样本。针对这些特点，Mulvey等（1999）提出一个更为有效的估计方法综合参数估计（integrated parameter estimation，IPE）模型。它不仅能较好地满足上述三个特征，同时可以应用目标规划的权重来控制各种目标的相对重要性。
三、 综合参数估计法（IPE）
综合参数估计法（integrated parameter estimation,IPE）是在模拟矩估计法的基础上发展起来的。它在两个方面对SME方法进行了改善：首先，增大了目标函数集，它的目标函数集中不仅包括矩向量，还包括相关描述统计量，如自相关、分布百分位模型等。其次，目标函数的适应类型广，使IPE方法具备较大的灵活性。
测度、可行域没有限制且Ψ仅包括矩统计量时，IPE等价于SME，这说明SME法仅仅是IPE的一种特殊形式。
IPE的目标函数集不仅包括低阶矩向量，还包括高阶矩向量，如峰度和偏度等。此外，还包含相关的描述统计量，如均值、方差、自相关、分布百分位模型等。一般来讲，90th和10th的百分位统计量就已经足够，当然也可以把两个百分位之间的距离纳入统计目标。总之，IPE的目标函数可以包括任何性质，只要他们能表示为参数集的函数就可以。
IPE要求随机模型产生的样本满足给定的描述统计量，要求以这些描述统计量作为目标，并用式（11）对样本偏离区间进行限制。参数的可行域由使用者直接设限，也可以通过约束施加限制。
IPE的目标函数值是各个单独目标函数值的加权平均，每个目标函数根据其自身相对重要性赋予相应的惩罚权重。权重的选择需要慎重，要充分考虑投资者所处的经济环境。例如，风险中性的投资者重视短期资产价格，因而比风险规避者重视长期资产配置，他们的目标函数的权重必然不同。此外，通过历史数据预测将来情景在一般情况下是可以接受的，但是未来毕竟不是过去的重复，在某些时候，政策变化带来的经济趋势的变动是很剧烈的，这就需要对模型进行及时调整，以反映变化中的情况，IPE方法可以加以适当的调整来适应这种情况。
Hull（1993）提出一种类似的调整方法。首先估计参数，然后基于估计出的参数对资产集进行定价，对其用市场价格进行评估。如果偏差较大，就需要继续改进参数估计，直至满足一定的条件。这种方法能保证得到符合市场波动的一套参数。比如，在给定的利率期限结构下，这种方法用来定价是必要和足够的，但用来预测长期经济环境时，就略显不足了。
下面用例子具体介绍IPE的使用方法。对于方程（7），传统的统计方法由于其含有一个随机波动项，难以估计。用IPE方法来估计参数，首先设置目标函数集Ψ，取Ψ={均值、方差、自回归、90th、75th、25th，和10th百分位}对于i∈Ψ，让Si和Ti表示第i个模型及目标统计量。相应的IPE模型具有以下结构：
其中，^Wi是对应的目标规划的权重。采用二次误差函数，在给定的点，用根据方程（7）模拟出的Si值来算出函数值。问题的求解难度取决于决策变量的数目和类型以及目标函数的类型。
IPE参数估计问题是一个非凸规划问题，Mulley等（1999）提出了一个适应性记忆规划（adaptive memory programming，AMP）方法，取得了较好的效果。具体过程为：首先找出局部最优解，然后利用拓扑法扩大寻找范围，试图找到全局最优解。与其他全局求解算子相比，其有以下优点：（1）对目标函数要求不高，可以通过短期或长期记忆来加速搜索过程；（2）通过拓扑法可以很方便地从当前局部最优点向潜在更优点移动，便于找到潜在的最优解，且可以很方便地处理多目标函数。
四、对三种参数估计方法的比较
参数 MLGMMIPE α0.19170.27050.21 β-0.0235-0.0312-0.02 r0.73880.32150.29 σ0.12400.32720.35分别用ML、GMM和IPE方法对模型（1）的参数进行估计。所用数据为英国1980年1月―1995年3月债券月收益数据。参数估计结果见表1。
IPE的目标函数是一个包括均值、方差、1，2，3阶自相关及90th-10th百分位及75th-25th区间函数集合。函数值是二次惩罚项的线性组合，所有的统计量赋予相同权重。方差的权重由于不稳定降为0.5。目标统计量设为英国的债券票面利率的历史数据（1980.1―1995.12）。结果如表2所示。
五、小结
本文介绍了用于随机情景生成系统的参数估计方法。试验表明，IPE方法比ML法和GMM方法具备较大的优势，表现出更大的灵活性和更小的偏差，而且适应面也更广泛，可广泛用在一些比较复杂的模型参数估算中。此外， IPE权重及惩罚函数的选择也较为灵活。
在以往的多阶段随机规划中，随机模型参数的估计和决策模型的优化是各自独立的两个部分。但这些问题往往是紧密联系在一起的，因为估计和抽样带来的误差同样会导致次优的决策建议。因此，把IPE参数估计和优化决策结合起来是一个需要进一步研究的方向。
参考文献：
于立勇. 2004. 基于随机规划的动态投资组合选择［D］. 北京：中国科学院数学与系统科学研究院博士学位论文:6.
COX J, INGERSOL J, ROSS S. 1985. A theory of the term structure of interest rates ［J］. Econometrica, 53:363-384.
DUFFIE D, SINGLETON K. 1993. Simulated moments estimation of Markov models of asset prices ［J］. Econometrica, 61:929-952.
HANSEN L. 1982. Large sample properties of generalized method of moments estimators ［J］. Econometrica, 50:1029-1054.
HULL J C. 1993. Options, futures, and other derivative securities［M］. 2nd edition: Englewood Cliffs, NJ Prentice Hall.
MERTON R. 1973. An intertemporal capital asset pricing model ［J］. Econometrica, 41:867-887.
MULVEY J M, RUSH R, SWEENEY J. 1998. Generating scenarios for global financial planning systems ［J］. International Journal of Forecasting, 14:291-298
MULVEY J M, ROSENBAUM D P, SHRTTY B. 1999. Theory and methodology-parameter estimation in stochastic scenario generation systems ［J］. European Journal of Operational Research, 118:563-577
MULVEY J M, GOULD G, MORGAN C. 2000. An asset and liability management system for Towers Perrin -Tillinghast ［J］. ABI/INFORM Global, 30(1):96-97.
MURTAGH B A, SAUNDER M A. 1982. A projected lagrange algorithm and its implementation for sparse nonlinear constraints ［J］ . Mathematical Programming, 14: 41-72.
The Parameter Estimation of Stochastic Scenario Generation Model
WEI Faming LIANG Dan CHEN Weizhong
(Institute of Modern Finance, Tongji University, Shanghai 200092)
Abstract：Stochastic programming model (SPM) is widely used in asset and liability management by many financial institutions and individual investors for its special merits. It's an essential step to describe future uncertainty (often named scenario generation) accurately to ensure its successful utilization. Stochastic differential equation is an important way to generate scenarios. Model parameters need to be estimated accurately for getting more representative scenarios. This paper briefly introduces some usual ways for estimating parameters and offers a better method named integrated parameter estimation (IPE) in detail. Then it empirically compares the effects of ML, GMM and IPE and offers the direction for further study of parameter estimation.
Keywords: programming model; ALM; stochastic scenario generation; parameter estimation
为了更准确地估计散焦模糊点扩散函数的散焦半径，提出了一种基于阶跃边缘的参数估计方法。首先利用Canny算子对散焦模糊图像进行边缘检测，然后利用Hough变换方法提取边缘图像中存在的边缘直线，最后对沿直线法线方向的各像素，计算其灰度级的二阶导数，导数值中的最大值则对应模糊阶跃边缘的边界，进而计算出点扩散函数的参数。实验表明，该方法能够快速准确地计算出散焦模糊参数。
【关键词】散焦模糊 参数辨识 阶跃边缘 点扩散函数
散焦模糊是因调焦不准确而造成的图像模糊，并丢失了一些重要的高频成分，而这些高频成分恰恰蕴含着图像中最重要的信息，使得人们对图像的辨识能力下降。散焦模糊广泛存在于图像应用的各个领域，它造成的图像信息丢失，严重影响了应用效果，制约了这些领域的进一步发展。因此，对散焦模糊点扩散函数参数的估计方法进行研究有着重要的实用价值和意义。
对散焦模糊点扩散函数参数估计的研究是图像复原的一个重要研究领域，并已提出多种估计散焦模糊PSF参数的算法，这些方法大概可分为三类：（1）基于空域的参数估计。（2）基于变换域的参数估计。（3）基于迭代技术的参数估计。随着人工神经网络和遗传算法等新兴技术的出现，人们将其应用到图像处理领域，提出了基于这些知识的参数估计算法。基于以上的理论研究，本文提出了一种新的估计散焦模糊点扩散函数参数的方法。
2 散焦模糊理论分析
2.1 散焦模糊点扩散函数
通过对成像的原理和过程的分析，通常认为成像系统具有空间移不变行，因此一幅降质图像的降质过程在空间域可用如下过程来表示：
（2-1）
式中g（x，y）为降质图像，f（x，y）为原始清晰图像，h（x，y）为点扩散函数，即成像系统对点光源的响应，n（x，y）表示加性噪声，*表示二维卷积操作。在上述表达式中，通常假设噪声为高斯白噪声，尤其是在噪声不明显的情况下，可忽略。那么在上式中，点扩散函数就是惟一未知项。
在散焦模糊点扩散函数的几种模型中，由于圆盘模型只需估计出散焦半径便可计算出PSF，在计算上更容易、更简便，因此在参数估计时通常选用圆盘模型。其表达式如下所示：
（2-2）
其中，r表示散焦半径，决定了散焦模糊的程度，即是参数估计方法所需要估计的参数。
2.2 阶跃边缘
图像中的边缘对应着相邻的两个类型区域的分界线，表示一个区域的结束和另一个区域的开始。设s（y）为一条沿x轴的理想阶跃边缘，可用下式来表示：
（2-3）
系统对s（y）的响应称为边缘扩散函数。
在计算出散焦阶跃边缘区域的左右边界LI和Lr后，可根据下式计算出散焦半径r，即散焦模糊点扩散函数的参数：
（2-4）
3 散焦模糊点扩散函数参数估计算法
3.1 模糊阶跃边缘图像分析
经过散焦模糊后的阶跃边缘，其在图像中呈现为一个模糊区域，称之为模糊阶跃边缘。基于阶跃边缘的散焦模糊图像的点扩散函数参数估计，其关键是根据直线边缘确定模糊阶跃边缘的模糊区域的边界宽度，进而计算出散焦模糊点扩散函数的模糊半径。在散焦模糊图像中，阶跃边缘的模糊区域与检测到的直线边缘的关系如图1所示：
图1中，设为检测到的模糊阶跃边缘的任意一条直线，长度为；定义以为高、区间为宽的区域为阶跃边缘的支撑区域；定义以为高、区间为宽的区域为模糊阶跃边缘的支撑区域；和分别为直线到模糊阶跃边缘左侧和右侧边界的距离。
3.2 确定模糊阶跃边缘的边界
以计算图3-1中直线的右侧边界为例，介绍利用二阶导数确定边界的方法。设散焦模糊图像为f（x，y），其沿x轴方向的一阶偏导为，在离散情况下可以用差分来表示。为了处理上的方便，本文对计算出的一阶导数值执行取绝对值操作。根据阶跃边缘的散焦模糊图像的特点可知，在范围内，直线右侧沿其法线方向的灰度值的变化率即导数值会在模糊阶跃边缘区域的边界两侧出现较大的变化，因此可根据相邻导数值出现较大变化的点的位置确定其右边界。
为了找到相邻两个点的差值的最大值，可以对计算出的一阶导数值再对进行一次一阶偏导操作，即对执行对的二阶偏导操作，对计算出的导数值依然取绝对值。同理，可以用二阶导数值确定直线的左边界。
3.3 计算散焦模糊PSF的参数
在计算出散焦阶跃边缘的左右边界和后，便可根据（2-4）计算出散焦半径r，即散焦模糊点扩散函数的参数。
4 结论
本文提出了一种基于阶跃边缘的散焦参数估计方法，对于存在阶跃边缘的散焦模糊图像，能够快速准确地估计出散焦半径。该方法不但对阶跃边缘有良好的估计效果，对于轻度平滑的边缘也同样适用。该方法的优点是计算简便、时间复杂度低，比较适合于实时应用系统。
作者单位
【关键词】 数字调制信号 参数 估计算法
前言：随着科学技术的发展，通信技术被应用到各个领域，但在通信技术研究中却发现，经常出现数字信号并不能在通信道中传输，这就需要应用到调制技术。而数字调制技术因便于集成，同时具有加密能力，因此被广泛应用，对于数字调制信号参数估计算法的研究也是在此基础上发展起来的。
一、数字信号处理模块概述
数字信号处理模块的应用主要是为了完成数字变频、调制以及处理等工作，对于其电路而言，构成部分有ARM、FPGA以及必要的接口电路等，在数字信号处理模块中，其基础带有无线电技术的通用数字平台，该平台不仅具有良好的开放性，还具有一定的通用作用，能够很好的吸纳各种软件，以便满足不同用户与环境的需求[1]。而软件无线电的使用则需要以现代通信为基础，其中心则是数字信号处理技术，同时还需要得到微电子技术的支持。由于这一通用数字平台有ARM、FPGA的支持，可以很好控制AGC。总的来说，数字信号处理模块对数字调制信号参数估计具有重要作用。
二、数字调制信号参数估计算法的实现
要研究数字调制信号参数估计算法的实现，具体来讲可以从以下几方面入手：
2.1高阶统计量载波估计
之所以选用高阶统计量载波作为实现数字调制信号参数估计算法，主要是由于这种算法具有很好的抗噪声作用，即便是在信号相对微弱的情况下依然具有良好的识别能力，但值得注意的是，在该算法计算量相对大的情况下，应通过一些措施适当减少计算量，并做好粗估计，确定好估计范围，然后再搜索，然而，在实际搜索中，还要关注搜索步长与计算量之间的关系，当搜索步长变小的情况下，估计精度便会提升，反之估计精度则会降低，这就需要联系实际情况提出一种既能减少计算量，又能确保估计精度且可以变化的搜索方法。
为实现这一目标，应先估计接收机噪声，不让接收机出现热燥声等情况，进而保证了算法性能，可以让信号源在有线方式的作用下完成信号接收，在信号源传输BPSK信号时，就可以从接收机中看到一些高低起伏的信号，在这一过程中要忽略热燥声的影响[2]。在研究中发现，当SNR为5时，信号会受到周围噪声的干扰，产生高低起伏变化，但由于该信号为调制信号，它的频域就会比其他部分高一些，然而，要避免直接估计频域，防止出现误差，只要选择粗估计即可，粗估计只是转移了频谱，还有效防止了信号丢失的发生。
2.2小波变换码元速率估计
对于小波变换来说，能够精准的确定盲信号码元速率，同时，利用小波变换码元速率完成数字调制信号参数估计，可以有效减少不利因素对算法的影响，信号预处理也可以减少频偏的出现，这些都为小波变换码元速率的应用奠定了基础。在利用小波变换码元速率出来前端滤波时，要先完成信号处理，这样不仅可以有效提高信噪比，还能减少外界因素所带来的不良影响[3]。通过研究发现，在SNR为5时，并没有出现信号频偏，且BPSK信号还出现了部分峰值，在既定周期中小波变换系数还存在相对稳定的情况，但有些BPSK信号也会受到一些影响，继而影响到峰值，一般来讲，当信噪比多大或多小的情况下都会影响信号，很可能还会出现算法失效的情况，面对这一现实，就需要通过频偏减少小波变换系数所带来的影响，进而强化算法性能。
2.3循环谱联估计算法
对于循环谱来说，整个计算过程计算量相对较大，这就需要应用FAM算法完成计算，在利用该算法中，应先将数据分成多个小部分，然后为各个数据加窗，同时为各个数据进行傅里叶转换，再完成相关处理，且再次实现傅里叶转换，最后将有价值的数据展示出来。在实际应用中发现，部分信号为摆脱带外辐射的影响，在正式发射前会形成一定的基带信号脉冲，而这些脉冲也会带来不利影响，这就需要应用汉明窗实现加窗，这也是提高算法性能的有效方法。
结束语：通过以上研究得知，数字调制信号的应用是为强化数字信号传输效率，保证传输效果，而数字调制信号的实现也需要得到一些参数估计算法的支持，针对这种情况，本文联系数字信号处理模块基本情况，重点研究了三种参数估计算法，并指出了这些算法在实现中容易出现的问题，同时也提出了合理解决办法，希望能为相关人士带来有效参考，做好数字调制信号设计工作，提高通信能力。
参 考 文 献
[1]许华，王爱粉，杨晓宇. 常规数字通信信号信噪比估计综述[J]. 信号处理，2013，06：723-733.
[2]林财永，张国毅. 基于AR功率谱和STFT的混合调制信号参数估计算法[J]. 雷达与对抗，2013，03：31-34+45.
关键词： 粒子群优化算法； 非线性系统； 参数估计； 优化
中国分类号：TP301.6 文献标志码：A 文章编号：1006-8228（2012）04-34-02
An algorithm of parameter estimation of nonlinear system model
Wei Zhengfang, Qi Mingjun
（Hebi Occupation Technology College， Hebi， Henan 458030， China）
Abstract： Aiming at the diversity of nonlinear system model, it is proposed in this article a parameter estimation method based on particle group optimization algorithm that is applicable to a variety of nonlinear models. The result shows that the particle group optimization algorithm for parameter estimation of nonlinear system model is an effective tool.
Key words： particle group optimization algorithm； nonlinear system； parameter estimation； optimization
0 引言
非线性系统广泛地存在于人们的生产生活中，但是，目前我们对非线性系统的认识还不够深入，不能像线性系统那样，把所涉及的模型全部规范化，从而使辩识方法也规范化。非线性模型的表达方式相对比较复杂，目前还很少有人研究各种表达方式是否存在等效关系，因此，暂时还没有找到对所有非线性模型都适用的参数模型估计方法[1]。如果能找到一种不依赖于非线性模型的表达方式的参数估计方法，那么，也就找到了对一般非线性模型系统进行参数估计的方法[2]。
粒子群优化算法[3](Particle Swarm Optimaziton，简称PSO)是由Kennedy博士和Eberhart博士于1995年提出的一种基于群体智能的优化算法，它源于对鸟群群体运动行为的研究，即粒子群优化算法模拟鸟群的捕食行为。设想这样一个场景：一群鸟在随机搜索食物，在这个区域里只有一块食物，所有的鸟都不知道食物在那里，但是他们知道当前的位置离食物还有多远，那么找到食物的最优策略是什么呢？最简单有效的方法就是搜寻目前离食物最近的鸟的周围区域。粒子群优化算法从这种模型中得到启示并用于解决一些优化问题。粒子群优化算法中，每个优化问题的解都是搜索空间中的一只鸟，我们称之为“粒子”。所有的粒子都有一个由被优化的函数决定的适应值(fitness value)，每个粒子还有一个速度决定他们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。粒子群优化算法将粒子解初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个"极值"来更新自己，第一个就是粒子本身所找到的最优解，这个解叫做个体极值pBest，另一个极值是整个种群目前找到的最优解，这个极值是全局极值gBest。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。其基本思想[4]是模拟自然界生物的群体行为来构造解的随机优化算法，即从一组初始解群开始迭代，逐步淘汰较差的解，产生更好的解，直到满足某种收敛指标，即得到了问题的最优解。假设在一个n维的目标搜索空间中，有m个粒子组成一个群落，其中第i个粒子在n维搜索空间中的位置表示为一个n维向量，每个粒子的位置代表一个潜在的解。设为粒子i的当前位置；为粒子i当前飞行的速度；为粒子i所经历的最好位置，也就是粒子i所经历过的具有最好适应值的位置，称为个体最优位置；为整个粒子群直至当前时刻搜索到的最优位置，称为全局最优位置。将带入目标函数计算出其适应值，根据适应值的大小可以衡量的优劣。每个粒子的位置和速度按下文中式⑶和⑷两个公式迭代求得。用j 表示粒子的第j维(j=1，2，…，n)，i表示第i个粒子(i=1,2，…，m)，t表示第t代，c1、c2为加速度常数，通常在0～2间取值，c1调节粒子向自身最优位置飞行的步长，c2调节粒子向全局最优位置飞行的步长。，为两个相互独立的随机函数。为了减小在进化过程中粒子离开搜索空间的可能性，vij通常限定于一定范围内，即。如果问题的搜索空间限定在内，则可设定。迭代中若粒子的位置和速度超出了限定范围，则取边界值。代表第i个粒子在t时刻位置到直至t时刻搜索到的最优位置的距离，代表第i个粒子在t时刻位置到整个粒子群直至t时刻搜索到的最优位置的距离。公式⑵用于计算粒子的速度，如当前是t时刻，则粒子在t+1时刻速度是由当前时刻的速度、当前位置与该粒子的局部最优位置的距离、当前位置与全局最优位置的距离共同决定的；公式⑶用于计算粒子速度更新后的位置，它由粒子当前位置和粒子更新后的速度决定。所有粒子的初始位置和速度随机产生，然后根据上述两个公式进行迭代，不断变化它们的速度和位置，直到找到满意解或达到最大的迭代次数为止(粒子的位置即是要寻找的解)。因此，粒子群优化算法具有多点寻优、并行处理等特点。而且粒子群优化算法的搜索过程是从初始解群开始，以模型对应的适应函数作为寻优判据，从而直接对解群进行操作，而与模型的具体表达方式无关。这就决定了粒子群优化算法可适用于一般非线性系统模型的参数估计。
1 基于粒子群优化算法的非线性系统模型参数估计方法
1.1 问题的提出
一般非线性系统模型可用式⑴表示。
⑴
式中，y(t)为系统输出向量；u(t')为系统输入向量，0≤t'≤t；，θ为待定参数向量。f的形式已知，且u(t')已知。现已知y(t)的一组实际测量的离散数据y0(t)，t=1,2，…，n。要求根据已知的y0(t)的值估计出θ的值。
为了能够进行辩识，式⑴所代表的非线性系统模型还必须满足以下假设：①y必须可测；② 每个参数必须与输出y有关，即参数可估计；③系统的信噪比足够大，以至噪声可忽略不计；④ 只要参数确定，通过系统仿真可得到确定的输出值；⑤系统在有限时间t内不发散，即y值不趋于无穷大。
1.2 基于粒子群优化算法的参数估计方法
本文用一种改进粒子群优化算法自动寻找θ。具体步骤如下。
⑴确定适应函数：在已知各参数值的基础上，基于式⑴，可通过仿真实验求得各个时间的系统输出数值y(t)。辨识的目的是要使求得的系统输出数值y(t)尽量接近已知的系统输出数值，越接近说明仿真的效果越好，也就证明仿真所用的一组参数更接近实际参数值，因此应使这组参数对应的粒子群个体具有更小的适应值。所以，我们取y(t)曲线与y0(t)曲线之间距离的为适应值，
即： ⑵
⑵随机产生n个θ。
⑶计算适应值fi，再根据式⑵中确定的适应函数计算出各个θ对应的适应值fi。
⑷计算每个粒子的适应值。
⑸对于每个粒子，将其适应值与所经历过的最优位置的适应值进行比较，若较好，则将其作为当前的最优位置。
⑹对于每个粒子，将其适应值与全局所经历的最优位置的适应值进行比较，若较好，则将其作为当前的全局最优位置。
⑺根据下面2个公式对粒子的速度和位置进行更新；
⑶
⑷ ⑻如未达到结束条件（通常为足够好的适应值）或达到一个预设最大代数Gmax，则返回步骤2 直至算法收敛，即所有个体基本相同，适应值很难进一步提高为止。
2 仿真研究
为了体现粒子群算法能适用于多种非线性系统模型的优点，我们分别以非线性系统的传递函数模型[5]，非线性系统的状态空间模型及在非线性系统研究中应用较为广泛的Hammerstein 模型[6]为例进行仿真研究。
传递函数模型的形式如下：
可以看出，这是一个惯性环节加纯时滞模型，待估计的参数是比例系数K，惯性系数T 和时滞系数τ。在仿真实验中， 参数设置如下：学习因子c1=1.5，c2=2.5，惯性权重,T为最大代数,t为当前进化代数，在这里w将随着迭代次数的增加而逐渐减小，当w小于0.4时，将令w=0.4，即不再减小，以保证迭代后期粒子能够在一定空间探索更好的解。它们的群体规模是100，其他参数不变。在搜索过程中，以100代为上限(实际上，迭代50～80次即可得到满意结果)。仿真结果如表1所示。
表1 例1 参数估计结果
[[\&K\&T\&τ\&真实值\&10\&5\&9\&估计值\&10\&511\&9\&]]
在例1的仿真实验中，因为模型结构简单，待定参数较少，应用粒子群算法搜索较为容易，所以为了提高运算速度，参数精度定得较底，仅为小数点后一位，但从搜索结果来看，参数估计是令人满意的。实验说明了以下几点：①用粒子群优化算法进行参数估计是有效的；②在模型较简单，需要估计的参数较少时，用粒子群优化算法进行参数估计可达到比较满意的精度。
3 结束语
本文在利用粒子群优化算法对非线性系统模型参数估计方面作了一些尝试，得到了比较满意的结果。仿真实验结果表明，粒子群优化算法切实可行，对非线性系统模型参数估计具有一定的实际价值和理论意义。
参考文献：
[1] 徐南荣,宋文忠, 夏安邦. 系统辨识[M].1991.
[2] Goldberg D E Genetic Algorithms In Search ,Optimization [M] and Machine Learning[M] . Reading ,MA :Addison2Wesley,1989.
[3] Kennedy J, Eberhart R C.Particle swarm optimization[C].In: IEEE
International Conference on Neural Networks.Perth, Piscataway, NJ, Australia:IEEE Service Center, 1995; IV: 1942～1948
[4] 张鸿宾,郭建军, 遗传算法在曲线多边形近似中的应用[J].计算机学报,1999.10:1100~1104
[5] 方菲等.基于测试执行的失效数据建模研究[J].软件学报,1999.12:1233~1237
关键词：参数估计 适线法 模糊加权优化 相对隶属度
1 问题的提出
P—Ⅲ型分布做为我国水文分析计算中规定的线型，长期以来得到广泛应用，我国水文工作者对其参数估计问题作了大量研究工作，先后提出了适线法、权函数法、数值积分权函数法、Monte—Carlo法、概率权重矩法、下限法、动点动线适线法、密度函数法等等.金光炎对适线法作了大量工作［1］.适线法实际是一种优选法，即优选一组参数使之相应的理论曲线与经验点据最好吻合.但由于目估适线没有一个明确定量的拟合优度标准，使结果因人而异，往往达不到优选目的，找不到最佳结果.为了有一个客观定量标准，并利用计算机适线，近年来产生了一种新的方法——模型搜索适线法，其原理同目估适线法一样，寻找一组参数，使这组参数所决定的理论曲线与经验点据拟合最优.由于不同适线准则与不同的经验频率计算公式得到的结果有时差异很大，为了确定合理的适线准则和相应的经验频率计算公式，丛树铮等［2］对此问题作了大量的水文统计试验研究.
其实，在采用优化适线法时，尚有一个经验点与所配曲线之间离差的加权问题，由于各经验点的精度不同、在频率曲线上的位置不同、配线的目的不同，其在优化适线时的重要程度理应不同.现有的优化适线法均不能考虑这些因素，只能将各经验点与所配频率曲线之间的离差视为等权，因此计算结果受个别异常点影响十分敏感，致使所配曲线受少数异常点控制而时常偏离大多数经验点的分布趋势，这正是现行优化适线法计算结果欠稳定、平差能力低的症结所在.如何合理的确定各经验点与所配曲线之间离差在优化适线中的权重，使参数估计具有更好的稳健性与良好的平差能力，是一个值得探讨的问题.本文根据模糊数学的基本原理［3］，首次提出了经验点对所配曲线的相对隶属度概念及该相对隶属度的分析计算方法，进而提出了以该隶属度为权重的P—Ⅲ型分布参数模糊加权优化适线法，以期克服现行优化适线法之不足.
2 模糊加权优化适线法
若随机变量X服从P—Ⅲ型分布，则其概率度函数和分布函数分别为：
(1)(2)
其中α，β，γ为3个未知参数，它们与常用的3个特征参数Ex,Cv,Cs有如下关系
(3)
设想自总体X中抽取K组容量为M的样本系列，则可用适线法选配出K条频率曲线，若有某条频率曲线恰好完全通过Exm(即X(pm)=Exm,m=1,2,…，M)，则称此条频率曲线为随机变量X的理想最优频率曲线.因为只有通过它外延得出的设计值才是无偏的，其它频率曲线都是有偏的.
在理想最优频率曲线概念的基础上，我们进一步给出经验点据(xm,pm)对理想最优频率曲线的相对隶属度定义及其计算方法.
若以随机独立方式从随机变量X的总体中抽取一组容量为M的样本，将其按大小为序排列，则序位为m的变量xm(相应的经验频率为pm)称为次序统计量.它也是一个随机变量，若xm恰好落在理想最优频率曲线上，则xm对该曲线的相对隶属度为1；若xm距离理想最优频率曲线愈远则相对隶属度愈小，若xm距离曲线无穷远，则相对隶属度为0.
次序统计量xm是原随机变量X衍生出来的，其概率分布特性可以通过条件分布密度函数f(xm｜pm)来描述，按陈守煜教授提出的概率分析与模糊集分析相结合的模糊水文学基本思想［4］，本文假定条件分布密度函数f(xm｜pm)近似服从均值为Exm的正态分布，即f(xm｜pm)=N(Exm,σm)，并按下式定义xm对理想最优频率曲线的相对隶属度.
定义：若随机变量X的理想最优频率曲线为f(x,Ex,Cv,Cs)，则某经验点据(xm,pm)对f(x,Ex,Cv,Cs)的相对隶属度为
(4)
f(x,Ex,Cv,Cs)为理想最优频率曲线密度函数，其它符号同前.
在f(xm｜pm)服从正态分布的条件下，可按文献［2，5］所述的水文统计试验方法推求不同适线准则下式(4)中的σm，金光炎对此问题作了深入细致的研究，并绘制了绝对值准则和最小二乘法适线的B值诺模图［1，5］.
本文采用了金光炎的以上成果，并利用下式计算σm
式中S为X系列的标准差，n为样本容量，B为B值模图中的取值，它是Cs和pm的函数.
因为样本均值对期望值Ex的估计具有无偏性且方差最小，因而不需通过适线法来确定，于是便可得到本文提出的模糊加权优化适线法的数学模型。
(5)
式中c为优化准则参数.
3 各种参数估计方法的对比检验
3.1 计算结果 任何参数估计方法，即使在理论上正确，也必须通过大量适当的数据充分检验合格以后，才有应用的价值，本文采用刘光文教授［6］提出的与水文安全因数有关的绝对精度标准，即x1%和x0.1%(作为代表)的真误差不应超过8%—10%，用本文提出的模糊加权优化适线法和现行的几种参数估计方法，对大量的理想样本资料同时进行了对比计算，以下仅将与文献［6］中20组理想样本对
比的计算结果列于表1.
表1 假想样本对各种P—Ⅲ型参数估计方法的对比检验 样本序号* 样本容量(n) 项 目 假想理论值
(TV) 数值积分*
双权函数法
(NIDWF) 模糊加权最小
平方准则
(FWLSC) 模糊加权绝对值
最小准则
(FWLAVC) 1 19 cv
cs
x1%
x0.1%
1
1/2
1
2.51128
3.26556 0.99344(-0.7%)
0.48154(-3.7%)
0.96948(-3.1%)
2.42987(-3.2%)
3.14018(-3.8% 0.97942(2.1%)
0.5046(+1.3%)
0.8573(14.3%)
2.4184(-1.2%)
3.18867(-2.4% 0.97942(2.1%)
0.5089(+1.8%)
0.7474(25.3%)
2.44141(2.4%)
3.11913(4.5%) 2 29 cv
cs
x1%
x0.1% 1
1/2
1
2.51128
3.26556 0.99561(-0.4%)
0.48796(-2.4%)
0.98480(-1.5%)
2.45920(-2.1%)
3.18634(-2.4%) 0.98454(1.5%)
0.5097(+1.9%)
0.9041(-9.6%)
2.51049(0.0%)
3.24205(1.0%) 0.98454(1.5%)
0.4953(-0.9%)
0.6418(35.8%)
2.38095(5.2%)
2.99158(8.4%) 3 19 cv
cs
x1%
x0.1% 1
1/2
3/2
2.66518
3.61676 0.99072(-0.9%)
0.47775(-4.5%)
1.47686(-1.5%)
2.56057(-3.9%)
3.45267(-4.5%) 0.97031(3.0%)
0.5119(+2.4%)
1.2878(14.1%)
2.63823(1.0%)
3.52551(2.5%) 0.97031(-3.0%)
0.5084(+1.7%)
1.1152(-25.7%)
2.57268(-3.5%)
3.38809(-6.3%) 4 29
cv
cs
x1%
x0.1% 1
1/2
3/2
2.66518
3.61676 0.909381(-0.6%)
0.48511(-3.0%)
1.48752(-0.8%)
2.59586(-2.6%)
3.50864(-3.0%) 0.97775(-2.2%)
0.5095(+1.9%)
1.3440(-10.4%)
2.65240(-0.5%)
3.55990(-1.6%) 0.97775(-2.2%)
0.5115(+2.3%)
1.2061(-19.6%)
2.61539(-1.9%)
3.47005(-4.1%) 5 19
cv
cs
x1%
x0.1% 1
3/4
3/2
3.49776
4.92515 0.98608(-1.4%)
0.72573(-3.2%)
1.47852(-1.4%)
3.56030(-3.9%)
4.71006(-4.4%) 0.95547(-4.5%)
0.7839(+4.5%)
1.1547(-23.0%)
3.44339(-1.6%)
4.74307(-3.7%) 0.95547(-4.5%)
0.7689(+2.5%)
1.0853(-29.4%)
3.3642(-3.8%)
4.57124(-7.2%) 6 29
cv
cs
x1%
x0.1% 1
3/4
3/2
3.49776
4.92515 0.99072(-0.9%)
0.73903(-1.5%)
1.48737(-0.8%)
3.42364(-2.1%)
4.80974(-2.3%) 0.96663(-3.3%)
0.7724(+3.0%)
1.2649(-15.7%)
3.46916(-0.8%)
4.70729(-2.6%) 0.96663(-3.3%)
0.7580(+1.1%)
1.2038(-19.7%)
3.39449(-3.0%)
4.66096(-5.4%) 7 19
cv
cs
x1%
x0.1% 1
1/2
2
2.83259
3.95388 0.98850(-1.2%)
0.47137(-5.7%)
2.00927(+0.5%)
2.67059(-4.7%)
3.74683(-5.2%) 0.96261(-3.7%)
0.5066(+1.3%)
1.8071(-9.6%)
2.77295(-1.6%)
3.86085(-2.4%) 0.96262(-3.7%)
0.5058(+1.2%)
1.5219(-23.9%)
2.68934(-4.0%)
3.66066(-7.4%) 8 29
cv
cs
x1%
x0.1% 1
1/2
2
2.80259
3.95388 0.99233(-0.8%)
0.48052(-3.9%)
2.00328(+0.2%)
2.71220(-3.2%)
3.81141(-3.6%) 0.97203(-2.8%)
0.5049(+1.0%)
1.8811(-5.9%)
2.79284(-0.3%)
3.90880(-1.1%) 0.97203(-2.8%)
0.5073(+1.5%)
1.6324(-18.4%)
2.73194(-2.5%)
3.75257(-5.1%) 9 19
cv
cs
x1%
x0.1% 1
2/3
2
3.40345
4.93850 0.98466(-1.5%)
0.64245(-3.6%)
2.00993(+0.5%)
3.26847(-4.0%)
4.72963(-4.2%) 0.94014(-5.0%)
0.6951(+4.2%)
1.5256(-2.37%)
3.32663(-2.3%)
4.83837(-2.0%) 0.95014(-5.0%)
0.6954(+4.3%)
1.5155(-24.2%)
3.31947(-2.5%)
4.65110(-5.8%) 续表1
样本序号* 样本容量
(n) 项 目 假想理论值
(TV) 数值积分*
双权函数法
(NIDWF) 模糊加权最小
平方准则
(FWLSC) 模糊加权绝对值
最小准则
(FWLAVC) 10 29 cv
cs
x1%
x0.1% 1
2/3
2
3.40345
4.93850 0.98978(-1.0%)
0.65060(-2.4%)
2.00250(+0.1%)
3.31215(-2.7%)
4.79620(-2.9%) 0.96271(-3.7%)
0.6940(+4.0%)
1.7077(-14.6%)
3.39963(-0.1%)
4.83837(-2.0%) 0.96271(-3.7%)
0.6930(+3.9%)
1.6239(-18.8%)
3.36387(-1.2%)
4.75381(-3.7%) 11 19 cv
cs
x1%
x0.1% 1
1
2
4.60517
6.90776 0.97699(-2.3%)
0.99364(-0.6%)
2.01400(+0.7%)
4.48373(-2.6%)
6.72992(-2.6%) 0.92521(-7.5%)
1.0652(+6.5%)
1.9021(-4.9%)
4.77998(+3.8%)
7.14573(+3.4%) 0.92521(-7.5%)
0.8599(-14.0%)
1.4378(-28.1%)
3.82762(-16.9%)
5.41929(-21.5%) 12 29
cv
cs
x1%
x0.1% 1
1
2
4.60517
6.90076 0.98466(-1.5%)
0.99481(-0.5%)
2.00226(+0.1%)
4.51724(-1.9%)
6.77452(-1.9%) 0.94406(-5.6%)
1.0472(+4.7%)
1.9021(+4.9%)
4.75414(+3.2%)
7.11889(+3.1%) 0.94406(-5.6%)
0.9684(-3.2%)
1.7022(-14.9%)
4.34908(-5.6%)
6.35447(-8.0%) 13 19
cv
cs
x1%
x0.1% 1
5/2
5/8
3.40337
5.09259 0.98316(-1.7%)
0.59299(-5.1%)
2.60188(+4.1%)
3.25111(-4.5%)
4.87426(-4.8%) 0.94549(-5.4%)
0.6308(+0.9%)
2.2220(-11.1%)
3.34111(-0.2%)
4.90459(-3.7%) 0.94549(-5.4%)
0.5990(-4.2%)
1.9462(-22.2%)
3.14045(-7.7%)
4.49276(-11.8%) 14 29
cv
cs
x1%
x0.1% 1
5/1/2
5/8
3.40337
5.09259 0.98873(-1.1%)
0.60223(-3.6%)
2.52765(+1.1%)
3.28902(-3.4%)
4.90818(-3.6%) 0.95920(-4.1%)
0.6354(-1.7%)
2.3285(-6.9%)
3.40037(0.0%)
5.03408(-1.1%) 0.95920(-4.1%)
0.6209(-0.7%)
2.0570(-17.7%)
3.26343(-4.1%)
4.72475(-7.2%) 15 19
cv
cs
x1%
x0.1% 1
3/4
3
4.03853
6.36426 0.97690(-2.3%)
0.72046(-3.9%)
3.14988(+5.0%)
3.86713(-4.2%)
6.13317(-3.6%) 0.92712(-7.3%)
0.7508(+0.1%)
2.5727(-14.2%)
3.90632(-3.3%)
5.97705(-6.1%) 0.92712(-7.3%)
0.6686(-10.9%)
2.4218(-19.3%)
3.54320(-12.3%)
5.30654(-16.6%) 16 29
cv
cs
x1%
x0.1% 1
3/4
3
4.03853
6.36426 0.98448(-1.6%)
0.73826(-2.9%)
3.07152(+2.4%)
3.90830(-3.3%)
6.17219(-2.0%) 0.94530(-5.5%)
0.7713(+2.8%)
2.7660(-7.8%)
4.06518(+0.7%)
6.31981(-0.7%) 0.94530(-5.5%)
0.7119(-5.1%)
2.4904(-17.0%)
3.74470(-7.3%)
5.66888(-10.9%) 17 19
cv
cs
x1%
x0.1% 1
1
3
5.05138
8.15235 0.96920(-3.1%)
0.99675(-0.3%)
3.16031(+5.3%)
4.97404(-1.5%)
8.16870(+0.2%) 0.90283(-9.7%)
1.1035(+10.4%)
2.2418(-25.3%)
5.09906(+0.9%)
7.84849(-3.7%) 0.90283(-9.7%)
1.0692(+6.9%)
2.3048(-23.2%)
5.00685(-0.9%)
7.72542(-5.2%) 18 29
cv
cs
x1%
x0.1% 1
1
3
5.05138
8.15235 0.97931(-2.1%)
0.99508(-0.5%)
3.08385(+2.8%)
4.95773(-1.9%)
8.04434(-1.3%) 0.92707(-7.3%)
1.0778(+7.8%)
2.4421(-18.6%)
5.13616(+1.6%)
8.01065(-1.7%) 0.92707(-7.3%)
1.0472(+4.7%)
2.4619(-17.9%)
5.02789(-0.5%)
7.83777(-3.9%) 19 19
cv
cs
x1%
x0.1% 1
3/2
3
7.07706
11.72853 0.95380(-4.6%)
1.55446(+3.6%)
3.18574(+6.2%)
7.06141(-0.2%)
11.87708(+1.3%) 0.85424(-14.6%)
1.6422(+9.5%)
2.7897(-7.0%)
7.49812(+5.9%)
12.30353(+4.9%) 0.85424(-14.6%)
1.3063(-12.9%)
2.43470(-12.9%)
5.96459(-15.7%)
9.41022(-19.8%) 20 29
cv
cs
x1%
x0.1% 1
3/2
3
7.07706
11.72853 0.96896(-3.1%)
1.52306(+1.5%)
3.08189(+2.7%)
6.99291(-1.2%)
11.66507(-0.5%) 0.89060(-10.9%)
1.5864(+5.7%)
2.8152(-6.2%)
7.35641(+3.9%)
12.07083(+2.9%) 0.89060(-10.9%)
1.3903(-7.3%)
2.6342(-12.2%)
6.46617(-8.6%)
10.3966(-11.4%) * 计算结果引自文献［6］
又以淮河蚌埠站平均流量系列为例，进行了计算，其结果见表2，图1.
表2 各种方法计算结果对照
估计方法 Ex Cv Cs 矩 法 827.5 0.65 1.20 极大似然法 827.5 0.72 1.88 适 线 法 827.5 0.75 1.88 概率权矩法* 827.5 0.68 1.69 模糊加权适线法 827.5 0.71 1.52 *计算结果引自文献［8］
3.2 分析与讨论 本文仅对模糊加权优化适线法进行讨论，关于其他方法的讨论请参考文献［6］.
表1的计算结果表明：模糊加权最小二乘准则，对20组理想样本资料全都能达到精度标准.x1%的真误差绝对值平均为1.6%；而且正负误差基本持平，其中误差小于5%的有19组，最大误差为-6.1%；x0.1%的真误差绝对值平均为2.4%，其中正误差1.0%.负误差1.4%，误差小于5%的有19组，最大一组误差+5.9%，没有随Cv,Cs变化而出现明显的系统偏大或偏小现象，对不同参数的资料适应性强，计算结果稳定.
图1 淮河蚌埠站年平均流量频率曲线
模糊加权绝对值和最小准则，设计值估值全是负误差.x1%的真误差平均为5.2%，误差小于5%的有12组，5%—10%的5组，大于10%的3组，最大误差16.9%；x0.1%的误差均值为8.7%，小于5%的有4组，5%—10%的10组，大于10%的6组，最大误差为21.5%，不能达到绝对精度标准，计算结果系统偏小，对工程设计偏于不安全.计算结果还表明：模糊加权最小二乘准则适线法所估计的参数Cv有系统偏大，Cs有系统偏小的现象，这主要是由于用矩法估计有系统偏小的原因所致.但最终的设计估计值x1%,x0.1%没有出现明显的系统偏差，计算精度很高，这正是由于本方法有在，Cv,Cs估计值偏小的情况下，使Cv适当偏大，在一定程度上补偿了,Cs估值偏小的影响，使最终设计估值有很高的稳定性和计算精度，及较强的资料适应能力的优点.而模糊加权绝对值和最小准则适线法在估计偏小的情况下，不能很好地协调Cv,Cs估值，补偿各参数的估计值偏差，使最终设计估值系统偏小明显，难以达到绝对精度要求.
淮河蚌埠站实测资料的计算也表明，模糊加权最小二乘准则适线法所获得的理论频率曲线与实测资料拟合良好.
4 结 语
P—Ⅲ型曲线作为我国水文计算规范所规定的频率分析线型，在我国水文频率计算中得到广泛应用.但对于P—Ⅲ型分布参数的估计问题长期以来没有得到很好的解决，不同学派众说纷云，难达一致.本文提出了理想最优频率曲线和经验点据对理想最优频率曲线的相对隶属度概念及该相对隶属度的分析计算方法.然后从充分利用样本信息，增强参数估计方法的稳健性和平差能力等方面考虑，提出了以经验点对理想最优频率曲线相对隶属度为权重的模糊加权优化适线法.
模糊加权优化适线法可以充分利用样本信息(包括经验点的大小，排位和可能存在的误差等信息)，区分不同样本点在配线分析中的作用，较好的反映了“参数估计应以精度精高的实测点据为主，不要过于崇信可能存在较大误差的历史洪水，但又要适当考虑历史洪水点据所带来的信息”等配线分析过程中的模糊性思想，并把这些思想模型化、规范化，克服了一般优化适线方法在配线过程中等同看待精度不同的各经验点据，最终导致计算结果受个别异常点影响甚大，计算结果不稳定的缺点.
利用与水文安全因素有关的绝对精度检验标准，以理想样本和实测样本对不同模糊加权优化准则(模糊加平方和最小、模糊加权绝对值和最小)进行计算的结果表明：模糊加权绝对值和最小优化适线法，Cv,Cs,x1%,x0.1%均系统偏小，最大误差为-16.9%—-21.5%，难以达到容许误差标准；模糊加权平方和最小优化适线法，虽,Cs系统偏小，Cv系统偏大，但x1%,x0.1%计算精度很高，计算结果稳定，分析参数有系统偏差的原因，主要因用矩法估计系统偏小所致，Cv偏大是为补偿因，Cs估值偏小，使最终设计估计值稳定.
模糊加权平方和最小优化适线法计算结果稳定，对不同参数资料适应性强且参数,Cv,Cs间有很好的补偿效果，设计估值精度能达到容许误差标准，能很好拟合实测资料，是一种较好的P—Ⅲ型分布参数估计方法.
参 考 文 献 1 金光炎. 论水文频率计算中的适线法.水文，1992，(4).
2 从树铮，等. 水文频率计算中参数估计的统计试验研究.水利学报，1980，(3).
3 陈守煜. 系统模糊决策理论与应用.大连：大连理工大学出版社，1994.
4 陈守煜. 模糊水文学与水资源系统模糊优化原理.大连理工大学出版社，1990.
5 金光炎，等. Γ分布保证率修正值参数B的确定.水文，1991，(6).
6 刘光文. 皮尔逊Ⅲ型分布参数估计.水文，1990，(4，5).
7 刘治中. 数值积分权函数法推求P—Ⅲ型分布参数.水文，1987，(7).
8 宋德敦，等. 概率权重矩法及其在P—Ⅲ型分布中的应用.水利学报，1988，(3).
9 邱 林，等. 水文水资源系统模糊集分析理论及应用.郑州：黄河水利出版社，1995.
Weighted optimum curve－fitting method for estimating the
parametersof Pearson type－Ⅲ distribution
Abstract The weighted optimum curve－fitting method for estimating the parameters of Pearson type－Ⅲ distribution is presented in this paper. This method can increase the stability of calculation and reduce the influnce of errors in hydrologic data on estimation of the parameters as well as utilizing information in real hydrologic data.The method is tested by safety standard of hydraulic engineering with a great deal of ideal and real hydrologic data. The results show that the method of least squares criterion is stable and its precision satisfies the safety standard for hydraulic engineering.
关键词：瞬时测频；多项式拟合法；脉内参数估计
雷达信号的脉内特征是雷达信号相位特征的重要体现，是电子侦察中对雷达信号分选和识别的重要参数。因此，要可靠的分选和识别雷达信号就必须对雷达信号进行脉内特征分析。
在各种脉内分析方法中，分析信号的瞬时频率变化规律是进行脉内调制分析的有效手段。由于信号的频率是信号相位的变化率，所以对信号相位的估计自然可以用来对信号瞬时频率进行估计，从而进行相关脉内参数的估计，因此该方法的前提是能够得到准确的瞬时频率。
文中以线性调频信号以及相位编码信号为例，从侦察接收机的信号相位入手，利用瞬时相位法实现对瞬时频率的快速、准确的推算，然后利用多项式拟合法[1]来估计信号的相位参数。在Matlab中进行了仿真验证，同时给出了在某电子侦察系统中工程实现的过程。
4 结束语
本文讨论了利用信号相位建模进行瞬时频率估计算法，将频率的估计转化为求解信号模型参数的问题。在满足白噪声的情况下，由于在参数的估计过程中使用了最小二乘方法，因此能够有效的消除噪声的干扰，获得对信号瞬时频率的准确估计。基于瞬时频率进行脉内信号的调制特征的识别，还能够实现对常规脉冲、线性调频信号、非线性调频信号以及相位编码信号等调制类型信号的识别，提取相应信号的调制参数。该方法在一定信噪比条件下有较高的正确识别率，对脉宽变化不敏感，更加符合雷达侦察数字接收机上高速实现，并且在某雷达侦察数字接收机系统中已经取得成功验证，但是该类算法的问题在于，必须精确的建立相位模型，否则估计性能会受到很大的影响。
参考文献
[1] 林琳. 多项式相位信号的参数估计.实验科学与技术. 2006年6月第三期.
[2] 张春杰，郜丽鹏，司锡才. 瞬时相位法线性调频信号瞬时频率提取技术研究，弹箭与制导学报 TN911 2006.
[3] 齐国清，贾欣乐. 插值FFT估计正弦信号频率的精度分析. 电子学报 Vol.32 No.4 Apr. 2004.
[4] 胡来招. 数字瞬时测频—相位推算法测频[ J] . 电子对抗，2005.
关键词：离散匹配傅里叶变换;线性调频;参数估计;低信噪比
中图分类号：TN971 文献标识码：A
文章编号：1004373X(2009)0100103
LFM Signal Parameter Estimation Based on Discret Match Fourier Transform
AN Weigang1,CHENG Shaoyun2
(1.College of Electronic&Information,Jiangsu University of Science & Technology,Zhenjiang,212003,China;2.The 723 Institute of CSIC,Yangzhou,225003,China)
Abstract：LFM signal is a common signal form of low probability of intercept radar.How to detect the LFM signal of low SNR is one of the focus for people to investigate.The arithmetic is ameliorated on the basis of Discrete Match Fourier Transform(DFMT),then moni-component LFM signal and multi-component LFM signal is simulated by using the improving arithmetic.The simulation results show that DMFT can truly estimate parameters of the LFM signal in the low SNR,and the problems of across-terms do not exist.DMFT is a valid way to estimate the LFM signal parameter.
Keywords：discrete match Fourier transform;linear frequency modulation;parameter estimation;low SNR
线性调频（LFM）信号是低截获概率雷达［1,2］常用的一种信号形式,对LFM信号检测和参数估计一直备受人们的关注。针对该信号的处理方法有短时Fourier变换、Wigner-Ville变换［3］、分数阶Fourier［4］、Hough-Wigner［5,6］等,都存在分辨率不够高,交叉项严重或者运算量太大的问题。而匹配傅里叶变换是一种线性变换,不存在多分量信号交叉项的影响,能在低信噪比条件下检测信号,而且分辨率很高,是一种针对线性调频信号有效地进行参数估计的方法。
1 LFM信号形式
LFM信号的复数形式表示为：
s(t)=A(t)ej2π(f0t+1/2k0t2)
(1)
式中,A(t)为信号包络函数,f0为中心频率,k0=B/T为调频斜率,B为调频带宽,T为信号持续时间。
对于实际需要处理的信号,都是经过采样的离散信号。LFM信号的离散形式为：
s(n)=A(nTs)ej2π\[f0nTs+1/2k0(nTs)2\]
(2)
式中,Ts为采样时间间隔,如果信号持续时间为T,那么采样点数N=T/Ts。
2 DMFT基本原理
LFM信号s(t)的匹配傅里叶变换有如下两种形式［7］:
S(f,k)=∫T0s(t)e－j2π(ft+1/2kt2)(f+kt)dt
(3)
S(f,k)=∫T0s(t)e－j2π(ft+1/2kt2)tdt
(4)
称式(3)和式(4)分别为二阶匹配傅里叶变换和二步匹配傅里叶变换,对应其离散形式为:
S(f,k)=∑N-1[]n=0s(n)e-j2π\[fnT2+1/2k(nTs)2\](f+knTs)=
∑N-1[]n=0e-j2π\[(f0-f)nTs+1/2(k0-k)(nTs)2\](f+knTs)
(5)
S(f,k)=∑N-1[]n=0s(n)e-j2π\[fnTs+1/2k(nTs)2\]nTs=
∑N-1[]n=0e-j2π\[(f0-f)nTs+1/2(k0-k)(nTs)2\]nTs
(6)
由式(5)计算得到的谱图可称为离散二阶匹配傅里叶变换谱,其中k不为零,它表示了不同基条件下的匹配傅里叶变换；由式(6)计算得到的谱图可称为离散二步匹配傅里叶变换谱,它表示在不同频率补偿条件下信号的匹配傅里叶变换。
无论对离散二阶匹配傅里叶变换谱还是离散二步匹配傅里叶变换谱,在对应于信号(f0,k0)的位置上,信号能量会发生聚集,在谱上表现为尖峰。在匹配傅里叶变换谱分布图上进行二维搜索,尖峰的坐标(f0,k0)即为该LFM信号的线性频率f0和线性调频斜率k0。
由于离散二阶匹配傅里叶变换和离散二步匹配傅里叶变换具有不同的分辨率,通过文献［8~10］表明二步匹配傅里叶变换总是有比二阶匹配傅里叶变换更高的分辨率,因此下面的分析都采用离散二步匹配傅里叶变换进行LFM信号的检测和参数估计。
3 算法改进
对离散匹配傅里叶变换的二维搜索求极大值可以在低信噪比条件下获得较高精度的信号参数。但是当信号带宽增加,采样频率提高时,采样点数增加,运算量增大。下面从减少运算量的角度进行算法改进。对离散之后的信号进行离散匹配傅里叶变换,借助傅里叶变换的快速算法思想,实现离散匹配傅里叶变换的快速算法。对于长度为N的线性调频信号序列x(n),其N点离散匹配傅里叶变换定义如下：
Xc(f,k)=∑N-1[]n=0x(n)Wfn+kn2N,
0≤f,k≤N-1
(7)
其实质是将一个输入一维时间序列x(n)变换为关于线性频率和调频斜率的二维序列Xc(f,k),其中f为线性调频信号的初始频率,k为调频斜率。从式(7)可以看出,对于每一个固定的调频斜率k来说,{Xc(f,k)}0≤f,k≤N－1是信号x(n)Wkn2N的DFT；当调频斜率k=0时式(7)就转变为DFT。对式(7)进行改进得：
Xc(f,k)=∑N-1[]n=0x(n)Wfn+kn2N=∑N-1n=0x(n)Wkn2NWfnN=
FFT(x(n)Wkn2N),0≤f,k≤N-1
(8)
{Xc(f,k)}0≤f,k≤N－1的计算可以通过x(n)Wkn2N的快速傅里叶变换得到。在式(7)中需要N3次复数运算,经过式(8)变换,运算量减小为N2/2log2 N,提高了运算速度。为了提高该算法的估计精度,还可以在搜索范围内多次估计,分为粗估计和精估计。即首先在搜索范围内选择大步长,估计出信号参数,然后再在估计值邻近的区域内改变搜索步长重新估
计,从而达到需要的精度要求。
4 仿真实验
4.1 单分量LFM信号仿真
先对单分量LFM信号s(t)进行参数估计,s(t)=exp［j2π(f0t+1/2k0t2)］,经过下变频的信号线性频率f0=200 MHz,信号时宽T=5 μs,以带宽200 MHz的信号进行仿真,比较在不同信噪比条件下信号参数估计的结果,如图1所示。
图1 信号带宽B=200 MHz,SNR=-15 dB参数估计
表1为B=200 MHz时不同信噪比情况下初始频率和调频斜率的测量值与其对应真值(f0=200 MHz,k0=4.0×1012 Hz/s)的绝对误差。
表1 测量值与真值的绝对误差情况表
绝对误差
SNR /dB
0-5-10-15-16-18
Δf /MHz0.3000.3000.3000.3001.17017.920
Δk /1012 Hz/s0.0100.0100.0100.0100.1900.850
从以上结果可以看出,该方法对信号参数的估计有较高的精度,在SNR=-15 dB的情况下还能估计信号参数,这是一般的时频分析方法不能比拟的。SNR低于-15 dB时,参数估计绝对误差将逐步增大,信号经过离散匹配傅里叶变换淹没在随机噪声中,无法正确检测信号。
4.2 多分量LFM信号仿真
离散匹配傅里叶变换是一种线性变换,所以在对多分量信号进行分析时不会产生交叉项。但是信号中强分量LFM信号的旁瓣可能大于弱信号的主瓣峰值,影响到多分量LFM信号的分辨和参数估计。为了解决这个问题,借助“Clean”的思想［11,12］：首先计算多分量LFM信号的离散二步匹配傅里叶变换,然后进行二维搜索找极大值。并根据峰值的位置和大小估计最强LFM信号分量的幅度、初始频率和调制斜率,然后由上述参数重构LFM信号并从信号之减去,最后将处理过的信号重复上述过程估计下一个LFM信号的参数。
多信号的参数估计仿真采用如下信号：
s(t)=a1ej2π(f1t+k1t2)+a2ej2π(f2t+k2t2)
其中,a1=3,a2=2,f1=200 MHz,f2=220 MHz,B1=200 MHz,B2=220 MHz,T=4 μs,SNR=-5 dB。
进行第一次DMFT之后信号频谱如图2所示,只出现强信号分量的一个峰值,弱信号的峰值淹没在强信号分量的旁瓣中。此时,在图2中搜索谱峰最大值,得出强信号的分量：f1=2.002×108,k1=4.96×1013,a1=s(f1,k1)/Ts∑N－1n=0n=2.98,构造LFM信号第一个分量s1(t),得到剩余信号s2(t),再进行一次二步DMFT,对其余LFM信号分量估计,得到如图3所示结果,估计得到第二个分量的参数：f2=2.197×108,k2=5.53×1013,a2=1.89。
图2 信号s(t)的MFT三维谱
图3 去除s1(t)后弱信号s2(t)的MFT三维谱
上面的仿真结果表明,离散匹配傅里叶变换结合“clean”思想是一种检测多分量LFM信号的有效的方法。仿真进一步表明,当较小分量的信噪比不小于-15 dB时,LFM信号的参数估计能
达到较高的精度。随着信噪比的进一步降低,参数估计精度将下降,无法正确估计信号的参数。
5 结 语
首先介绍了LFM信号的形式以及DMFT的基本原理,然后从减小运算量的角度对DMFT算法进行改进,最后分别对单分量和多分量LFM信号进行Matlab仿真,结果表明,DMFT能够在低SNR情况下估计出LFM信号的参数,不存在多分量信号交叉项问题,而且运用本文改进的算法运算量较小,在对低截获概率雷达信号的处理中将有广阔的应用前景。
参考文献
［1］Guosui L,Hong G Weimin.The Analysis and Design of Modern Low Probability of Intercept Radar［J］.IEEE.2001:120-124.
［2］史林,彭燕,张毓峰.一种低截获概率雷达信号及其信号处理［J］.现代雷达,2003,25(6):26-28.
［3］孙泓波,顾红.基于互Wigner-Ville分布的SAR运动目标检测［J］.电子学报,2002,30(3):347-350.
［4］齐林,陶然,周思永.基于分数阶Fourier变换的多分量LFM信号的检测和参数估计［J］.中国科学(E辑),2003,33(8):749-759.
［5］ Liu Jiancheng,Zhou Jie,Wang Guoyu.Detection Perfor-mance of Linear Frequency Modulated Signals Based on Wigner-Hough Transform［J］.Signal Processing,2005,21(4A):546-549.
［6］孙晓昶,皇甫堪.基于Wigner-Hough变换的多分量LFM信号检测及离散计算方法［J］.电子学报,2003,31(2):241-244.
［7］刘爱芳,朱晓华,刘中.基于离散匹配傅里叶变换的多分量LFM信号检测和参数估计［J］.信号处理,2002,18(6):539-542.
［8］王盛利,李士国,倪晋麟,等.一种新的变换――离散匹配傅里叶变换［J］.电子学报,2001,29(3):403-405.
［9］王盛利,倪晋麟,张光义.匹配傅里叶变换的分辨力［J］.系统工程与电子技术,2002,24(4):29-32.
［10］Wang Shengli,Zhu Li.Changeable Sampling Processing for Chirp Signals in the SAR ［A］.Fourth International Conference on Signal Processing［C］.Beijing,China,1998:1 473-1 476.
［11］张贤达,保铮.非平稳信号分析与处理［M］.北京:国防工业出版社,1998.
［12］王宏禹.非平稳信号处理\[M\].北京:国防工业出版社,1999.
作者简介 安伟刚 男,1983年出生,硕士研究生。研究方向为电子对抗。
关键词: 图像复原；运动模糊；点扩散函数；倒谱；Radon变换
中图分类号:TP 391
文献标志码:A文章编号:1672-8513(2012)02-0150-04
Parameter Estimation of Motion Blur Based on Cepstrum Analysis
GUO Hongwei1, ZHU Jiaxing1, CUN Ning2
(1.Engineering College, Honghe University, Mengzi 661100, China; 2.Kunming Experiment Station, Yunnan Provincial Administration of Radio and Television, Kunming 650031, China)
Abstract: Motion blur is one of the most common blurs that degrades images，estimating the parameters of the PSF is necessary for image restoration. For uniform linear motion blurred images, the motion direction and blur extent are two key parameters determining the point spread function. The paper analyzed the characteristics of the spectrum and cepstrum of motion blurred images, and an algorithm which estimated the point spread function in the cepstrum domain was proposed. The coordinate of two negative peak points in cepstrum was used to estimate the blur extent and to obtain the absolute value for cepstrum, and then the Radon transform was employed to find motion direction. The experiments for the parameter estimation of motion blurred images showed the validity of the proposed method.
Key words: image restoration; motion blur; point spread function; cepstrum; Radon transform
在图像拍摄记录的过程中，由于被摄物与成像系统产生相对运动造成图像降质而导致的图像模糊称为运动模糊.经典图像复原技术都是以图像退化的某种先验知识比如点扩散函数(Point Spread Function，PSF) 已知为基础的.运动模糊图像的点扩散函数由模糊方向和模糊尺度共同确定.
匀速直线运动所造成的图像模糊更具有一般性和代表性，对它的成因进行分析，出现了多种辨识模糊参数的方法，主要的有空域法、频域法和倒谱法.空域法直接在图像空域利用微分、相关等方法计算估计退化参数，文献［1］提出了一种基于方向微分的运动模糊方向鉴别方法，文献［2］利用对退化图像进行差分自相关的方法检测图像模糊尺度.频域法利用匀速运动模糊图像频谱具有规则明暗条纹的特性分析估计退化参数［3-7］，如文献［3］ 的Radon变换法，文献［4］的相关系数法等.倒谱法对图像频谱取对数，然后再进行傅里叶反变换，可以分离退化图像模糊信息和原始信息，进而获取图像模糊参数［8-11］.文献［8］ 通过对退化图像的倒谱图实施灰度变换,运用Canny边缘检测提取出倒谱，实施Radon变换确定模糊方向，用人工方式在边缘检测后的倒谱中判读亮条纹两端亮点坐标来计算模糊尺度.由于运动模糊图像倒谱中只存在1条亮条纹，其它像素为低灰度值（黑色），文献［8］中的灰度变换和边缘检测增加了算法的计算量，而未从本质上提高参数估算精度，其模糊尺度的估算需人工测量，增大了判别误差，且不便于实现算法自动化.文献［9-10］利用倒谱具有左右对称的性质，取倒谱的右半平面自动检测灰度极小值（为负数）点的坐标，采取几何运算方式估算模糊方向和模糊尺度.这种方法计算量较小，但由于图像像素的离散性，仅利用一个点来估算模糊方向，除某些特殊角度（0°,45°,90°）外，容易出现较大估算误差，模糊尺度越小，对模糊方向的估算误差会越大.文献［11］对倒谱图进行Radon变换鉴别出运动模糊方向后，通过倒谱旋转将模糊方向旋转到水平轴方向，再进行模糊尺度的鉴别.本文通过对模糊图像频谱和倒谱的分析，结合以上文献算法的经验，提出一种基于倒谱分析的运动模糊参数估计方法.利用倒谱中的灰度负峰值点坐标，采取几何运算方式得到模糊尺度；对倒谱取绝对值后用Radon变换鉴别出模糊方向.实验数据表明，本文算法估计的模糊方向和模糊尺度误差较小.
1 运动模糊图像的分析
图像的退化过程如图1所示.
对于匀速直线运动模糊图像来说,点扩展函数可以描述为：
h(x,y)=1L,x2+y2≤L且yx=tan θ;0, 其它 . (2)
式中，L为模糊尺度，θ为运动方向与x轴正向夹角.若θ=0°，即水平匀速直线运动，则点扩展函数变为：
h(x,y)=1L,0≤x≤L;
0,其它 . (3)
1.1 运动模糊图像的频谱分析
对式(3)做Fourier变换可得：
H(u,v)=sin(πuL)πuLe－jπuL .(4)
sin xx称作辛格函数，用sin c(x)表示，是偶函数.在x为0时函数取得最大值1，x为nπ（n是不为0的整数）时函数值为0.所以，匀速直线运动模糊图像的频谱G(u,v)有一系列的平行暗条纹，这些暗条纹与水平方向垂直,且位置与sin(πuL)πuL函数的零点对应.图2所示为运动模糊图像和其频谱图，可看出频谱中暗条纹方向与运动方向垂直，且条纹间距与模糊尺度L成反比关系.一些文献［3-7］就是基于这一频谱特性估算运动模糊参数的.
1.2 运动模糊图像的倒谱分析
倒谱是对数倒频谱的简称，图像g(x,y)的倒谱定义［12］如下：
Cg(p,q)=F－1［lg｜G(u,v)｜］ .(5)
式中，G(u,v)是图像g(x,y)的傅里叶变换，F－1［•］表示傅里叶逆变换.为使G(u,v)=0时，函数有意义，在实际工程应用中，一幅图像的倒谱通常用下式计算:
Cg(p,q)=F－1{lg［1+|G(u,v)|］} .(6)
无噪声影响时，图像退化的倒谱描述为：
Cg(p,q)=Ch(p,q)+Cf(p,q) .(7)
式中各项为式(1)中对应的倒谱.可见,空域卷积在倒谱域变成了加法,因此可以比较容易地分离出模糊信息.为了便于说明运动模糊图像倒谱的特性，选取大小为65×65像素，中心点灰度为256，其余像素灰度为0的仿真图像进行模糊处理（模糊角度45°、模糊尺度20 像素）.图3是其倒谱的平面图和三维图，图3(a)显示倒谱中沿运动模糊的方向上有1条亮带，因此检测出亮带与水平方向的夹角即为运动模糊角度.图3(b)的倒谱三维图由2部分组成，在运动模糊方向上，一部分是正峰值成分，反映未退化图像特性；另一部分是负峰值成分，反映模糊系统的特性.2部分在图上占有的区域不同,2个负峰值点之间的距离即为运动模糊尺度的2倍.
2 Radon变换
二维空间中函数f(x,y)的Radon变换定义为：
R(β,ρ)=Sf(x,y)δ(ρ－xcosβ－ysinβ)dxdy .(8)
其中，S是被积函数所在区域，β是旋转角度，ρ是原点到一直线的距离，δ是冲激函数.式(8)表示f(x,y)在一条直线上的投影，即沿该直线的线积分.这条直线斜率为tan(β+π2)，截距为ρsinβ.图4所示给出了Radon变换的示意.对于一幅图像，当β一定，ρ取遍所有值时，就得到其在β方向上的投影；再改变β值，就可以得到图像在不同方向上的投影.
Radon变换与Hough变换一样，都可以将图像平面上的线转换成参数平面上的点.Radon变换可以理解为图像在ρβ平面上的投影，ρβ平面上的每一点都对应着xy平面上的1条线.因此，图像中1条高灰度值的线就会在ρβ平面上形成1个亮点；低灰度值的线则会形成1个暗点.图像平面中直线的检测就转化为在参数平面中对亮、暗点的检测.对图像做1°~180°的Radon变换，结果可用矩阵R表示.该矩阵有180 列，每一列向量是图像在1个角度上沿1族直线积分的投影值.对模糊图像倒谱做1°~180°的Radon变换，由其特征可知,在β=(θ+90°)时, 由于积分直线束与倒谱中的亮带平行，所得投影向量中的极大值即为矩阵R中的最大值.
图5所示为对图3(a)的倒谱做1°~180°的Radon变换，取每个角度上Radon变换的极大值形成的曲线.最大值对应角度β=135°,所以用Radon变换检测到的运动模糊角度θ=β－90°=45°,这一检测结果与实际模糊角度吻合.
3 算法实现
根据以上的分析，提出估计运动模糊参数（模糊尺度L和模糊方向θ）的算法过程如下：
1) 由式(6)计算退化图像倒谱Cg(p,q),并且移位使p=0,q=0位于中心位置;
2) 检测出倒谱Cg(p,q)中2个负峰值点的坐标Cg(p1,q1)和Cg(p2,q2),由下式计算得到模糊尺度：
L=12(p1－p2)2+(q1－q2)2 .(9)
由于倒谱具有对称性，检测2个负峰值点坐标时，可以先设置Cg(p,q)右半平面灰度全为0，检测此时倒谱的极小值，即得左半平面的负峰值点Cg(p1,q1)；然后再设置Cg(p,q)左半平面灰度全为0，检测出右半平面的负峰值点Cg(p2,q2).
3) 对倒谱的绝对值Cg(p,q)做1°~180°的Radon 变换，找出变换矩阵R中最大值对应的β角，则模糊角度：
θ=β－90° . (10)
4 实验与结果分析
在Matlab平台下对算法进行验证.表1为在不添加噪声情况下，对图6所示的Bridge图像做不同尺度和不同运动方向的模糊后，采用文献［10］及本文方法估计的点扩散函数.
由于运动模糊图像倒谱具有对称性，文献［10］采用倒谱右半平面负峰值点坐标和倒谱中心正峰值点坐标计算模糊尺度，与本文采用2个负峰值点坐标计算结果一致，模糊尺度在5～30像素时的估算误差均不超过0.6像素，估计精度较高.相比文献［10］用1个点计算运动方向，本文的Radon变换法由于利用了倒谱中整条亮带上的点，所以估算精度更高，几乎没有误差.而文献［10］估计的运动方向除某些特殊角度（0°,45°,90°）外，出现较大估计误差，模糊尺度越小，对模糊方向的估计误差越大.
图7(a)是一张贴在白色墙面上的A4纸打印广告，用相机拍摄时由于相机抖动造成运动模糊.使用本文方法检测出退化图像模糊方向为7°，模糊尺度为15.6像素，使用维纳滤波进行复原，并对复原图像进行了剪切去除边缘的振铃得图7(b).复原图像中，广告标题“红河学院送水部”字样已能辨识，说明参数的估计结果是正确的.
5 结语
分析运动模糊图像傅里叶频谱和倒谱的特征，根据倒谱中有2个负峰值点和亮条纹方向与运动方向一致的特点，提出用负峰值点坐标估算模糊尺度，用Radon变换检测模糊方向.通过对仿真运动模糊图像估计PSF参数的实验证明，该方法可以很准确地识别运动模糊尺度和方向.由于倒谱信息很微弱，在图像存在噪声时，倒谱中的模糊系统特性极易被掩盖，针对含噪声运动模糊图像的参数估计，算法还需做进一步改进.由于绝大多数实拍运动模糊图像不是理想的均速直线运动，对其参数估计和复原都比仿真图像困难，本文算法对大多数实拍图像的参数估计结果不理想.下一阶段将对实拍运动模糊图像的参数估计和复原算法开展研究.
参考文献:
［1］陈前荣,陆启生,成礼智.基于方向微分的运动模糊方向鉴别［J］.中国图像图形学报,2005,10(5):590-595.
［2］NOROUZI N,MOGHADDAM M E.Motion blur identification using image derivative［C］//IEEE International Symposium on Signal Processing and Information Technology. Yugoslavia: Sarajevo,2008:380-384.
［3］邓泽峰,熊有伦.基于频域方法的运动模糊方向识别［J］.光电工程,2007,34(10):98-101.
［4］TANAKA M,YONEJI K,OKUTOMI M.Motion blur parameter identification from a linearly blurred image［C］//International Conference on Consumer Electronics.US: Las Vegas, 2007:1-2.
［5］ARASHLOO S R,AHMADYFARD A.Fine estimation of blur parmeters for image restoration［C］//15th Intenational Conference on Digital Signal Processing.UK: Cardiff,2007:427-430.
［6］乐翔,程建,李民.一种改进的基于Radon变换的运动模糊图像参数估计方法［J］.红外与激光工程,2011,40(5):963-969.
［7］李宇成,贾宝华,杨光明.运动模糊图像的参数估计与恢复［J］.计算机工程与设计,2010,31(19):4247-4249.
［8］郭永彩,郭瑞瑞,高潮.运动模糊图像点扩展函数的参数鉴别［J］.仪器仪表学报,2010,31(5):1052-1057.
［9］邹文洁.一种新的基于倒谱法和Richardson-Lucy算法的运动图像盲复原方法［J］.南华大学学报:自然科学版,2010,24(2):61-67.
［10］谢伟,秦前清.基于倒频谱的运动模糊图像PSF参数估计［J］.武汉大学学报:信息科学版,2008,33(2):128-131.
［11］石明珠,许廷发,张坤.运动成像混合模糊的全变分图像复原［J］.光学精密工程,2011,19(8):1973-1981.
【关键词】 SHIBOR 利率期限结构 GMM
一、前言
SHBIOR是上海银行同业拆借利率的简称。随着以SHBIOR为基准的浮息债券、远期利率协议和利率互换的推出，SHBIOR作为我国货币市场基准利率的地位得到了进一步的加强。
我国货币市场上共有三种基准利率，分别是一年定期存款利率、银行间七天回购利率和上海同业拆借利率。
基于一年定期存款的产品，其优点是与人民币贷款的相关性比较高，而且市场上也有不少基于一年定存的浮息债，适合企业进行债务管理以及银行资产负债管理；其缺点是利率受到管制不随市场变化而变化。特别是在目前紧缩货币政策背景下，市场单边需求现象严重，交投不活跃。另外，从利率市场化的目标来看，存贷款利率会逐渐放开，一年定存产品也会随之逐渐淡出市场。
七天回购利率是市场化程度比较高的货币市场利率，被国内很多机构作为交易的成本，再加上人民币远期外汇市场进行投机或套利的需求，所以基于FR007的互换产品是目前交易最活跃的产品。但是，由于七天回购利率是短期资金利率，只反映了货币市场的资金供求情况，与金融机构和企业资产负债的相关性并不高，起不到套期保值的作用。
SHIBOR利率的推出，对构建我国货币市场基准利率具有重要的意义。目前拆借利率、银行票据利率和部分浮息金融债已经与SHIBOR利率进行挂钩，SHIBOR利率已经具有“基准利率”的雏形。SHIBOR的推出意味着我国金融市场建立起了真正的价格中枢，将为人民币汇率产品如远期汇率、汇率期权，人民币利率衍生产品如利率互换、利率期权及利率期货等提供定价基准，极大地推动了人民币衍生市场的发展。
二、理论基础
1、模型介绍及模型转变
利率期限结构主要讨论金融资产到期时的收益与到期期限这两者之间的关系，本文要研究的就是利率期限结构里服从均值回复的三个单因子利率模型。
Vasicek模型的利率期限结构：dr=a（b-r）dt+?滓dw。
Brennan-Schwartz模型的利率期限结构：dr=a（b-r）dt+?滓rdw。
其中，a、b、?滓为常数，三个模型均考虑了均值回复，a其中表示短期利率向均值回复的速度，b为均值水平，?滓反映利率的波动情况。
由于三个模型均为连续的形式，难以直接用GMM进行估计，故根据欧拉方程将上述连续时间模型化成离散形式。
上述就是本文要研究的模型。通过估计?琢、?茁来估计连续型中的a、b，离散型中的?滓即为连续型中的?滓。本文采用的方法是现在使用比较广泛的广义矩估计方法。
2、GMM方法介绍
三、实证研究
1、数据选择
本文选取自2006年10月8日SHIBOR试运行以来到2008年5月9日共398个O/N SHIBOR利率数据（数据来自www.省略），对数据进行连续复利处理，得数据的基本情况如下图1。
2、实证与结果
编写SAS程序，得上述3个模型的GMM参数估计如下。
（注：括号内为P值。）
从上述SAS检验的数据来看，在置信水平为94%的情况下，三个模型的?琢、?茁、?滓参数估计均相当显著。总体看来，三个模型均能反映SHIBOR的波动情况，其中又以Brennan-Schwartz模型最具解释力，且波动幅度最小 。因此，Brennan-Schwartz模型是最符合短期SHIBOR市场利率的。
下面把研究结果与基于R007的利率模型参数估计进行比较，文献[4]中GMM参数估计结果如下表2。
通过上述比较可知，SHIBOR的市场利率比R007的波动要大些，主要是因为SHIBOR推出的时间较短。另外，文献[4]的数据选取自1997年6月20日到2006年4月10日，而本文数据的选取是从2006年10月8日到2008年5月9日，且包含了央行连续多次加息的过程，所以没有文献[4]的稳定。
四、政策与建议
虽然SHIBOR利率没有R007稳定，基于SHIBOR的衍生品交易量也没有基于R007的衍生品的交易量大，但是我国的SHIBOE利率体系已初步显示出一些基准利率的特性，在未来成为我国货币基准利率也已成为大势所趋。所以，对SHIBOR体系进行研究具有重要意义。但是要真正实现这种跨越，则需在保证其稳定运行的情况下，提高它的市场效力。
1、深化同业拆借市场，完善SHIBOR利率形成机制
一方面，在完善公司治理结构、内控制度及信息披露制度和加强外部监管的基础上，逐步放宽对非银行金融机构的交易限额和期限管理，以促进资金供给和需求主体多元化，把SHIBOR利率培育成高度市场化的货币市场基准利率。另一方面，提高SHIBOR利率的成交性。
2、提高SHIBOR利率与金融市场利率的相关性
目前，转贴现利率、部分浮息金融债已经与SHIBOR利率进行挂钩，但从市场所占比例和市场影响力来看，这种相关性仍非常小。下一步首先应多鼓励政策性银行和商业银行多发行以SHIBOR为基准利率的浮动利率债券。在此基础上，推动短期融资券、公司债以及企业债利率与SHIBOR利率进行挂钩。
3、加强宣传，引导金融机构正确认识SHIBOR的地位
进一步推进SHIBOR与区域性金融产品定价的联系，加大市场基准利率SHIBOR的宣传力度，逐步加深社会公众对SHIBOR的认识，引导金融机构产品定价逐步以SHIBOR为基准，以应对未来利率市场化挑战，防范产品定价风险。
【参考文献】
[1] 谢赤、吴雄伟：基于Vasicek 和CIR 模型中的中国货币市场利率行为实证分析[J].中国管理科学，2002（6）.
[2] 傅曼丽、董荣杰、屠梅曾：动态利率模型估计方法的一个实证检验[J].华中科技大学学报（自然科学版），2005（4）.
[3] 李子奈、潘文卿：计量经济学[M].北京：高等教育出版社，2006.
[4] 赵峰：利率期限结构理论、实证与应用[D].上海：上海社会科学院，2006.