核密度计

核密度计（精选7篇）

核密度计 篇1

引言

随着国民经济的发展、科技的进步, 对生产过程中自动化程度的要求越来越高, 核密度计能够实时地对生产过程中的密度变化进行检测, 能够随时掌握测量料液的密度、浓度, 提高工艺水平, 节约能源, 保证生产的安全性, 提高机组的安全性, 减少运行人员工作量, 通过自动控制回路, 实现工艺过程的自动控制, 有效的排除故障。如今核密度计已经在国外的炼油, 化工等技术中广泛的应用, 近几年来, 国内的一些石油化工厂也开始采用核密度计这种装置。而且国内的核密度计的研制也取得了突破的进展。文章将介绍一种基于γ射线原理研制的数字智能型核密度测量装置, 它的最大优点是采用非接触测量、数字显示、稳定可靠, 不破坏、连续自动地测量密度。且在安装、调试、维护时, 不影响正常生产, 特别适合恶劣环境下密度的测量。为厂家节省大量的人力、物力, 具有很好的社会效益和生态效益。核密度计对密度的精确可靠监测, 保持密度在正常范围内运行, 合理控制密度的波动幅度, 并达到闭环控制, 对保证设备安全经济运行有着十分重要的意义。

1γ射线辐射技术的原理

γ射线是一种电磁波、光子流。与无线电波, 红外光, 紫外线和可见光一样。当放射性物质的原子核能级跃迁时就会产生γ射线。其能量高, 波长短, 穿透力极强。

γ射线与物质发生相互作用后射线的强度会减弱, 放射性同位素γ射线源衰变过程中稳定地释放出一定强度的射线, γ射线与物质发生相互作用, 在其通过被测物质时其强度将进行衰减, 当被测物质的密度和厚度越大, 强度衰减越大。当放射源固定后, 射线在传播时的衰减就只取决于介质的密度和传播的距离。当一束准直的单能γ射线, 沿水平方向垂直通过被测物质时, 其透射强度和被测物质的关系为[1]:

其中:

I0为未穿过被测物质前的射线强度;I为穿过被测物质后的射线强度;μ为质量吸收系数, 单位为cm2/g;d为厚度, 单位为cm;ρ为密度, 单位为g/cm3。

通过推导可以得出:

由 (2) 式可见, 只要能准确测量出被介质吸收后的射线强度, 即可随时计算出被测物质的密度ρ。

2核密度计的工作原理

2.1核密度计检测系统的组成

采用铯-137辐射源为发射源, 该射源具有半衰期长 (30年) , 在满足系统需求的同时便于屏蔽防护;利用核物理测量原理进行设计, 以点状辐射源发射、闪烁探测器进行接收, 采用射线准直方式测量, 并以微处理器为核心的数字智能核密度测量装置进行运算, 在监测过程中加入温度、放射源衰减补偿, 采用LED显示, 脉冲在线监测及RS485接口。核密度的测量系统如图1所示。有放射源铯-137, 闪烁探测器, 微处理器、和控制回路[2]。

2.2核密度计的工作原理

铯-137辐射源放射的射线穿过被测物, 被闪烁探测器接收, 并将接收到的γ射线转换成电脉冲信号, 对探测器得到的信号通过放大电路放大进行预处理, 并滤除掉噪声和干扰信号。将输出信号整形为矩形脉冲信号, 以适应微处理器的测量要求。接收到的射线强度, 决定闪烁探测器的输出脉冲的计数频率, 用计数频率来测量强度时, 如果测量时间为T, 脉冲数为n, 则平均计数率N为[3]:

式 (3) 无法满足工业参数连续快速测量的需要。因此, 采用定时计数方法来优化原有算法, 将测量时间分割成c个小定时单元Δt, 在第i个定时单元Δti, 其计数脉冲为Δni, 则式 (3) 可表示为:

铯-137辐射源与探测器之间的距离以及它们之间的工艺介质都是固定不变的, 被测物所吸收射线的能量就取决于介质的密度[4]。介质所吸收的射线越多, 闪烁探测器接收到的就越少, 介质的密度就越小, 反之, 介质的密度越大。通过将介质的密度与探测器的读数一一对应, 就可以方便的得到介质密度的读数[5]。

按照工艺介质的流向, 核密度计的安装位置最佳是在环管反应器的垂直段入口处, 距离轴流泵较近。介质自下而上的流动可以防止在管壁上挂垢, 使得较重的颗粒能够处于悬浮的状态。

3结束语

该核密度测量装置相比传统的测量仪表, 具有如下的优势:

(1) 利用核物理测量原理γ射线辐射技术, 进行设计, 采用射线准直方式测量, 减少了散射误差影响;它对介质 (特别是固体和浆液介质) 工艺变量的感应和控制能力非常灵敏。采用非接触测量方式, 被测物料液态、固态均可, 稳定可靠, 维护量小, 适用于高温、高压、腐蚀等环境下的密度测量。

(2) 采用铯-137辐射源为发射源, 一般强度为1.85×108Bq, 该射源具有半衰期长 (30年) , 射源在国家标准范围内合理使用, 不会对人体造成伤害。放射源装置和探测器都固定安装在器壁外, 在满足系统需求的同时便于屏蔽防护。

(3) 以微处理器为核心进行运算, 在监测过程中加入温度、放射源衰减补偿, 提高产品的测量精度。通过LED显示仪器的输出信号, 脉冲在线监测及RS485接口, 对工艺过程实行自动控制。

(4) 设计的自动稳峰技术, 集成电路, 以及贴片工艺设计, 可提高提高抗干扰性, 补偿由于温度和元件老化等造成的漂移。

(5) 在信号处理与程序设计过程中, 利用连续滑动平均、比较剔除等措施降低统计涨落对测量精度的影响。

随着人们对γ射线的逐渐认识, 核密度计等先进的仪表将会得到越来越广泛的应用。

参考文献

[1]熊建平.基于单片机的核辐射监测仪[J].核电子学与探测技术, 1999, 19 (1) :63-65.

[2]孙普男.ZHM-1型智能核密度计的研制与应用[J].轻金属, 2002 (9) :57.

[3]姚海峰.同位素密度计的等效标定法及安装中应注意的问题[J].核技术, 1997, 20 (9) :235.

[4]何为民.智能放射性勘察仪器[M].北京:原子能出版社, 1994.

[5]蒋自力.核料位计和核密度计在石油化工装置中的应用[J].石油化工自动化, 2000 (2) :58-62.

基于核密度估计的喷枪沉积建模 篇2

喷涂机器人现已广泛应用于现代制造业,并且工作方式也由人工示教-再现式向离线编程等方式发展。喷枪沉积建模是对涂料离开喷枪后在目标表面沉积成膜的过程进行建模和分析。在喷涂机器人离线轨迹规划中, 喷枪沉积建模是离线编程的前提和重要依据,对最终喷涂质量、效率等具有重要作用。

涂料沉积的过程复杂,一般涂料通过空气喷枪、高压无气喷枪或旋杯等装置雾化成大量液滴,通过气流场喷射到待涂目标表面,部分液滴沉积在目标表面,经流平、 固化,成为涂层;另一部分液滴不能沉积,成为漆灰等废物。由于缺乏准确实用的建模手段,喷涂机器人的离线编程较焊接、搬运、码垛等其他工业机器人更加复杂[1]。

喷涂时涂料沉积过程如下:涂料借助喷枪雾化后, 高速冲击待涂目标表面,部分沉积生成涂层。在喷涂过程中,涂料通过管道输送至喷枪,从喷嘴高速喷出进入空气,其后迅速的减速、雾化成不同尺寸的液滴,直径从几微米到几十微米不等,通常流量的喷枪,每秒产生的液滴数在109量级[2]。当液滴冲击待涂目标表面,液滴要么沉积要么飞溅,这与液滴大小、速度、角度、表面粗糙度以及涂料的粘度等因素有关。

1常规建模方法

建立准确的喷枪沉积理论模型困难。目前,常见的模型主要有无限范围模型和有限范围模型两大类:

1)无限范围模型

无限范围模型多使用于喷涂面为平面的情况,对于曲面的模型精度不佳,典型无限范围模型如高斯分布[3和柯西分布模型[4,5]等,只适应喷枪垂直于工件表面的情况,应用较少。

2) 有限范围模型

有限范围模型更接近实际物理模型,都是以喷枪喷雾形状为圆锥的基础建立的,对于实际喷涂有一定的指导意义;典型的有椭圆分布模型、抛物线分布模型、分段函数模型[6,7],均匀厚度模型[8]、变厚度模型[9]、β分布模型[9]。

为降低喷枪沉积建模过程的计算复杂度,针对液滴的一个有限子集进行仿真。建模方法应具有灵活性、良好的边缘效应、鲁棒性和可接受的计算量。

2基于核密度估计的建模方法

核密度估计(kernel density estimation)是数理统计学中构造未知密度函数的方法。由Rosenblatt(1955) 和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window)。Ruppert和Cline基于数据集密度函数聚类算法提出修订的核密度估计方法。核密度估计方法,不利用有关数据分布的先验知识,对分布函数不附加任何假定,是一种从样本本身出发研究分布特征的方法,因而受到高度重视[10]。

具体如下:

核密度估计是从一组有限样本构建一个未知的概率密度函数估计的方法。从一个未知的概率密度函数测得一组样本X =(X1,X2,...,XN),单变量核估计可表示为:

其中h是一个带宽参数,K是正值、归一化、对称核函数,即:

尽管核函数的选择对建模结果并不是至关重要的, 但更平滑的核函数往往能得到更好的结果。通常人们更愿意采用定义在有限范围的内核函数,比如x定义在 [-1;1]区间上。易见到,这能降低计算量。表1给出了常见的核函数范例。在本文中,我们选择使用Biweight内核由于是无论是在顺利和有限的性能[-1,1]。

带宽参数h的选择对结果同样重要。h值的大小和模型误差相关。目前已有一些能够自动选择带宽的方法, 如cross-validation方法[11],插件法[12]等。这些方法在本文中不做讨论,通过手动选择带宽,实现误差最小化。

喷枪沉积厚度函数可以看作体积为Vi的液滴的密度函数。代入核密度函数,有:

其中,SK是K的归一化常数:

假设如下:X点附近为平面,法线方向为nx。沉积点Xi在表面上被投影沉积方向vi上(如图1所示)。这些投影的表面点通过归一化和带宽选择比较,以确定他们是否应该包括在估计,使用:

通过式(6)和式(7)可以排除喷涂表面背面液滴以及与法线夹角较小液滴对x点沉积的影响。

使用核密度估计算法所得到的单层厚度曲线示意如图2所示,相比经典的直方图方法,本方法所得的结果更平滑。

考察指标:厚度方差:

其中, 表示核密度沉积模型仿真得到的厚度,T表示实际测量得到的厚度函数。需要指出,作为输入,算法生成的各点被取样,因此不能够完全代表厚度的函数。厚度的绝对误差不仅取决于建模方法,同样取决于样本的数量。

核密度估计沉积模型的结果在独立于网格划分的类型和质量,因此可以直接使用CAD网格,从而满足灵活性和鲁棒性的双重要求。

3试验验证

试验设备 采用A B B I R B 5 4 0 0 - 1 2型喷涂机 器人、GRACO AL型喷枪及配套高压泵、马口铁板 (100cm×100cm)、磁性测厚仪(误差±5μm)、某型涂料等。喷涂时分四道水平喷涂马口铁板,每道间隔15cm,如图3所示。喷涂运动速度设置为90cm/s,喷枪距离马口铁板表面距离恒定为12cm,高压泵工作压力150Bar,涂料粘度60s。测量涂层厚度时,测量L线上15个样本点涂层厚度,每个样本点间间隔2cm。

在获得实际喷涂数据的基础上,将仿真核密度模型预测值和测量值进行了比较。仿真计算时,n取120000,带宽h取3mm。结果如图4所示。

从图4中不难发现:建模得到的厚度曲线与测量数据的趋势基本相符,同时,各点测量值与模型期望值误差也在较低水平上。

4结论

在本文中,通过使用核密度估计算法进行喷枪沉积建模,得到一种与现实情况较符合的方法,此种方法不依赖于区块划分,适合直接在CAD中计算估计,不需要任何重新网格化,更适合指导离线轨迹规划。试验结果证明了本文模型的正确性与实用性。

摘要：针对喷枪的涂层沉积建模问题进行研究,通过分析喷涂系统和涂层形成过程,并引入核密度估计方法,构建一种能够适应涂层不同分布情况的喷枪沉积模型,以获得更好的建模精度和应用适应性。开展喷涂工艺试验验证,结果表明,新模型具有高精度和广泛适应性。

核密度计 篇3

在雷达辐射源信号识别中,雷达辐射源基本参数如载频、脉冲重复周期和脉冲宽带都是采用直方图法进行参数分析与特征提取。在直方图法中分组数的确定,一直没有一个确定的规则,影响了直方图的有效性。在分组数适当情况下,直方图的形状应能反映参数的概率密度分布情况[1]。因此,根据参数的概率密度函数来确定直方图分组数是一可行的途径。若参数概率密度函数未知,则可预先采用核密度估计方法估计出样本的概率密度函数。

设X1,X2,…,Xn是从一维总体X中抽出的独立同分布样本,X具有未知的密度函数p(x),x∈R,则p(x)的密度核估计为:

式中n为样本集包含的样本数目;K(u)为核密度估计函数,简称核函数;h称为窗宽,决定了核函数的方差。常用的核函数有Epanechnikov核、三角核、高斯核和均匀核。窗宽h的选取是核函数密度估计能否成功的关键。理论上选择最优窗宽是从密度估计与真实密度之间的误差开始[2,3],如使密度核估计与真实密度函数p(x)之间的均方误差式达到最小的h作为最优窗宽值[4]。

通过极小化MSE(h0),得到最优窗宽:

其中。然而真实密度是未知的,否则无需进行估计。因此,本文先采用改进核密度估计方法,使窗宽的选择不依赖于真实概率密度,而只依赖于给定的核函数和所抽取的样本点;然后根据估计出的密度核估计,来搜索直方图的最优分组数;最后,在最优分组数直方图下进行参数的分析与提取。

1 改进的核密度估计方法

对于来自某未知总体分布的一组样本X1,X2,…,Xn,根据选取核函数的不同,由式(1)给出样本的两个不同的核估计[5]:

和之间的积分平方误差ISE为:

由于与为同一样本的密度核估计,二者之间的误差越小,说明选取的窗宽越好。故最优窗宽应使ISE(h)取最小值。

将式(4)和式(5)代入式(6),通过计算可得:

其中:

改进的核密度估计方法就是在ISE最小准则下,寻求最优窗宽,此时的最优窗宽不是理论上的最优窗宽,而是接近理论最优窗宽的一个值。具体过程如下:

(1)设定初始窗宽h,门限δ和步长Δ;

(2)选取两个不同核函数,由式(1)给出样本的两个不同的核估计与(x);

(3)根据式(6)计算(x)与(x)的积分平方误差ISE(h);

(4)如果ISE(h)>δ,令h=h+Δ,返回步骤(2);否则最优窗宽,执行步骤(5);

(5)取样本的核密度函数为(x)和(x)的平均值。

2 基于改进核密度估计确定最优分组

设T是某一实参数观测数据的集合,区间A是该参数的取值范围,δi(i=1,2,…,N)是区间A的子区间,fi是T集合中落入子区间δi的数值个数,称为频数。其中δi和fi分别满足下面的条件:

(1)对任意i和j(i=1,2,…,N;j=1,2,…,N;i≠j),δ1和δj的交集为空,即δi∩δj;Φ,Φ表示空集;

(2)。

δi(i=1,2,…,N)的宽度可以相等也可以不相等,视具体情况而定。以子区间δi的宽度和fi为边长,做出N个长方形,称为直方。把参数取值范围作为横轴,频数作为纵轴,将每个直方按相应的δi放置在横轴上,这样做出的图形就称为直方图[6]。

将直方图进行归一化处理,令,令归一化处理后直方图的上部轮廓线函数为f(x),如图1所示。参数服从的概率密度函数为p(x),定义f(x)和p(x)的贴近度为:

贴近度越小,说明直方图的上部轮廓线与参数概率密度函数越接近,此时的分组就越好;贴近度越大,说明直方图的上部轮廓线与参数概率密度函数相差较远,此时的分组就不是最优分组。

设参数观测数据为xi(i=1,2,…,n),n代表该批观测数据的总个数,则基于改进核密度估计确定最优分组的直方图算法如下:

(1)采用改进核密度估计算法估计出样本的概率密度函数;

(2)计算样本的平均值和标准差,取参数的下限值可取(xmin-w'×Sx),上限值可取(xmax+w'×Sx),其中xmin、xmax分别为xi(i=1,2,…,n)中的最小值和最大值;

(3)令分组数N=5,作出直方图,将直方图进行归一化处理,令,令归一化处理后直方图的上部轮廓线函数为f(x);

(4)根据式(13)计算f(x)和的贴近度σ,如果σ小于门限σth,则令N=N+2,返回步骤(3);否则此时的分组数N即为最优分组,执行步骤(5);

(5)根据最优分组下的直方图进行参数分析和提取,算法结束。

3 仿真实验

实验一为了验证改进核密度估计算法的有效性,选取n个来自正态分布N(0,2)的样本点即:Xi=normrnd(0,2,1,n),分别利用式(3)计算出理论最优窗宽h0和改进核密度估计算法得到最优窗宽,结果如表1所示。样本集数目为n=1 000时,真实概率密度函数p(x)和核密度估计(x)的曲线如图2所示。

实验二为了验证改进直方图算法的有效性,将其用于分析雷达辐射源信号的基本参数脉冲重复周期(PRI)。设定重复周期固定类型的取值为50ms;重复周期抖动类型的中心值为50ms;重复周期三参差类型的取值分别为50ms、100ms、150ms;重复周期捷变或滑变取值分别为40ms、42.5ms、45ms、47.5ms、50ms、52.5ms、55ms、57.5ms、60ms。每种类型仿真500个样本,并设定1%的测量误差。

实验中取w'=0.25,σth=0.8。直方图法中的分组有等距分组和不等距分组,实验采用等距分组,将直方图的分组数从小到大进行搜索,当归一化处理后直方图的上部轮廓线函数与估计出概率密度函数之间的贴近度小于门限σth时,搜索停止。此时,计算出的贴近度与分组数如表2所示,概率密度曲线及最优分组下的直方图如图3至图6所示。从图中可以看出,在取最优分组数时的直方图可以清晰地反映出参数工作类型和取值情况。

4 结语

为解决直方图法中的最优分组问题,本文采用改进的核密度估计算法来估计参数的概率密度函数,该方法得到窗宽可以接近于理论上的最优窗宽。然后,引入贴近度的概念,将参数直方图的上部轮廓线与参数服从的概率密度函数之间的接近程度作为最优分组的判决准则。利用改进核密度估计确定最优分组的直方图算法分析雷达辐射源信号基本参数之一PRI,结果表明,该方法可以定量获得最优的分组数,在此过程中不需要人参与判断搜索过程是否结束,从而实现算法的自动化。

参考文献

[1]刘甸瑞.作直方图中的最优分组问题[J].物探化探计算技术, 1995,17(1):62-67.

[2]Ibrahim A Ahmad,Fan Yanqin.Optimal Bandwidth for Kernel Density Estimators of Functions of Observations[J].Statistics&Probability Letters,2005,51(3):245-251.

[3]摆玉龙,杨志民.基于Parzen窗法的贝叶斯参数估计[J].计算机工程与应用,2007,43(7):55-58.

[4]黎运发,黄名辉.核密度估计逐点最优窗宽选择的改进[J].统计与决策,2011(14):28-32.

[5]郭照庄,霍东升,孙月芳.密度核估计中窗宽选择的一种新方法[J].佳木斯大学学报:自然科学版,2008,26(3):401-403.

核密度计 篇4

点云配准是计算机视觉领域的基本问题之一,广泛应用于距离数据融合、医学图像对准、目标定位、跟踪及识别[1]。ICP算法是经典的点云配准算法,然而,ICP算法存在如下不足:①代价函数的不可微性要求点云之间有足够的重叠区域,否则算法容易陷入局部极小;②ICP算法需要一个紧的初始条件。上述两点要求在应用ICP算法之前必须通过一定的预处理方法估计点云之间的重叠区域,并进行粗配准,以获得ICP算法的初始值。点云配准的难点在于异常值、背景和噪声等对配准精度的影响。为了消除上述因素对配准算法精度和鲁棒性的影响,文献[2,3]相继对ICP算法进行了改进。除此之外,人们还提出了两类方法对配准算法的鲁棒性进行改善。第一类方法通过对异常值、背景和噪声进行统计学建模,进而采用ML/MAP估计的方法进行配准。如Hasler等[4]研究了图像配准中的异常值建模问题;Ying等[5]研究了基于统计模型的局部重叠对象识别问题。然而,这类方法的问题在于,在许多实际应用场合,对干扰因素的精确统计学建模是不可行的。第二类方法采用变换不变的鲁棒特征结合代价函数进行配准,以获得针对异常值、背景和噪声等干扰鲁棒性更好的配准算法。如Viola等[6]提出了一种基于融合信息的配准算法;Tsin等[7]提出了一种基于统计相关的鲁棒点云配准算法(KC算法);Huttenlocher等[8]提出了一种基于局部Hausdorff距离的配准算法;Chui等[9]提出了一种基于混合模型的特征配准算法;Jian等[10]提出了一种基于高斯混合模型的鲁棒配准算法(GMM(gaussian mixture model)算法)等。杨静等[11]提出了一种基于简单Schur凹函数的图像配准测度;李晖等[12]提出了塔式分解和模糊梯度场的图像配准算法。这类算法的问题在于不同的代价函数和优化方法对算法性能影响较大。因此,设计合适的代价函数并根据代价函数的特点探寻全局合适的优化方法成为此类算法的关键,本文即属于此类方法。

本文基于核密度估计理论提出了一种新的测度函数,对不同参数下测度的性能进行分析,论证了利用该测度进行寻优配准时存在的局部极值现象和极值漂移现象。在此基础上提出了一种基于BFGS拟牛顿法的变尺度优化配准算法,解决了算法鲁棒性与配准精度之间的矛盾。

1 点云的核密度分布模型建立

设M={mi,i=1,2,…,Nm}表示模型点云,S={si,i=1,2,…,Ns}表示场景点云,两个点云可以视为对目标所对应的实向量空间Rd进行观测所得到的样本。

设X为d维随机变量,S为对随机分布X进行独立观测得到的随机样本集,$\widehat{f}$X(x|S)表示由样本集S估计的随机变量X的概率密度函数,则由核密度估计理论,$\widehat{f}$X(x|S)可以表示为[13]

$\widehat{f}{}_{\mathit{X}}(\mathit{x}|\mathit{S})=\frac{1}{{\rm N}{}_s|\mathit{{\rm H}}|}{\displaystyle \sum _{i=1}^{{\rm N}{}_s}{\rm K}}(\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{s}{}_i))$ (1)

其中,H为d×d阶正定带宽矩阵,为简化计算,取H=hId,Id为d×d阶单位矩阵;K:Rd→R为满足如下条件的非负核函数:

(2)

对于三维点云,d=3。出于一般性考虑,选取高斯函数作为概率估计的核函数,即

${\rm K}(\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{s}{}_i))=\exp [-\frac{(\mathit{x}-\mathit{s}{}_i){}^{\text{Τ}}\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{s}{}_i)}{2}]$ (3)

由此,可以通过上述过程分别建立离散点云M和S的连续核密度分布模型。

2 相似性测度的构造

核密度函数相似性测度对配准算法性能的提高具有重要作用,通常用于描述两个核密度分布相似性的测度有积分绝对值误差测度、Kullback-Liebler测度、MI(mutual information)测度、欧氏测度等[10]。不同的测度函数优化过程具有不同的鲁棒性和计算效率,如欧氏距离计算方便、意义明确,在相似性评价方面得到了广泛应用,但相比于Kullback-Liebler测度,其鲁棒性稍差。在文献[10]提出的GMM测度函数的基础上,本文提出了如下测度函数:

$d{}_{\alpha }(\widehat{f},\widehat{g})={\displaystyle \int (}\frac{1}{\alpha }\widehat{f}{}^{1+\alpha }+\alpha \widehat{g}{}^{1+\alpha }-\frac{1}{\alpha }\widehat{f}\widehat{g}{}^{\alpha }-\alpha \widehat{f}{}^{\alpha }\widehat{g})\text{d}\mathit{x}$ (4)

当α→0时,有

$\begin{array}{l}d{}_0(\widehat{f}‚\widehat{g})={\displaystyle \int \underset{\alpha \to 0}{{\displaystyle \lim }}}(\frac{1}{\alpha }\widehat{f}{}^{1+\alpha }+\alpha \widehat{g}{}^{1+\alpha }-\\ \frac{1}{\alpha }\widehat{f}\widehat{g}{}^{\alpha }-\alpha \widehat{f}{}^{\alpha }\widehat{g})\text{d}\mathit{x}={\displaystyle \int \widehat{f}}\ln \frac{\widehat{f}}{\widehat{g}}\text{d}\mathit{x}(5)\end{array}$

此时,d0($\widehat{f}$,$\widehat{g}$)为$\widehat{f}$与$\widehat{g}$之间的Kullback-Leibler距离。

当α→1时,有

$d{}_1(\widehat{f},\widehat{g})={\displaystyle \int (}\widehat{f}-\widehat{g}){}^2\text{d}\mathit{x}$ (6)

此时,d1(f,g)为$\widehat{f}$与$\widehat{g}$之间的欧氏距离。因此,利用本文提出的测度函数可以通过单一参数α可以调整算法的鲁棒性和渐进效率。

当配准对象为刚体时,空间变换矩阵T可以表示为:T(M)=Rmi+t,(i=1,2,…,Nm),其中R为3×3阶旋转矩阵,t为3×1阶平移向量。此时,

$\begin{array}{l}{\displaystyle \int \widehat{f}}{}_{\mathit{X}}^{1+\alpha }(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\text{d}\mathit{x}=\\ {\displaystyle \int [}\frac{1}{n|\mathit{{\rm H}}|}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}{\rm K}}(\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{R}\mathit{m}{}_i-\mathit{t}))]{}^{1+\alpha }\text{d}\mathit{x}(7)\end{array}$

令x=Rx′+t,则

$\begin{array}{l}{\displaystyle \int \widehat{f}}{}_{\mathit{X}}^{1+\alpha }(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\text{d}\mathit{x}=\\ {\displaystyle \int [}\frac{1}{n|\mathit{{\rm H}}|}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}{\rm K}}(\mathit{{\rm H}}{}^{-1}\mathit{R}({\mathit{x}}^{\prime }-\mathit{m}{}_i))]{}^{1+\alpha }\text{d}(\mathit{R}{\mathit{x}}^{\prime }+\mathit{t})(8)\end{array}$

由式(3)可得

$\begin{array}{l}{\rm K}(\mathit{{\rm H}}{}^{-1}\mathit{R}({\mathit{x}}^{\prime }-\mathit{m}{}_i))=\\ \exp [-\frac{(\mathit{x}-\mathit{m}{}_i){}^{\text{Τ}}\mathit{R}{}^{\text{Τ}}\mathit{{\rm H}}{}^{-1}\mathit{R}(\mathit{x}-\mathit{m}{}_i)}{2}]=\\ \exp [-\frac{(\mathit{x}-\mathit{m}{}_i){}^{\text{Τ}}\mathit{R}{}^{\text{Τ}}\mathit{R}(\mathit{x}-\mathit{m}{}_i)}{2h}](9)\end{array}$

由于R为刚体变换矩阵,故,RRT=I3,detR=1,则

$\begin{array}{l}{\rm K}(\mathit{{\rm H}}{}^{-1}\mathit{R}({\mathit{x}}^{\prime }-\mathit{m}{}_i))=\\ \exp [-\frac{(\mathit{x}-\mathit{m}{}_i){}^{\text{Τ}}(\mathit{x}-\mathit{m}{}_i)}{2h}]={\rm K}(\mathit{{\rm H}}{}^{-1}({\mathit{x}}^{\prime }-\mathit{m}{}_i))(10)\end{array}$

将式(10)代入式(8)可得

$\begin{array}{l}{\displaystyle \int f}{}_{\mathit{X}}^{1+\alpha }(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\text{d}\mathit{x}=\\ {\displaystyle \int [}\frac{1}{n|\mathit{{\rm H}}|}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}{\rm K}}(\mathit{{\rm H}}{}^{-1}\mathit{R}[{\mathit{x}}^{\prime }-\mathit{m}{}_i)]{}^{1+\alpha }\det \mathit{R}\text{d}{\mathit{x}}^{\prime }(11)\end{array}$

式(11)表明,表达式∫f1+αX(x|T(M))dx与配准参数无关,故为了提高计算效率在刚体点云配准中该项予以忽略。又因为${\displaystyle \int \alpha }\widehat{g}{}^{1+\alpha }(\mathit{x})\text{d}\mathit{x}$也与配准参数无关,故对刚体对象配准测度函数可表示为

当α=1时,测度函数F(T)的解析表达式为

$\begin{array}{l}F{}_1^1({\rm T})=-{\displaystyle \int \widehat{f}}{}_{\mathit{X}}(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\widehat{f}{}_{\mathit{X}}(\mathit{x}|\mathit{S})\text{d}\mathit{x}=\\ \frac{1}{{\rm N}{}_s{\rm N}{}_m|\mathit{{\rm H}}|{}^2}{\displaystyle \int {\displaystyle \sum _{i=1}^{{\rm N}{}_m}\exp }}[-\frac{(\mathit{x}-\mathit{{\rm T}}(\mathit{m}{}_i)){}^{\text{Τ}}\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{{\rm T}}(\mathit{m}{}_i))}{2}]\cdot \\ {\displaystyle \sum _{j=1}^{{\rm N}{}_s}\exp }[-\frac{(\mathit{x}-\mathit{s}{}_j){}^{\text{Τ}}\mathit{{\rm H}}{}^{-1}(\mathit{x}-\mathit{s}{}_j)}{2})\text{d}\mathit{x}=\\ \frac{\sqrt{\text{π}}h{}^{-\frac{3}{2}}}{{\rm N}{}_s{\rm N}{}_m}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}{\displaystyle \sum _{j=1}^{{\rm N}{}_s}\exp }}[-\frac{\parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{s}{}_j\parallel {}^2}{4h}](13)\end{array}$

由于系数$\frac{\sqrt{\text{π}}h{}^{-\frac{3}{2}}}{{\rm N}{}_s{\rm N}{}_m}$与配准过程无关,式(13)可等价于

$F{}_1^2(\mathit{{\rm T}})={\displaystyle \sum _{i=1}^{{\rm N}{}_m}{\displaystyle \sum _{j=1}^{{\rm N}{}_s}\exp }}[-\frac{\parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{s}{}_j\parallel {}^2}{4h}]$ (14)

当α=2时,测度函数F(T)的解析表达式为

$\begin{array}{l}F{}_2(\mathit{{\rm T}})={\displaystyle \int [}\frac{1}{2}\widehat{f}{}_{\mathit{X}}(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\widehat{f}{}_{\mathit{X}}^2(\mathit{x}|\mathit{S})+\\ 2\widehat{f}{}_{\mathit{X}}^2(\mathit{x}|\mathit{{\rm T}}(\mathit{{\rm M}}))\widehat{f}{}_x(\mathit{x}|\mathit{S})]\text{d}\mathit{x}=\\ \sqrt{\frac{2\text{π}}{3h{}^3}}\frac{1}{{\rm N}{}_m{\rm N}{}_s}\{\frac{1}{2{\rm N}{}_s}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}}{\displaystyle \sum _{j=1}^{{\rm N}{}_s}}{\displaystyle \sum _{k=1}^{{\rm N}{}_s}}\\ \exp [-(\parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{s}{}_j\parallel {}^2+\parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{s}{}_k\parallel {}^2+\\ \parallel \mathit{s}{}_j-\mathit{s}{}_k\parallel {}^2)/6h]+\frac{2}{{\rm N}{}_m}{\displaystyle \sum _{i=1}^{{\rm N}{}_m}}{\displaystyle \sum _{j=1}^{{\rm N}{}_m}}{\displaystyle \sum _{k=1}^{{\rm N}{}_s}}\\ \exp [-(\parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{{\rm T}}(\mathit{m}{}_j)\parallel {}^2+\\ \parallel \mathit{{\rm T}}(\mathit{m}{}_i)-\mathit{s}{}_k\parallel {}^2+\parallel \mathit{{\rm T}}(\mathit{m}{}_j)-\mathit{s}{}_k\parallel {}^2)/6h]\}(15)\end{array}$

3 相似性测度的性能分析

首先,随机生成一个3×100的服从正态分布的随机矩阵作为场景点云M。其次,对场景点云M沿x轴方向平移3个单位,并随机选取其中50列作为模型点云S。h分别取0.1,10,20,50,100时,测度函数随偏移量tx的变化关系如图1所示(说明:本文除角度的量纲取(°)或rad,其他均为量纲一量)。再采用同样的方法测试测度函数对旋转参数的性能,设定模型点云和场景点云之间的旋转关系为θx=π/3。此时,测度函数随旋转量θx的变化关系如图2所示。

1. h=0.1 2.h=10 3.h=20 4.h=50 5.h=100

1. h=0.1 2.h=10 3.h=20 4.h=50 5.h=100

由式(14)可知,α=1的情况下本文测度等价于文献[7]和文献[14]所采用的测度。图1表明,针对平移变换,通过搜索测度函数的极大值可以确定平移参量,实现点云的平移配准。图1a表明:α=1时,尺度参数越大测度函数曲线越宽,曲线越光滑,峰值位置偏离真实值,采用大尺度参数的测度函数进行配准时将存在残余误差,我们把这种现象称为过尺度参数导致的极值漂移现象;尺度参数越小测度函数曲线越窄,曲线可能出现局部极值点,采用小尺度参数的测度函数进行配准时算法可能因陷入局部极小而导致失配,我们把这种现象称为欠尺度参数导致的局部极值现象。图1b表明,当选取α=2时,测度函数曲线在极值漂移和局部极值方面性能更好,采用此测度函数的配准算法具有更好的鲁棒性和配准精度。图2表明,在进行旋转参数的配准时,选取参数α=1也会导致极大值漂移和局部极值的问题,而选取参数α=2可以改善这两种现象,提高算法的鲁棒性和准确性。

为了检验测度函数极值点的统计稳定性,将实验重复进行100次,统计各种测度函数全局极大值所在位置的均值μ和标准差σ,实验结果如表1和表2所示。

对比三种测度可以发现,与MI测度相比,本文提出的测度和GMM测度对于平移参数和旋转参数的估计均具有更好的准确性和有效性。通过对比GMM测度实验结果和本文测度实验结果可以发现,当α=1时,本文测度的估计性能与GMM测度估计性能相同,而当α=2时,测度对参数估计的性能优于GMM测度的估计性能。

4 基于BFGS优化的变尺度配准算法

通过表1和表2可以发现,应用本文提出的测度可以改善配准参数估计的准确性和有效性。然而,对比不同尺度参数下GMM测度和本文测度的性能可以发现,随着尺度的增大,对平移参数和旋转参数估计的偏倚现象明显,导致参数估计准确度降低。图1和图2也表明,尺度参数过小可能会导致局部极值现象的产生,从而导致点云失配。为此,本文提出了一种基于BFGS优化的变尺度配准算法,其优化过程如图3所示。

首先,随机选取初始参数T0作为极值点搜索的初始条件,同时,为了保证算法的收敛域和平滑性,选取一个较大的初始参数h0。再利用BFGS拟牛顿法搜索该尺度下测度函数极大值所在位置T1,在通过某种机制减小尺度参数h并以T1为搜索初始条件,获得T2依次类推,直至达到预定的配准精度。

基于BFGS优化方法的变尺度配准算法流程如下:

(1)初始化。h0=3λmax和T0=[0 0 0 0 00]T,并令k=0,其中,λmax为两个点云自相关矩阵的最大特征值,I4表示4×4阶单位矩阵。

(2)更新变换矩阵T。以Tk为初始条件,以hk为尺度,利用BFGS拟牛顿法搜索测度函数F(Tk,hk)的极大值。并将此极大值所在位置记为Tk+1。

(3)更新尺度参数h。hk+1=μ hk,选取μ∈(0.9,0.96),k←k+1。

(4)循环步骤(2)、(3),直至‖Tk+1-Tk‖<ε。

采用这种衰减尺度的BFGS拟牛顿搜索法进行配准主要有三个优点:①拓展了算法的收敛域;②克服了欠尺度参数导致的局部极值现象,减小了算法失配的可能;③克服了过尺度参数导致的极大值漂移现象,提高了配准精度。

5 实验研究

为了验证算法的性能,本节分别通过仿真信号和实测信号的配准实验对算法的可靠性和鲁棒性进行分析和比较。

5.1 受白噪声干扰点云配准实验(实验1)

(1)模型点云M的构造。通过如下公式构造一条三维空间曲线,并在该曲线均匀选取200个点作为模型点云M:

$\{\begin{array}{l}x=\sin t\\ y=\cos tt=0~10\text{π}\\ z=t\end{array}$

(2)场景点云S的构造。对模型点云M进行旋转和平移变换,变换参数

$[\theta {}_x\theta {}_y\theta {}_{z}t{}_{x}t{}_{y}t{}_z]{}^{\text{Τ}}=[\frac{\text{π}}{10}000.100]$

再增加信噪比为20dB的白噪声作为干扰,由此构造场景点云S。原始点云以及通过ICP、GMM和本文方法的配准结果如图4所示。

图4表明,在给定实验条件下ICP算法、GMM算法和本文算法均可实现点云的精确配准,说明这几种方法相对与白噪声干扰都具有鲁棒性。

5.2 配准算法收敛性能实验(实验2)

为了检验算法的收敛性能,设计了实验2,实验中模型点云M的生成方法和实验1相同。在此基础上分别改变绕x轴方向转角θx∈(-π,π)和沿x方向平移tx∈(-20,20),并附加20dB白噪声作为场景点云S。分别用ICP算法,GMM算法和本文算法(α=1,μ=0.9)进行配准,求取各种算法的配准误差。配准误差随变换参数变化关系曲线如图5所示。

1. ICP 2.GMM 3.OUR

观察图5a可以发现,在针对单纯平移参数的配准中,ICP算法和本文算法在(-20,20)范围内可获得全局收敛的配准效果,而GMM算法的收敛范围为(-4.3,4.3)。这说明在处理单一平移参数的配准时,GMM算法要求两个点云之间的相对平移量较小,限制了算法的实际应用。由图5b可知,针对单纯旋转变换,ICP算法的收敛区间为(-85°,85°),GMM算法的收敛区间为(-90°,105°),而本文算法的收敛区间为(-180°,170°)。收敛区间越宽,说明使用该算法时两个点云之间的相对变换参数的允许范围就越大,从而算法的适用性就越好。综合分析图5可以发现,不论是针对旋转参数的配准还是针对平移参数的配准,相对于ICP算法和GMM算法,本文提出的算法都具有更好的收敛性能。

5.3 局部重叠点云配准实验(实验3)

为了验证点云局部重叠情况下本文算法与传统算法的性能,设计了实验3。实验原始数据为一组人脸激光扫描散乱点云,记为PtsFace,扫描点数NPtsFace=392,点云PtsFace的坐标分布范围如下:(x×y×z)∈([-1.18,1.11]×[-0.58,1.15]×[-1.82,2.56])。取散乱点PtsFace的前80%的点作为模型点云M。再对PtsFace进行变换,参数[θxθyθztxtytz]T=[π/5 0 00.2 0 0]T,取变换后点云的后70%作为场景点云S。实验结果如图6所示。

图6表明,在给定实验条件下,ICP算法陷入局部极小,不能实现给定点云的配准。图6c表明,当参数选取为σ=0.56时,GMM算法存在残余误差。由前面的理论分析可知,该残余误差是由过尺度参数引起的极值点漂移引起的,不能通过增加迭代次数和降低门限值的方法加以改善。

图6d表明,当参数选取为σ=0.02时,GMM算法配准失败,由上文分析可知,这是由于欠尺度参数导致的局部极值现象所引起的,而这种现象主要存在于点云之间存在非重叠区域的情况。然而,在实际应用中点云之间往往存在明显的非重叠区域。因此,GMM算法在实际应用中也受到限制。图6e和图6f表明,本文方法可以实现具有非重叠区域点云的精确配准。

6 结束语

针对传统点云配准算法鲁棒性差、收敛区间窄的缺点,基于核密度估计的思想,提出了一种新的鲁棒的点云配准算法。本文从点云的概率分布模型建立入手,将点云配准问题表述为核密度函数相似性寻优问题,提出了用于描述核密度分布相似性的测度函数,研究了该测度函数的性能及其与其它测度函数之间的关系。针对算法实际应用中过尺度参数和欠尺度参数两种情况造成的配准残余误差和配准失败的情况,提出了一种变尺度BFGS寻优方法,实现了点云的鲁棒、精确配准。

核密度计 篇5

随机潮流（probabilistic load flow,PLF）计及电力系统运行中的不确定性，可综合分析系统的薄弱节点与运行风险，是静态安全分析的有力工具[1,2,3,4,5,6,7]。

蒙特卡洛法，又称为随机抽样方法或统计试验方法，目前在随机潮流计算中应用广泛[5,6,7]。随着研究的深入，蒙特卡洛的扩展问题凸显出来。扩展蒙特卡洛方法可将随机潮流的误差转化为可控量，并在扩展样本规模时保留已知的潮流计算结果。简单随机抽样的样本之间相互独立，实现扩展较为方便，但其计算效率较低，获得精确结果的计算代价大[8]。以拉丁超立方抽样[5,6,7]为代表的伪随机采样方法根据输入变量的累积分布函数进行分层采样，可提高蒙特卡洛的精度与计算效率。针对其扩展问题，文献[9]提出扩展拉丁超立方抽样，目前已应用于随机潮流计算[7]，但样本仅能由N扩展到2 N，因此无法兼顾算法效率与扩展规模。此外，由于该方法不能保证序列的低差异性，因此文献[10]指出这限制了蒙特卡洛的计算精度及收敛速度，而以低差异序列（low-discrepancy sequence,LDS）为基础的准蒙特卡洛方法在精度和收敛速度上均优于基于拉丁超立方抽样的方法，可以任意维度扩展样本序列，已应用于电子电路设计中[10]。而在低差异序列中，Sobol序列可由布尔代数方法快速生成，故应用最广泛。

除了样本产生方法，扩展蒙特卡洛方法的另一问题是算法收敛判据的选取。基于简单随机抽样的蒙特卡洛方法中，由于样本之间的独立性，一般可采用方差系数作为收敛判据[11]。为弥补其在扩展拉丁超立方抽样方法中应用的缺陷，文献[7]中提出改进方差系数定义随机潮流的收敛判据，获得了置信概率较高的结果。但由于方差系数推导的假设前提是输出变量为正态分布，因此这类判据在输出变量偏离正态分布时会带来较大误差，影响扩展蒙特卡洛方法的准确收敛。

本文针对扩展拉丁超立方抽样方法的缺陷，提出扩展准蒙特卡洛方法，可克服蒙特卡洛方法收敛性的瓶颈，并能够实现任意步长的扩展，通过算例分析证明了其在随机潮流计算中的正确性与高效性。此外，针对方差系数应用于非正态变量时的问题，本文提出基于非参数核密度估计结果的收敛判据，改进了现有判据的不足，其有效性通过算例分析得以证明。

1 基于扩展准蒙特卡洛方法的随机潮流

1.1 准蒙特卡洛方法及其误差

当节点注入功率向量S的分布特性已知时，节点电压向量V和支路潮流向量B的概率分布如式（1）所示。

式中：h（·）为概率密度函数分别为节点注入功率对节点电压和支路潮流的Jacobian矩阵。

蒙特卡洛方法利用经验概率分布函数代替h(S），因此可得输出变量概率分布如式（2）所示。

式中：Si为根据节点注入功率分布产生的序列；n为蒙特卡洛的模拟次数；δ（S-Si）为注入功率样本Si的Dirac测度[12]。

为便于分析，将式（2）简化为式（3):

式中：Qn为n次模拟得到的数值结果；pi表示式（2）中的

根据大数定律，Qn会以概率1收敛到输出变量的真实分布函数Q[12]，即

大数定律的意义在于，其指出蒙特卡洛方法随着采样规模无限增大时收敛的必然性，但却不能给出收敛速度的信息。而文献[13]证明，序列的差异性（discrepancy）是影响蒙特卡洛收敛速度的主要因素。

定理1[13]蒙特卡洛方法的误差ε可定义为：

式中：Dn*为L∞-差异性[13],Dn*越小，代表序列在空间中分布更均匀；C为由系统本身性质决定的常数[14]。

由式（5）可知，当系统参数确定时，蒙特卡洛方法的效率取决于序列的差异性。文献[15]给出拉丁超立方抽样等抽样方法产生s维均匀分布序列的Dn*为：

结合式（5），基于拉丁超立方抽样的蒙特卡洛方法收敛速度为。文献[16]指出，任意维度均匀分布序列最小Dn*如式（7）所示。满足该式的序列称为低差异序列，而基于低差异序列的蒙特卡洛方法称为准蒙特卡洛方法。

对比式（6）和式（7）可知，序列的L∞-差异性由变为O(n-1），再结合式（5），引入低差异序列后蒙特卡洛方法的收敛于真实值的速度由变为O(n-1），从而大大加快蒙特卡洛方法的收敛。

1.2 准蒙特卡洛的扩展

在准蒙特卡洛方法中，Sobol序列可以利用布尔代数方法快速生成[17]，过程参见附录A。本文在此基础上提出其扩展方法，基本思路是已有规模为N的满足独立均匀分布的样本矩阵Uold，通过新增加规模为Nstep的样本矩阵Uadd，获得拥有N+Nstep个样本的满足独立均匀分布的样本矩阵Unew。过程如下。

已知Uold、原始多项式P(x）、方向数v=[v1,v2，…，vj](j=1,2，…，N）。

1）采用生成Uold相同的多项式：

2）按照递推关系继续生成Nstep个方向数vj(N<j≤N+Nstep):

式中：表示异或运算。

3）由此可生成Nstep个Sobol序列的新元素：

式中：下标ci表示N+k的二进制表示“…b3b2b1”中最右非零位的下标；Uadd,k和Uold,k表示Uadd和Uold中的第k个样本。

二维Sobol序列长度由100扩展到150的效果分别如图1(a）和图1(b）所示。可以看出，扩展后的Sobol序列保持了序列的差异性，因此不会影响算法的效率。对比文献[7]中的扩展拉丁超立方抽样方法，本文提出的扩展准蒙特卡洛有以下优势：(1)扩展方便，无需额外运算，因此效率高；(2)生成样本均匀性好，可提高随机潮流效率；(3)可实现任意维度扩展，因此可兼顾算法效率与计算精度。

1.3 输入样本生成方法

文献[18]中的Nataf变换可将均匀分布的输入变量转化为任意分布，过程见附录A。在输入变量相关性的处理上，由于输入变量相关系数矩阵可能为非正定或非满秩阵，此时其Cholesky分解不存在，无法采用文献[18]的方法对样本进行排序。而由于相关系数矩阵均为对称阵，存在奇异值分解。针对这一问题，本文提出采用奇异值分解的方法改进文献[18]的相关系数矩阵处理方法，本文通过证明定理2说明奇异值分解排序方法的过程。

定理2设K为n维独立标准正态分布向量，矩阵由向量z的相关系数矩阵ρz的奇异值分解生成：

式中：Uρz为酉矩阵；Λρz为由矩阵奇异值构成的对角矩阵，并按从大到小排列。

则由式（12）定义的z*相关系数矩阵等于ρz:

证明：对，其期望值E(zi*）满足式（13）。

相关系数矩阵为：

因此，可根据定理2的步骤生成相关系数为ρz的正态分布向量后可再通过Nataf反变换，得到服从任意分布的随机向量。

2 基于非参数核密度估计方法的收敛判据

扩展随机潮流的另一个问题是其收敛判据的选取。基于方差系数的扩展随机潮流收敛判据以正态随机变量为假设，当实际情况偏离假设时会存在较大误差，这将影响扩展随机潮流的准确收敛。非参数核密度估计对数据分布不附加任何假定，直接根据原始数据构造其分布。本文提出基于此方法的收敛判据，进一步拓宽扩展随机潮流的适用范围。简要介绍如下。

设Sn={X1,X2，…，Xn}为原始数据，其满足未知分布f。文献[19,20,21]提出通过求解原始数据构造的傅里叶热传导方程如式（15）所示，可得到连续的分布函数估计式g，估计误差可由如式（16）所示的均方积分定义，称为均方积分误差（mean integrated squared error,MISE），本文表示为M。

为得到符合条件的解，必须事先给定方程的初始值和Neumann边界条件。其中初始值由数据集Sn的经验概率密度函数D(x）给定，如式（17）所示：

式中：δ（x-Xi）为数据集中元素Xi的Dirac测度；L为数据集维度。

Neumann边界条件可由式（18）给出。

基于以上条件，式（15）的解析表达式如式所示：

式中：x∈[0,1]；为均值等于2k±Xi、方差为t的高斯分布。方差t又称为尺度参数，是式中唯一的变量，需要通过选择合适的尺度参数使得式（16）中的M最小化。

通常M求解复杂，不适合于快速计算。本文提出采用M的一阶近似M(1）衡量单个变量的精度，并定义其平均值CM(1）为随机潮流的收敛判据，当其小于阈值ε时算法停止：

式中：;m为输出变量维度。

文献[21]中给出了求解t最优值topt的解析方法，如式（22）所示：

‖f″‖2中f未知，无法直接求解。文献[21]提出了近似方法可求解‖f″‖2。

3 扩展随机潮流方法

扩展准蒙特卡洛方法的流程图如图2所示。蒙特卡洛方法可以任意步长扩展，设初始样本规模为N0，最大扩展次数为Nmax，扩展步长为Nstep。

4 算例分析

4.1 算例设置

本文采用IEEE 30和IEEE 118节点系统验证所提算法的计算效果。仿真平台基于Intel Core dual i7-3820 3.6GHz,RAM 8GB。假设风速满足威布尔分布W(c,k)=W(10.7,3.97），风电场功率因数恒定，有功出力与风速关系满足：

式中：PWR为风电场额定功率；vci,vr,vco分别为切入风速、额定风速和切出风速，大小分别为2.5m/s,13m/s和25m/s。

对于IEEE 30节点系统，节点3,4,6,28分别接入容量为15MW的风电机组；对于IEEE 118节点系统，节点93,94,95,96分别接入200 MW的风电机组，风电机组出力的互相关系数均为0.9。负荷均为恒功率因数，有功负荷服从期望为额定有功功率，标准差为期望5%的正态分布。互相关系数均为0.5。发电机有功、无功调度方法见附录A图A1。

4.2 收敛判据变化趋势

为分析准蒙特卡洛方法对随机潮流的全局加速性能，IEEE 30和IEEE 118节点系统CM(1）随着采样规模增大的收敛趋势分别如图3(a）和图3(b）所示，N0=25。本节中，为方便与扩展拉丁超立方方法比较，采用加倍扩展的方式。

分别采用扩展拉丁超立方抽样[7]方法和本文方法分析特定变量的变化趋势，部分IEEE 30节点系统支路有功潮流及IEEE 118节点系统的线路无功潮流M(1）随样本规模扩大的变化分别见表1至表4。表中：MP(1）为有功潮流部分M(1);MQ(1）为无功潮流部分M(1）。

综合图3、表1至表4可得到以下结论：(1)由图3可知，采用准蒙特卡洛方法可以提高随机潮流的全局收敛速度；(2)准蒙特卡洛方法可以提高特定变量的收敛速度，同样本规模下可得到精度更高的结果；(3)扩展随机潮流不会影响最终的结果，两种方法所得结果相近。

4.3 收敛判据的应用分析

扩展随机潮流无需事先给定样本规模，可根据误差情况决定算法是否终止，具有较高灵活性。本文以简单随机抽样（SRS)50 000次的结果作为标准，采用式（21）作为收敛判据，设定ε=0.05。对本文提出的扩展准蒙特卡洛方法和扩展拉丁超立方抽样方法的算法效率进行比较。参数及方法组合设置如下：本文方法和扩展拉丁超立方抽样方法[7]起始步长均为100，本文方法扩展步长为100，扩展拉丁超立方抽样方法为规模加倍扩展。

采用IEEE 30节点系统和IEEE 118节点系统分别进行计算，算法计算情况见表5。

由表5中结果可知，本文提出的扩展准蒙特卡洛方法的效率高于扩展拉丁超立方抽样方法，算法收敛用时少，这一结果与4.2节中结论相吻合，在此证明了准蒙特卡洛方法的高效性。

计算收敛时IEEE 30节点系统支路6-7有功潮流和IEEE 118节点系统支路15-19无功潮流收敛时的概率分布分别如图4和图5所示。对比图4和图5可见，三条曲线较为接近。这不仅可以说明本文定义的收敛判据可以有效指导扩展随机潮流的收敛，得到较精确的结果，而且也说明了本文判据的适用性。

为进一步说明本文判据相对文献[7]方差系数的优势，利用本文方法分别采用两种判据应用于IEEE 118节点系统，结果如图6所示。此处本文方法采用规模加倍的扩展方式，方差系数置信概率为90%。

从图6中可以看出，方差系数作为收敛判据有明显的“早熟”现象，收敛判据下降快，在还未到达精度要求时判据数值已下降到很低的水平，这对准确指导算法收敛不利。而本文判据随着样本规模的扩大下降均匀，可较好地指示扩展随机潮流的收敛。此外，本文的判据还具有如下优势：(1)本文判据对输出变量分布没有附加条件，当输出变量偏离正态分布时不会带来误差；(2)不同于方差系数这类概率型指标，本文的判据是基于确定型的，不存在置信区间，因而没有误判的可能。

5 结论

1）本文提出的基于Sobol序列的扩展准蒙特卡洛方法保证了生成样本的差异性，因此保持了准蒙特卡洛方法的高效性，分析证明本文方法效率高于扩展拉丁超立方抽样的方法，同样本规模下可获得更精确的结果。

2）针对输入变量相关性矩阵非正定的情况，本文提出并证明基于奇异值分解的改进方法。

3）两种扩展方法最终均得到相似的结果，因此扩展技术不会影响随机潮流的精度。

4）本文提出的基于非参数核密度估计的收敛判据可以有效指示扩展随机潮流的收敛，可获得较精确的概率特性。

5）本文提出判据对输入变量概率分布没有附加条件，具有直观性和普适性，分析证明该判据没有方差系数的“早熟”现象，因此适合于作为扩展随机潮流的判据。

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx）。

摘要：扩展蒙特卡洛方法将随机潮流的误差转化为可控量,并可在样本数增加时保留已知的潮流计算结果。基于扩展拉丁超立方抽样的方法效率比简单随机抽样高,但仍存在两方面的问题:第一,扩展拉丁超立方抽样无法保证序列的差异性,这成为计算效率提高的瓶颈;第二,当输出变量偏离正态分布时,扩展拉丁超立方抽样方法缺乏性能良好的收敛判据。针对以上两个问题,采用基于Sobol序列的扩展准蒙特卡洛方法进行随机潮流计算,并提出基于非参数核密度估计方法的收敛判据。对IEEE 30和IEEE 118节点系统的仿真结果表明,所述方法比扩展拉丁超立方抽样方法更加方便、准确,同时效率更高、收敛更快;而基于非参数核密度估计的收敛判据直观、适应性强,对变量的概率分布没有附加条件,可准确指导扩展随机潮流的收敛。

核密度计 篇6

巴塞尔新资本协议明确将影响某金融资产风险的主要因素界定为三个：一是违约概率PD;二是违约回收率(RR)或违约损失率(LGD);三是违约风险暴露(EAD)。对于违约概率的研究已较为成熟，我们熟悉的三大信用风险模型，结构模型，简约模型，VaR模型均是围绕此因素展开，而在这些信用风险模型中，对于另一项重要因素LGD的处理却相对简单，或将其定为一个与PD不相关的固定值或随机数，或用某种简化的函数形态加以描述。然而随着越来越多的研究关注到LGD，这些假设的不足开始显现出来。诸多研究结果表明[1]，回收率通常呈现双峰分布，大量数据集中在较低和较高的两端，因此单纯的用一个固定数值，如均值，或一个单峰的beta分布来刻画就可能造成偏差。这样，一类不需事先对变量分布形态做出假设的非参数统计方法开始受到广大学者的关注。

非参数估计是根据样本资料对总体的某种性质或密度形态做出估计，它不要求总体分布已知或对总体分布做任何限制性假设，仅以数据点作为概率密度估计的依据。非参数统计问题中对总体分布的假定要求的条件很宽，因而针对这种问题而构造的非参数统计方法不致因为对总体分布的假定不当而导致重大错误，所以它往往有较好的稳健性。非参数方法有诸多种类，如直方图、核函数、K-临近方法等，其中核函数因能最终给出概率密度函数的显式形式而被广泛研究和应用。但是，正基于核估计的原理，它在大样本处的估计效果较好，在样本点稀疏处就可能产生较大误差，这类点多产生于有界区间的边界，这就出现了非参数估计的边界偏差问题(boundary bias problem)。针对这一问题，研究学者们进行了深入研究探索，其中以非对称核法、镜像法、边界核法为主要对策。

二、非对称bet a核

非对称核顾名思义就是指使用的核函数不以某个中间值对称，最先使用此类核的学者要数Chen (1999)，他创造了一个非对称beta核来估计一些随意抽取的点的概率密度[2]。他认为这类核可以改善传统对称核在应付有界区间上密度估计的边界偏差问题。

（一）定义

Beta核定义如下：X1, X2…Xn为取自某分布函数F上的样本点，F具有一未知密度函数f, f分布在闭区间[0, 1]上且二阶连续可导。对于x[0, 1]，f的beta核估计密度函数为：

其中非对称核K(·)是一个beta密度函数：

其中。式中参数b为光滑参数，也叫做窗宽，b需满足如下两个条件：

（二）性质

要了解核函数的性质和评价其估计的优劣，我们需从两个方面入手，其一是在估计点处的性质，叫做局部性质，由核函数f赞(x)在该点处的偏差(bias)和方差(variance)刻画;其二是在整个估计区间上的性质，叫做整体性质，由观察积分均方误差(MISE)得出。

1. 由定义知，的偏差为

Ch en (1999)给出了beta核偏差的推导结果并给出了证明[3]，这里我们引用其结果：

2. 由定义知，f赞(x)的方差为

方差化简结果的推导较为复杂，Chen (1999)也给出了一个推导结果并予以证明[3]，我们在此直接引用其结果：

其中Mb是一个足够大的常数，c是一个常数。对于x的定义域我们给予这样的解释：对于x/b>Mb， (1-x)/b>Mb，那么满足的x∈(b·Mb, 1-b·Mb)的x为在区间内部，对于x/b=c≤Mb， (1-x)/b>Mb或x/b>Mb， (1-x)/b=c≤Mb，那么满足和的x∈[0, Mb·b]和x∈[1-Mb·b, 1]的x处于区间边缘。内部和边缘的x的密度估计具有的方差不同也体现了我们上文提到的边界问题。

由(1)式可以看出，的偏差在整个[0, 1]区间内的阶数保持b的一阶不变，这说明以偏差衡量，beta核不产生边界偏差问题。但是方差在整个区间内有所不同，在内部是n-1b-1/2的同阶，在左右边界处是n-1b-1的同阶，也就是说，随着n的增大，方差的收敛速度在边界处慢于内部，从方差角度讲是存在一定边界问题的[4]。

3.积分均方误差 (MISE) 定义为

其中，[0, c1]为左边界，[c2, 1]为右边界，是方差计算中产生边界问题的区域;[c1, c2]为内部区域。Ch en (1999)认为[2]，我们可以找到一个合适的ε，，使c1=b1-ε，c2=1-c1。上述结果可以进一步利用积分中值定理简化，以上式等号右边第二项为例：

其中θ是函数中的某一点。这样，上式就为n-1b-ε的同阶无穷小。同理可证等号右边第三项第四项。最终可得MISE表达式：

最小化(3)式对b的一阶导数可以求得最优窗宽hoptimal，带回MISE中可求得MISE为n-4/5的同阶无穷小，记为MISE～O (n-4/5)，也就是说，随着n的增大，MISE以n-4/5的速度收敛于0。

三、数值模拟

为了测试非参数法在测量未知形状样本点密度函数上的优劣，我们使用蒙特卡洛抽样方法对比三种拟合方法在事先给定的四种密度函数上的误差大小。三种拟合方法分别沿袭学术界对于回收率的研究发展脉络，选择均值、beta函数拟合、beta核函数拟合三种方法;四种密度函数选择可以代表单峰、倾斜、双峰、平坦的四种函数形式[5];误差大小用积分均方误差(MISE)来表示。

代表单峰、倾斜、双峰、平坦的四种函数形式我们分别选择：

第一, 对数正态分布:, x∈ (0, +∞) , 取参数均值和方差分别为μ=0.8, σ=0.5。

第二，指数分布：f (x)=λe-λx, x∈(0，+∞)，λ=6取参数。

第三，两个正态函数相加：

x∈ (-∞, +∞) , 取四个参数分别为μ1=1, μ2=-1, σ1=σ2=0.2。

第四，均匀分布：f (x)=1, x∈[0, 1]。

为使我们的模拟更加符合我们所要研究的LGD的实际情况，即LGD分布在[0, 1]区间上的事实，我们需将以上(0，+∞)，(0，+∞)，(-∞，+∞)上的x值变换到[0, 1](或(0, 1))上。为使变换不改变原函数的特征，我们使用一个单调变换函数y=g (x)。对于[0，+∞)(或(0, 1))，选用，对于(-∞，+∞)，选用

在抽样之前我们需要说明的一点是, 使用核方法进行估计一个不可忽视的因素是对于窗宽的选择, 窗宽过大导致估计结果过于平滑, 掩盖重要特征形态, 窗宽过小导致结果不够光滑, 放大各个估计金融Finance NO.07, 2012 (CumulativetyNO.486) 细节淹没整体趋势, 因此一个估计的好坏很大程度上取决于窗宽是否合理。本文依据“拇指法则”将窗宽选为, 其中是样本标准差, n是样本数[6], 视情况再做调整。

抽样的具体做法为：

第一，首先我们从真实函数形式中抽取1000个点;

第二，用均值、beta函数、beta核分别拟合这1000个点;

第三，重复抽样100次，对三个拟合结果取平均;

第四，画图，并计算三种方法的MISE。

模拟的结果分别展示在图1和表1中。

图1中实线代表真实密度函数, 虚线代表beta分布模拟结果, 点状线代表beta核模拟结果。由图1我们看出，对于单峰的违约数据，beta分布与beta核的拟合结果相差无几，这主要是因为两种方法中包含的beta分布本身就可根据参数的调整做出不同的单峰情形。MISE值中核函数略好于beta分布。偏斜分布中两者更加几乎无异，MISE也几乎相同。然而在双峰的拟合中，结果却差异明显，虽然beta分布的形状可根据样本数值调整，但单个beta分布却永远无法创造出两个峰的图形，而beta核却不受形状的局限，可以最大化的模仿原函数。两者的MISE结果自然差别很大。在均匀分布的测试中，单纯从MISE结果来看，两者几乎无差，但是从图形上看，在内部点，beta核较为平坦，但是边界处有明显而且较大的偏差，beta分布在边界处有一个急剧的跌落，这与beta分布本身的形状有关。最后，我们给出了每一种函数的均值，显然若以均值方差等简单统计量来描绘某一分布，是无法体现这种分布的重要特征。可以注意到MISE统计中，beta分布有表现良好的时候，双峰及以上情形除外，但不论对于何种情形，beta核函数均表现良好。

考虑到实际研究中，回收率为0和1的两点是非常重要的两种情形，并有诸多学者和机构研究表明，在这两点附近的数据是非常多的，因此我们考虑总体而言，在我们无法准确了解我们所测试的回收率是何种分布形态的前提下，beta核估计是优于直接用beta分布拟合的，这种优势在真实形状特别复杂，尤其是多峰状态下更加明显。

四、结论

本文的研究表明, 在理论上, 非对称beta核函数在整个区间保持偏差阶数不变, 方差阶数在边界有所增大, 从偏差角度来说克服2012年第7期中旬刊 (总第486期) 时代Times了边界问题, 从方法角度来看依然存在一定问题。整体的MISE以n-4/5速度收敛, 收敛速度较快。实际数值模拟上, 在较常规的图形, 如单峰、偏斜情况下, beta核与beta分布差别不大, 边界处略优, 都明显优于均值, 对于较复杂图形, 如双峰或多峰, beta核明显优于有固定函数形式的beta分布, 同时更优于均值。只在均匀分布的前提下, 这两者的效果较均值处于劣势。考虑到近似均匀分布的违约回收率在实际中是很难成立的, 我们认为总体而言, 以函数形式去拟合要明显优于以均值、方差这种简单的统计指标去刻画。这充分体现了对于未知分布形态的数据研究中, 非参数估计方法的优势。考虑到实际数据的参差性, 我们认为在今后的实际研究中, 可将这种方法用于实际回收率的建模上。

参考文献

[1]Edward I.Altman.Default R ecovery R ates and LGD in Credit R isk Modeling and practice[M].London:Cambridge University Press, 2008.

[2]Song Xi Chen.Beta kernel estimators fordensity functions[J].Computational Statistics&Data Analysis, 1999, 31:131-145.

[3]Song Xi Chen.Beta kernel smoothers for regression curves[J].Statistica Sinica, 2000, 10:73-91.

[4][45]Shunpu Zhang, R ohana J.Karunamuni.Boundary perfor-mance of the beta kernel estimators[J].Journal of Nonparametric Statistics, 2010, 2 (1) :81-104.

[5]T.Bouezmarni, O.Scaillet.Consistency of Asymmetric Kernel DensityEstimators and Smoothed Histograms withApplication to Income Data[J].Econometric Theory, 2005, 21 (2) :390-412.

核密度计 篇7

近年来, Copula方法在金融市场相关性的研究越来越广泛。由于金融市场是动态发展的, 时变相关能够更好地描述金融市场间联系的动态结构变化。但在实际应用中, 人们往往更关注的是尾部相关性, 金融资产收益通常表现出尾部相关性, 而尾部相关又往往并不对称。韦艳华, 张世英 (2005) [1]构建了具有尾部变结构特性的R S-Copula-GAR CH模型, 对中国股票市场的研究表明, R S-Copula-GAR CH模型对市场间相关结构, 特别是对尾部相关结构的刻画能力要强于相应的静态Copula模型。本文运用非参数核密度估计技术来刻画沪深股市的对数日收益率, Copula函数采用Patton (2006) 中提出的Symmetric-Joe-Clayton (SJC) Copula函数, 用一个ARMA (1, 10) 过程描述尾部相关的时变性, 结合沪深股市的分析与传统基于Garch模型的估计方法进行了比较, 结果表明基于非参数核密度估计时变相关Copula方法能更好的估计尾部相关性。

二、Copula函数和边缘分布

(一) Copula函数定义

Copula理论最早由Sk lar (1959) [5]提出, 通过将一个联合分布分解成各个变量的边际分布和一个Copula函数, 其中的Copula函数描述变量之间的相关结构。即Copula函数实际上是一类将变量联合累积分布函数同变量边缘累积分布函数连接起来的函数, 因此也有人称其为“连接函数”。常用的Copula函数有正态分布Copula函数、t分布Copula函数、阿基米德Copula函数等 (Gumbel、Clay ton、Fr ank) 。所以要构建Copula模型核心在于首先确定边缘分布, 再选取一个适当的Copula函数, 以便能很好地描述出随机变量之间的相依结构。

(二) 时变非对称Copula函数

GJC (Gener alized Joe-Clay ton) copula表达式为:

本文使用了Patton提出的GJC (广义Joe-Clayton) copula。其表达式为:

三、实证研究

此部分主要讨论非参数核密度估计时变相关Copula在股市相依性上的应用与比较分析。选取上证指数和深成指为研究对象, 本文选择上证综合指数和深成指自1996年12月16日到2012年2月3日间的共3659个日收盘价格数据并计算对数形式的收益率, 对收益率序列进行建模分析。

本文选取正态核函数, 由经验法则选取的带宽分别为0.0034和0.0038, 进行积分变换后可得序列u与v, 其散点图可以看出沪深股市指数有很强的正相关关系。

采用GJC-Copula函数, 其中需要预先设定的是虚拟变量的取值, 本文选取的断点是2005年, 在1996年至2004年取值为0, 2005年至2011年取值为1, 即研究2005至2011年与之前几年我国股市与其它股市的相关性。根据核密度估计变换后的边缘分布序列去估计Copula函数参数结果如下所示:

上尾和下尾相关系数如下图所示:

间断点前后上尾和下尾相关系数均值如下:

从中可以看出, 从总体上来看, 沪深股市在下跌时的相关性要强于上升时的相关性, 即恐慌期强于利好期, 这可能是由于人们在熊市时的反应较为一致, 股市存在很强的“追涨杀跌”效应造成的。同时可以看到上尾和下尾相关系数在股权分置改革后有一定程度的下降, 反映股权分置改革后逐渐有人开始并不一味追逐上涨行情, 并且在股市下跌行情中寻找套利机会, 这可以从下尾相关系数下降强于上尾相关系数看到。

四、结论

本文详细介绍了Copula函数理论和估计方法, 应用非参数核密度估计刻画序列概率分布, 并引入时变相关非对称Copula函数 (GJC-Copula) 研究了沪深股市的相关性, 结果表明股市处在熊市时相关系数更大, 在加入虚拟变量后, 发现股权分置改革后相关系数有所变化, 结果表明05年股权分置改革后沪深股市上尾和下尾相关系数均有一定程度下降, 下尾相关系数下降强于上尾相关系数。Copula函数可以动态的去测量多个时间序列的相依性, 今后应将其更多的应用于汇率、债券等金融序列动态相关研究中。

摘要：近年来, Copula方法在金融市场相关性的研究越来越广泛。由于金融时间序列具有时变非对称相关的特性, 本文引入时变相关非对称Copula函数 (GJC-Copula) 应用非参数核密度估计刻画序列概率分布, 再采用极大似然估计方法估计尾部相关参数, 研究了沪深股市之间的相关性, 结果表明股市在下跌时相关系数更大, 并通过加入虚拟变量方式分析了股权分置改革后相关系数的变化, 结果表明05年股权分置改革后沪深股市上尾和下尾相关系数均有一定程度下降。

关键词：GJC-Copula,Monte Carlo模拟,非参数核密度估计

参考文献

[1]韦艳华, 张世英.金融市场非对称尾部相关结构的研究.管理学报, 2005 (5) :第601-605页.

[2]韦艳华, 张世英.多元Copula-GARCH模型及其在金融风险分析上的应用.数理统计与管理, 2007 (3) :第432-439页.

[3]任仙玲, 张世英.基于核估计及多元阿基米德Copula的投资组合风险分析.管理科学, 2007 (5) :第92-97页.

[4]任仙玲, 张世英.基于非参数核密度估计的Copula函数选择原理.系统工程学报, 2010 (1) :第36-42页.