变参数回归模型

变参数回归模型（精选8篇）

变参数回归模型 第1篇

近年来,并网型风电发展迅猛,风电的波动性已给电网调度带来严峻的挑战。风电功率预测是解决风电波动、实现风电与电力系统传统电源联合优化运行的关键技术之一。

风电功率预测按照预测的时间尺度划分一般分为超短期、短期和中长期预测[1]。超短期预测一般指6 h以内的预测,预测结果用于电力系统的在线优化运行,常采用基于历史风电功率数据的时间序列分析方法进行预测,例如自回归滑动平均(ARMA)模型[2,3]、Kalman滤波[4]等。短期预测一般指对未来6 h～48 h风电功率输出的预测,预测结果是电网安排日发电计划或进行电力市场交易的基础。中长期预测一般指未来几天的预测,预测结果主要用于安排风电机组的检修计划等。在实际应用中,短期预测和超短期预测应用较多。本文的研究对象为短期预测。

由于天气状况在未来6 h～48 h内一般有较大的变化,因此,短期预测主要依赖于数值天气预报(numeric weather prediction,NWP),通过建立NWP的气象信息与风电功率输出之间的关系模型,将预测时段内的气象信息转换为风电功率输出。按建模方法的不同,短期风电功率预测可进一步分为物理方法和统计方法[5]。本文研究方法属于统计方法的范畴。

经过多年的积累,欧洲和美国已经有多款商业化的风电功率预测软件[6],如丹麦的WPPT和Prediktor、西班牙的SIPREOLICO等。由于中国的气候条件与欧美相比差异较大,因此有必要研究适合中国风电场的风电功率预测方法。近几年,中国的风电功率预测研究也在逐步发展,但受气象服务条件的影响,预测方法大多基于历史数据和时间序列方法[7,8,9],对超短期预测较为有效,但对短期(如日前24 h)风电功率的预测效果往往较差。随着风电的大规模接入,为电网安排发电计划服务的短期风电功率预测亟需展开。中国电力科学研究院开发的基于NWP的短期风电功率预测软件[10,11],采用的预测方法为反向传播(BP)神经网络,是一种在风电功率预测中应用广泛的典型方法。但是,神经网络方法对模型训练的时间较长,并且需要不断调试合适的隐含层神经元个数、合适的隐含层输出函数及合适的输出层输出函数等,才能得到收敛性较好的神经网络。非参数回归方法也是模型估计的典型方法之一,在国外已有采用基于统计模型的风电功率预测方法的范例[12]。非参数统计模型只需调整合适的窗宽即可应用模型进行预测,实用性比神经网络模型更佳。

本文以内蒙古某风电场为例,研究将非参数回归方法应用于国内短期风电功率预测的有效性。内蒙古气象局引进了美国国家大气研究中心(NCAR)和美国宾州大学(PSU)开发研制的第5代中尺度模式MM5(Mesoscale Model 5),直接提供风机轮毂高度的NWP信息,如风速、风向等。本文采用内蒙古气象局提供的NWP数据,建立NWP与风电功率输出之间的转换模型,得到风电功率的点预测值;基于经验分布模型和非参数回归方法,建立风电功率预测误差的概率分布函数,进而得到风电功率预测值的概率区间,辅助电网运行决策。

1 非参数估计的基本原理

统计与计量的前沿研究领域是半参数与非参数方法[13,14]。相对于参数估计,非参数估计方法并不假定函数的形式已知,也不设置参数,函数在每一点的值都由数据决定,从而避免模型分布形式选择不当带来的误差。非参数估计方法在天气预报领域应用较为广泛[15,16,17]。使用非参数估计方法进行预报时,不需要建立预报方程,而是直接根据训练数据(历史样本)建立非参数估计模型,利用训练数据中蕴含的输入输出关系进行预报。核函数估计是非参数回归模型中的基本方法之一,其主要思想是在大量历史数据的基础上,应用核函数和一定窗宽范围内的历史数据对某一数据点对应的函数值进行估计或预报。

非参数估计的基本原理[14]如下:对n个给定样本(Y1,X1),(Y2,X2),…,(Yn,Xn),其中,Y为被解释变量,X为解释变量,假定Y1,Y2,…,Yn独立同分布,可建立如下非参数回归模型:

$Y{}_i=g(X{}_i)+\epsilon {}_{i}i=1,2,\cdots ,n(1)$

根据核函数估计的思想,X=x对应的Y值按下式进行估计:

$g(x)=\frac{{\displaystyle \sum _{i=1}^n{\rm K}}(\frac{x-X{}_i}{h})Y{}_i}{{\displaystyle \sum _{i=1}^n{\rm K}}(\frac{x-X{}_i}{h})}(2)$

式中:K(·)为核函数,用以确定样本点Yi(i=1,2,…,n)在估计g(x)中的权重,常用的核函数有均匀核K0(u)=0.5I(|u|≤1),高斯核${\rm K}{}_1(u)=(2\text{π}){}^{-\frac{1}{2}}\exp (-\frac{1}{2}u{}^2)$和抛物线核K2(u)=0.75(1-u2)+等;h为控制局部邻域大小的窗宽,是控制估计精度的重要参数。

最佳的窗宽应既不过小也不过大。窗宽过小会放大随机误差,窗宽过大则会得到过分光滑的曲线,使估计失去意义。常采用交错鉴定法[13]选择最佳窗宽。

2 基于NWP的短期风电功率预测

前已提及,NWP是短期风电功率预测的关键信息,其中,风速预报是影响风电功率输出的关键因素,本文主要考虑将风速预报作为风电功率预测的输入,采用统计方法进行预测。本文的短期风电功率预测分为预报数据校正、风电功率点预测和风电功率概率区间预测3个环节。短期风电功率预测基本流程如图1所示。

2.1 NWP数据校正

针对风电场的气象服务是风电大规模发展和并网运行对气象服务部门提出的新挑战。目前国内外的NWP服务中,针对风电场轮毂高度的风速预报精度较差。因此,使用NWP进行风电功率预测时有必要对预报数据进行校正。由于风速对风电功率预测的影响较大,本文以风速预报作为风电功率预测的输入,采用线性回归方法对NWP提供的风速预报数据进行校正,校正模型如下:

$v{}_{\text{Ν}\text{W}\text{Ρ},t}^{\text{c}}=v{}_{\text{Ν}\text{W}\text{Ρ},t}-e{}_{\text{Ν}\text{W}\text{Ρ},t}(3)$

式中:vNWP,t为校正前的t时刻NWP风速;vcNWP,t为校正后的t时刻NWP风速;eNWP,t=a+bvNWP,t,a和b为参数,可采用最小二乘法,由历史NWP风速预报及其误差样本进行估计,

$\begin{array}{l}a=\overline{e}{}_{\text{Ν}\text{W}\text{Ρ}}-b\overline{v}{}_{\text{Ν}\text{W}\text{Ρ}}\\ b=\frac{{\rm N}{}_{\text{c}}{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{c}}}e}{}_{\text{Ν}\text{W}\text{Ρ},i}v{}_{\text{Ν}\text{W}\text{Ρ},i}-{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{c}}}e}{}_{\text{Ν}\text{W}\text{Ρ},i}{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{c}}}v}{}_{\text{Ν}\text{W}\text{Ρ},i}}{{\rm N}{}_{\text{c}}{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{c}}}v}{}_{\text{Ν}\text{W}\text{Ρ},i}^2-({\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{c}}}v}{}_{\text{Ν}\text{W}\text{Ρ},i}){}^2}\end{array}$

Nc为样本数量;eNWP,i=vNWP,i-vmeas,i,为历史NWP风速预报误差;vmeas,i为风电场实测风速。

2.2 风电功率点预测

风电功率点预测的主要思想是通过风速—风电功率转换模型将预报风速转换为风电功率值。根据非参数估计的思想,建立风速—风电功率转换模型如下:

${\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}=f{}_{\text{p}}(v{}_{\text{Ν}\text{W}\text{Ρ},t+k}^{\text{c}})(4)$

式中:Ppredt|t+k为t时刻得到的t+k时刻的风电功率预测值;k为预测尺度,即待预测的风电功率提前于预测执行时刻t的时间;fp为表征风速与风电功率之间关系的函数;vcNWP,t+k为经校正的t+k时刻的预报风速。

采用样本点(vmeas,i,Pmeas,i) (i=1,2,…,Np)作为式(4)的样本数据,其中,Pmeas,i为风电场实测风电功率数据,Np为样本数量。

2.3 风电功率概率区间预测

概率区间预测[18,19]是描述真实值相对于预测值的不确定性的常用方法之一。由于针对风电场的NWP精度较差,风电功率预测的误差较大,因此有必要对风电功率预测的不确定性进行描述,辅助电网运行决策。

定义风电功率预测误差为某一时间点的风电功率预测值与风电功率实测值之间的偏差,即有:

$e{}_{t|t+k}={\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}-{\rm P}{}_{t+k}^{\text{m}\text{e}\text{a}\text{s}}(5)$

式中:Pmeast+k为t+k时刻的风电功率实测值。

设预测误差et|t+k的概率分布函数为Ft|t+k(ξ),其中,ξ为表示预测误差的随机变量,则真实值Pmeast+k的一个1-α概率预测区间为:

$[{\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}+\widehat{G}{}_{t|t+k}(\alpha {}_1),{\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}+\widehat{G}{}_{t|t+k}(\alpha {}_2)](6)$

式中:α2-α1=1-α,本文取对称概率区间,即α1=α/2,α2=1-α/2;$\widehat{G}{}_{t|t+k}(q)$为概率分布函数Ft|t+k(ξ)的反函数,有Pr$\{\xi \le \widehat{G}{}_{t|t+k}(q)\}=q$。

风电功率预测误差的分布特性是建立风电功率预测区间的基础。关于风电功率预测误差的分布特性已有多篇文献进行探讨。文献[20]指出风速预测误差服从正态分布,但风电功率预测误差不服从正态分布;文献[21,22]认为采用Beta分布描述风电功率预测误差更加合理;Pinson在文献[23]中采用经验分布函数对风电功率预测误差的特性进行描述。经验分布函数[24]不对模型的概率分布函数形式作任何的假设,而是基于历史值计算得到变量的概率分布模型。由于风电功率预测误差受多种因素的影响,尚无一种特定的分布形式可对其进行准确描述,因此,本文采用经验分布模型建立风电功率预测误差的概率分布函数,并在此基础上采用非参数回归技术建立风电功率预测区间。

设Ωt(j)表示时刻t之前的最近n个历史预测误差值的集合,对∀j,即有:

$\Omega {}_t(j)=\{e{}_{t-i|t-i+k}|i=1,2,\cdots ,n;k=1,2,\cdots ,{\rm K}\}(7)$

式中:K为n个历史预测误差中的最大预测尺度。

集合Ωt(j)的经验分布函数$\widehat{F}{}_t(\xi )$定义如下:设集合Ωt(j)的各元素的概率为1/n,则

$\widehat{F}{}_t(\xi )=\frac{1}{n}\text{Ν}\text{u}\text{m}\{e{}_i\le \xi |e{}_i\in \Omega {}_t(j),i=1,2,\cdots ,n\}(8)$

式中:函数Num用于求取集合Ωt(j)中满足给定条件的元素的数量。

设$\widehat{F}{}_t(\xi )$的反函数为$\widehat{G}{}_t(q)$,按式(6),t+k时刻的风电功率Pt+k对应的一个1-α概率预测区间为:

$[{\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}+\widehat{G}{}_t(\frac{\alpha }{2}),{\rm P}{}_{t|t+k}^{\text{p}\text{r}\text{e}\text{d}}+\widehat{G}{}_t(1-\frac{\alpha }{2})〗(9)$

采用非参数估计方法分别对$\widehat{G}{}_t(q)$进行模型估计,建立估计模型如下:

$\xi ={\rm H}(q)(10)$

式中:H为$\widehat{G}{}_t(q)$的非参数估计模型。

3 算例分析

以内蒙古某风电场(以下简称风电场A)为例,采用该风电场的NWP数据、实测风速和功率数据验证非参数回归方法在风电功率预测中的有效性。以2010年1月1日—1月27日的数据作为模型训练数据,2010年1月28日—1月31日的数据作为模型测试数据。风电场容量PA=100 MW。

分别采用平均绝对误差(MAE)和均方根误差(RMSE)对预测效果进行评价,计算公式如下:

$\begin{array}{l}V{}_{\text{Μ}\text{A}\text{E}}=\frac{1}{n}{\displaystyle \sum _{k=1}^{{\rm K}}|}\widehat{x}-x|(11)\\ V{}_{\text{R}\text{Μ}\text{S}\text{E}}=\sqrt{\frac{1}{n}{\displaystyle \sum _{k=1}^{{\rm K}}(}\widehat{x}-x){}^2}(12)\end{array}$

式中:x为预测量的真实值;$\widehat{x}$为预测量的预测值。

3.1 NWP数据校正结果

采用式(3)对NWP风速预报数据进行校正,图2为2010年1月28日—1月31日校正前后的NWP风速预报数据与风电场实测风速数据的对比曲线。

校正前NWP风速预报数据的MAE和RMSE分别为VMAE=1.98 m/s,VRMSE=2.44 m/s;校正后分别为VMAE=1.85 m/s,VRMSE=2.30 m/s。可见,对NWP风速预报数据进行校正只能在一定程度上降低预报误差,真正提高风速预报精度仍有待于NWP精度的提高。

3.2 风电功率点预测结果

选择抛物线核函数,采用式(4)进行风电功率点预测,分别以如下3种不同的风速作为输入:①风电场实测风速;②未经校正的NWP风速;③经校正的NWP风速。不同输入风速下的风电功率预测曲线如图3所示。

不同输入风速下风电机组功率预测的MAE与RMSE对比如表1所示。表中的MAE与RMSE均以风电场容量的百分比表示。

由表1可见,风电场实测风速输入对应的风电功率预测误差最小,未经校正的NWP风速输入对应的风电功率预测误差最大,即NWP误差是风电功率预测的主要误差来源,决定风电功率预测水平。风速向风电功率转换的模型误差较小。采用统计方法对NWP风速进行校正后,风电功率预测精度有所改善,但改善幅度较小。

3.3 风电功率概率区间预测结果

由于风电功率预测误差与风电功率预测值大小、风速等多种因素相关[20,23],因此,本文以预测功率水平作为预测误差的影响因素,将风电场A的风电功率按预测值等分为13个功率等级,分别建立不同预测功率水平下的风电功率预测误差分布。图4为2010年1月1日—1月27日所有处于功率段的风电功率预测值对应的功率预测误差经验概率密度函数。图4中,直方图为误差样本直方图,曲线为应用非参数回归技术拟合得到的预测误差经验概率密度函数。

分别取α=10%与α=30%,由式(9)得到风电场A在2010年1月28日—1月31日的风电功率90%预测区间和70%预测区间分别如图5和图6所示。

可见,采用经验分布函数和非参数回归方法建立的风电功率预测区间反映了其对实际风电功率的覆盖概率。实际应用中,对于预测区间的覆盖概率,一般取较为适中的值,如80%,避免因预测区间过大而失去参考意义和因预测区间过小而难以包含较大的风电功率预测误差。

4 结语

随着风电的大规模发展,短期风电功率预测成为调度运行的关键环节之一。本文将非参数回归技术用于短期风电功率预测,建立了短期风电功率点预测模型和概率区间预测模型,并采用内蒙古某风电场的NWP数据和风电功率实测数据验证了所建立的非参数回归模型在短期风电功率预测中的可行性和有效性。分析结果表明,除风电功率点预测外,采用非参数回归模型和经验分布函数得到的风电功率预测区间可以描述风电功率预测的不确定性,进一步辅助电网运行决策。此外,NWP精度是影响短期风电功率预测精度的主要因素。采用统计方法对NWP数据进行校正对于改善风电功率预测精度的作用十分有限。提高针对风电场的NWP精度对于提高短期风电功率预测精度有重要意义。

目前,基于本文模型研发的短期风电功率预测系统已在内蒙古部分风电场投入使用。随着风电场运行数据的积累和NWP水平的提高,进一步分析风向、气温、天气类型等气象因素对风电功率预测及预测误差的影响将在后续工作中进行。

变参数回归模型 第2篇

对于纵向数据下半参数回归模型,基于广义估计方程和一般权函数方法构造了模型中参数分量和非参数分量的`估计.在适当的条件下证明了参数估计量具有渐近正态性,并得到了非参数回归函数估计量的最优收敛速度.通过模拟研究说明了所提出的估计量在有限样本下的精确性.

作 者：田萍 薛留根 Tian Ping Xue Liugen 作者单位：田萍,Tian Ping(许昌学院数学系,许昌,461000;北京工业大学应用数理学院,北京,100022)

薛留根,Xue Liugen(北京工业大学应用数理学院,北京,100022)

变参数回归模型 第3篇

随着以GPS(Global Position System)为代表的GNSS(Global Navigation Satellite System)技术在大地测量领域的广泛应用,其全天候、高精度、测站间无需通视及高效低成本等特点,使大地测量领域发生了巨大变革。但是,由于GPS高程测量的基准(WGS-84参考椭球面)与传统地区及国家高程测量的基准(大地水准面或似大地水准面)之间的不同,一般我们需对GPS点进行水准联测,这一直制约着GPS三维定位应用的进一步推广。

因此,GPS高程系统(大地高系统)与传统地区及国家高程系统(正高或正常高系统)之间的转换方法一直是GPS数据处理的研究热点之一,即如何精密确定大地水准面高或高程异常(由于我国采用正常高系统,下文主要分析高程异常)的大小。关于不同高程系统的基本原理及其相互关系已广为熟知,本文不再赘述,具体可参考文献[1]。

目前为止,确定GPS高程异常的方法主要包括重力法和几何法。前者主要是综合局部或全球高分辨率、高精度的重力数据、重力场模型、数字高程模型及GPS水准数据,通过Stokes/Molodensky、最小二乘配置等理论确定该地区的似大地水准面模型[2]。而后者则是利用一定数量的GPS水准控制点建立数值逼近模型,以此插值获得未知点高程异常,根据模型类型不同,传统几何法[3]可分三类,函数逼近模型、统计逼近模型和综合逼近模型。

此外,还有神经网络法[4]、支持向量机法,Fourier分析法、小波分析法,有限元法[5]、球冠谐级数法[6],以及结合现有局部或全球重力场模型或地形改正的“移去—拟合—恢复”[7,8]等特殊方法。

其中,重力法具有高精度、高分辨率特点,但其建设成本高,已覆盖的范围有限,且使用保密性强,对于一般单位和工程需要而言,该方法并不现实。上述特殊方法介于重力法和几何法之间,但计算普遍较繁琐,且要求使用者具有测绘以外的专业知识,在应用推广上有一定的困难。几何法虽然没有考虑地球物理特性,但在局部范围的工程应用中,能以较低的成本获得具有较高精度和可靠性的结果,具有较广泛的应用价值。

但是,传统几何法中函数或统计逼近模型受到模型单一性的限制,仅适用于特定地区或范围;综合逼近模型相对单一性质模型而言,适用性和精度都有一定的改进,但尚无法弥补模型本身物理特性的不足,在地形复杂地区仍然效果不佳;而分区或移动拟合逼近由于理论依据不够严密,也存在一定的局限性[9]。

本文将主要利用半参数回归模型及其虚拟观测值解法,建立一种局部GPS高程异常的逼近和分析方法,并通过与传统几何方法的实例对比分析,重点突出该方法在逼近效果和物理解释方面的优势。

1 半参数回归模型

20世纪80年代,数理统计界提出了一种重要的统计模型———半参数回归模型,该模型的特点是既有参数分量又含有非参数分量[10]。参数分量部分可以用来描述函数关系明确的那一部分,如测量问题中的确定性参数;而非参数分量部分可以用来描述函数关系或规律不明确的那一部分。该模型可以看作是参数回归和非参数回归模型的综合,更能充分利用数据中所提供的信息,也具有更强的实际问题解释能力。

一般情况下,半参数回归模型可表示为:

式(1)中,L表示观测向量,B为系数矩阵,X为参数向量,S表示变化规律复杂的非参数向量,Δ表示偶然误差,σ02P-1表示为观测误差的方差。

1.1 补偿最小二乘解法

显然,(1)式是一个不适定模型,在最小二乘估计准则下不具有惟一解。为此,人们引入了补偿最小二乘估计准则,该准则函数既考虑估值和观测值之间的残差拟合,又顾及非参数函数的光滑性,即:

其中,J(S)为在残差拟合基础上对非参数部分的光滑性限制,其矩阵形式可表示为:

其中,矩阵R称为正则化矩阵,是一个适当给定的正定或半正定方阵,用来描述非参数函数的光滑性,α称为正则化参数,来调节残差拟合和非参数光滑效果之间的平衡。

此时,半参数回归模型在补偿最小二乘估计准则下的误差模型为:

其中,该估计方法的关键是如何合理确定正则化矩阵R和正则化参数α。目前为止,确定正则化矩阵R的思路很多,应用最广泛的主要有[11~13]:

(1)若假设非参数分量满足平稳随机信号性质时,可采用先验方差-协方差、自协方差或经验协方差函数来确定;

(2)若假设非参数分量为非随机参数时,一般采用光滑性较好的自然样条函数确定;

(3)若假设观测值为平稳时间序列,则可利用其相邻时刻的非参数分量变化不大的先验假设来确定正则化矩阵。

正则化参数α的选择方法也较多,较实用的有[11~13],交叉核实(CV)法或广义交叉核实(GCV)法、L曲线法、均方误差法等。

1.2 虚拟观测值解法

由于补偿最小二乘解法在合理确定正则化参数和正则化矩阵方面的不足,本文采用基于Helmer方差分量估计的半参数虚拟观测值解法[14],理论依据充足,且更接近传统平差方法。

当实际观测值是一个平稳时间序列时,可假设相邻时刻的非参数分量值变化不大,并将该先验信息作为一类虚拟观测,即:

则其向量形式为

联合(1)式建立虚拟观测和实际观测联合误差模型:

将上式按广义最小二乘准则估计:

最后,可得参数和非参数的估值为:

一般情况下,我们假设虚拟观测与实际观测互相独立,但他们的初始单位权方差不一致。因此,需采用验后方差分量估计方法来解决,一般采用较简单有效的Helmert方差分量估计,首先,利用初始单位权按(9)式进行平差计算;其次,利用(9)式结果按(7)式获得观测残差,并计算出两类观测的验后单位权方差,然后利用验后单位权方差估值重新对两类观测进行定权。按上述重复迭代运算,直至两类观测的单位权方差之比接近1,具体见文献[14]。

2 实例应用与分析

在本节中,我们将给出GPS水准测量实例,来验证上述方法的可行性和有效性。实验数据来源于我国南方某市的GPS控制网,选取了其中联测了国家四等水准的143个GPS水准点,来进行GPS高程异常拟合逼近试验,总面积达(160×120)km2

首先,为了避免因坐标差异引起计算过程不稳定的问题,对该测区内的所有点进行了中心化处理,即使测区坐标原点与其重心重合,如图1所示。

其次,分别设计三种不同方案对测区内所有GPS水准点的高程异常进行数值逼近分析。

方案a:函数逼近模型(二次曲面拟合);

方案b:半参数回归模型(补偿最小二乘解),利用文献[13]的方法确定正则化矩阵和参数;

方案c:半参数回归模型(虚拟观测值解),本文的新方法。

最后,分别给出各种方案的试验结果,主要包括拟合后残差的大小、分布及其精度分析,详见表1、图2、图3及图4。

3 结论

分析上述计算结果可以看出:

(1)总体上看,由于该测区范围较大,地形存在一定的起伏,函数模型逼近(二次曲面拟合法)的结果显然无法满足实际应用的精度要求,而半参数回归模型能明显改善拟合精度;

(2)相对二次曲面拟合法而言,半参数回归模型的补偿最小二乘解法将拟合精度提高到了厘米级,且残差明显变小,残差分布趋于合理,已经适用于一般的GPS高程拟合;

(3)虚拟观测值解法则将精度提高到了毫米级,残差变化甚小,其分布直方图更加陡峭。其拟合精度都优于半参数补偿最小二乘解,且该方法挖掘先验信息作为虚拟观测参与建模,同时利用Helmert方差分量估计理论合理确定两类观测的权阵,与传统平差理论和方法紧密结合。这充分证明了半参数回归模型虚拟观测值解法在GPS高程异常拟合逼近中的可行性和有效性。

(4)从图4中可以看出二次曲面拟合的残差曲线呈一定的系统性,而半参数回归模型中的非参数能有效地估计出这一系统效应,从而明显地提高了拟合精度,同时也说明了半参数回归模型能有效地分离系统误差,该模型具有较强的实际解释能力。

摘要：本文在分析现有GPS高程异常确定方法的基础上,提出了基于半参数回归模型的GPS高程异常逼近和分析的新方法。本文还通过实例对比和分析,证明了新方法在逼近效果和物理解释方面的可行性和有效性。

变参数回归模型 第4篇

滚动轴承是旋转机械中应用最广泛的零部件之一,轴承工作状况的好坏决定着机器能否正常工作。统计显示,旋转机械设备的功能失效有30%是由轴承故障引起的,因此,对滚动轴承进行故障诊断具有十分重要的意义。

轴承运行时的振动信号是典型的非线性非平稳时间序列,很难用一个完全确定的数学函数来表达。因而对轴承进行故障诊断常通过提取振动信号的特征参数并建立其与运行状态之间的关系来实现。时间序列分析(time series analysis)通过将观测数据拟合为一个参数模型,实现对系统动态特征与内在结构关系的近似描述[1]。常用的时间序列参数模型,如自回归(autoregressive,AR)模型、滑动平均(moving average,MA)模型和自回归滑动平均(autoregressive moving average,ARMA)模型等,都是在假定数据序列为平稳的条件下建立的,在实际应用中存在一定的局限性。时变参数模型因其具有较高的非平稳信号时频分布分辨率而受到普遍关注[2,3,4,5,6]。时变自回归 (time-varying autoregressive,TVAR) 模型是目前应用最多的一种时变参数模型,如文献[5]通过对滚动轴承振动信号的TVAR模型求取的时频谱进行奇异值分解,提取奇异值作为特征参数,文献[6]对转子系统振动信号的TVAR模型,提取其基函数的组合权值作为特征参数,实现了对旋转机械的状态监测和故障诊断。目前,直接提取振动信号TVAR模型的时变参数作为反映机械设备运行状态特征参数的研究还不多见。

通过对轴承振动信号的TVAR模型的时变参数进行大量实验研究分析,发现时变参数能有效利用信号的时频分布信息,较好地表征非平稳信号的动态特征。为此,本文提出一种基于时变自回归参数模型的滚动轴承智能故障诊断方法。首先对滚动轴承的振动信号建立时变自回归参数模型,提取时变参数的均值作为反映轴承运行状态的特征参数,然后采用支持向量机(support vector machines,SVM)分类器实现对滚动轴承的智能故障诊断。

1 TVAR建模原理

对离散时间序列{x(t),t=1,2,…,N}可建立如下TVAR模型:

$x(t)=-{\displaystyle \sum _{i=1}^{p}a}{}_i(t)x(t-i)+\epsilon (t)$ (1)

式中,ai(t)为时变参数;p为模型的阶次;ε(t)为均值为零、方差为σ2的白噪声。

时变参数ai(t)可用一组基函数的线性组合来表示[7],即

$a{}_i(t)={\displaystyle \sum _{j=0}^{m}a}{}_{ij}g{}_j(t)$ (2)

式中,m为基函数维数;ai j为基函数的组合权值;gj(t)为一组基函数。

常用的基函数有时间多项式基函数、傅里叶基函数、离散余弦基函数、勒让德多项式基函数和离散长球序列基函数等。由于轴承的振动信号可近似看作是循环周期信号[8,9],而傅里叶基函数比较适于周期变化情况[10],因此,本文采用傅里叶基函数,其表达形式为

将式(2)代入式(1)可得

$x(t)=-{\displaystyle \sum _{i=1}^p{\displaystyle \sum _{j=0}^{m}a}}{}_{ij}g{}_j(t)x(t-i)+\epsilon (t)$ (4)

令

AT=[a10 … a1m … ap0 … apm] (5)

XTt=[x(t-1)g0(t) … x(t-1)gm(t) …x(t-p)g0(t) … x(t-p)gm(t)] (6)

可将式(4)表述为最小二乘形式:

x(t)=-XTtA+ε(t) (7)

根据最小二乘原理,未知参数A的最佳估计值$\widehat{\mathit{A}}$需使模型的残差ε(t)的平方和达到最小:

$\min [{\displaystyle \sum _{t=p+1}^{{\rm N}}\epsilon }{}^2(t)]=\min \{{\displaystyle \sum _{t=p+1}^{{\rm N}}[x(t)+\mathit{X}{}_t^{\text{Τ}}\mathit{A}]}{}^2\}$ (8)

由此得到A的最小二乘估计$\widehat{\mathit{A}}$:

$\widehat{\mathit{A}}=-({\displaystyle \sum _{t=p+1}^{{\rm N}}\mathit{X}}{}_t\mathit{X}{}_t^{\text{Τ}}){}^{-1}{\displaystyle \sum _{t=p+1}^{{\rm N}}\mathit{X}}{}_{t}x(t)$ (9)

进而可以求得模型残差ε(t)的方差σ2的最小二乘估计值:

$\widehat{\sigma }{}^2=\frac{1}{{\rm N}-p}{\displaystyle \sum _{t=p+1}^{{\rm N}}[x(t)+{\displaystyle \sum _{i=1}^{p}a}{}_t(t)x(t-i)]}{}^2$ (10)

利用式(9)求解模型参数时存在矩阵求逆的问题,当矩阵较大时,求解式(9)需要很大的存储空间和很长的时间。可以利用递推算法来求解式(9)。令$\mathit{{\rm P}}{}_{{\rm N}}=({\displaystyle \sum _{t=p+1}^{{\rm N}}\mathit{X}}{}_t\mathit{X}{}_t^{\text{Τ}}){}^{-1}$,则有如下递推算法[11]:

其中,初始值$\widehat{\mathit{A}}{}_0=0,\mathit{{\rm P}}{}_0=\mu \mathit{{\rm I}}$(选μ远大于1,I为单位矩阵)。

估计出参数a10,…,a1m,…,ap0,…,ap m后,就可根据式(2)求出模型参数ai(t)在各个时刻的值。

在实际应用中,模型的阶次通常是未知的,需要根据观测数据来适当地判定。AIC准则(akaike information criterion)是目前应用比较广泛的定阶方法,通过使平均对数似然函数为最大或Kullback信息量最小来确定AR模型的最佳阶次,其定义如下[12]:

$A{\rm I}C(p)={\rm N}\ln \widehat{\sigma }{}^2+2pm{}^2$ (12)

当p、m增大时,$\widehat{\sigma }$2单调下降,AIC(p)将在某个p=popt、m=mopt处达到最小,则popt、mopt就是最佳模型阶次。

2 基于TVAR模型的智能故障诊断方法

本文提出的滚动轴承智能故障诊断方法流程见图1,首先对滚动轴承运行时的振动信号建立时变自回归参数模型,用模型参数的均值构建特征向量$\mathit{X}{}_i=(\overline{a}{}_1^{(i)}(t)‚\overline{a}{}_2^{(i)}(t)‚\cdots ‚\overline{a}{}_p^{(i)}(t))$来表示第i个样本信号(i=1,2,…,N),然后采用支持向量机实现状态的智能识别。

具体实现过程如下:

(1)分别在轴承正常、内环故障、滚动体故障和外环故障状态下,按一定的采样频率进行采样,获得一定数量的振动信号作为样本(1024个采样点为一个样本)。

(2)对来自传感器的含有大量噪声的原始信号进行低通滤波等预处理,得到待分析信号x(t),利用式(12)确定p、m,然后由式(11)估计出参数A,建立时变自回归参数模型。

(3)构建特征向量$\mathit{X}{}_i=(\overline{a}{}_1^{(i)}(t),\overline{a}{}_2^{(i)}(t),\cdots ,\overline{a}{}_p^{(i)}(t))$,并由Xi构造出训练样本集{(X${}_j^{\text{t}\text{r}\text{a}\text{i}\text{n}}$,yj)}和测试样本集{(X${}_k^{\text{t}\text{e}\text{s}\text{t}}$,yk)},其中,yj,yk∈{1,2,3,4},yj、yk为分类器目标输出值(1表示正常状态,2表示内环故障,3表示滚动体故障,4表示外环故障)。

(4)利用训练样本对SVM分类器进行训练,SVM的建立详见文献[13]。

(5)用训练好的SVM多类分类器对测试样本进行故障模式识别。

3 实验分析

3.1实验Ⅰ

本部分采用的数据全部来自美国 Case Western Reserve大学滚动轴承数据中心[14]。测试平台如图2所示,由1492W(2马力)三相交流电机(左)、转矩传感器(中)、测力计(右)和电子控制装置(未显示)组成。电机轴由测试轴承支撑,通过放电加工(electro-discharge machining,EDM)技术在测试轴承中植入单一局部故障缺陷,故障点直径分别为177.8μm(7mil)、355.6μm(14mil)和533.4μm(21mil),深度为279.4μm,测试轴承为6025-2RS JEM SKF型深沟轴承。

振动信号通过加速度传感器采集,采样频率为48kHz。在不同电机负载/转速工况工作条件(0、746W、1492W和2238W)、4种不同故障类型(正常状态、滚动体故障、内环故障、外环故障)和3种不同故障程度情况下记录振动加速度信号数据,共获得40组数据(正常状态无故障大小,只需采集不同工况下的4组数据),并将每组数据以1024个采样点为一个样本构成一组样本集。按相同负载/转速、相同故障程度分类,可将40组样本集分成12个数据集,每个数据集由在相同负载/转速和相同故障点直径条件下的4种不同故障类型样本集组成,如表1所示,表中,N表示正常状态,I表示内环故障,B表示滚动体故障,O表示外环故障(下同)。以相同负载/转速、相同故障类型分类(不包含正常情况),亦可将36组样本集分成12个数据集,每个数据集由在相同负载/转速和相同故障类型条件下的3种不同故障点直径的样本集组成,如表2所示,表中,FD7表示故障直径为177.8μm,FD14表示故障直径为355.6μm,FD21表示故障直径为533.4μm(下同)。按第2节所述方法对各数据集进行仿真实验研究,在此以数据集D007_0为例,对具体过程加以详细说明。

选取4种故障类型的样本各一个,进行模型阶次(p×m)的确定。得到正常状态的最佳模型阶次为19×4,内环故障的最佳模型阶次为18×3,滚动体故障的最佳模型阶次为18×3,外环故障的最佳模型阶次为20×3。为便于进行SVM分类,需对数据集中的各样本建立相同阶次的TVAR模型,为此,选取最佳模型阶次为18×3。然后求出数据集中各样本的时变参数,图3所示为数据集D007_0中每类故障各一个样本的部分时变参数情况。可见,不同故障类型的时变参数之间具有较大的区分度,如正常和内环故障的时变参数a1(t)与滚动体和外环故障的时变参数a1(t)之间有较大的差异等。因此,时变自回归模型的参数可以用作反映轴承运行状态的特征参数。

如果直接将时变参数作为表征轴承运行状态的特征向量输入到SVM,则一个样本的输入即为18×1024=18 432维的高维向量,不利于SVM的处理,且影响分类效率。为此,在求出轴承振动信号的TVAR模型参数的基础上,需进一步研究时变参数的特性来获取低维特征参数,用以表征轴承的运行状态。本文在建立轴承振动信号的TVAR模型后,分别对每一个时变参数求平均,这样18个时变参数可获得18个平均值,则特征参数的维数降至18维。在大量实验研究分析的基础上发现,在轴承的不同运行状态下,对应的TVAR模型参数的均值亦有着较大的可分性。以数据集D007_0为例,各样本的部分TVAR模型参数的均值如图4所示,显然,不同故障类型所对应TVAR模型参数的均值之间仍具有较大的区分度。选取表2中的数据集DINN_2进行分析,结果如图5所示,可见,不同故障程度所对应TVAR模型参数的均值之间同样具有较大的区分度。大量实验研究分析表明,时变自回归模型参数的均值可以较好地刻画运行在不同负载条件下不同故障类型、不同故障程度轴承的动态特性,因而可以用作描述轴承振动信号的特征参数。

因此,对轴承的振动信号建立18×3阶TVAR模型,采用模型参数的均值作为反应轴承运行状态的特征参数,即构建特征向量$\mathit{X}{}_i=(\overline{a}{}_1^{(i)}(t),\overline{a}{}_2^{(i)}(t),\cdots ,\overline{a}{}_p^{(i)}(t))$,用以反映轴承的运行状态。由特征向量Xi构造出训练样本集{(X${}_i^{\text{t}\text{r}\text{a}\text{i}\text{n}}$,yi)}和测试样本集{(X${}_j^{\text{t}\text{e}\text{s}\text{t}}$,yj)},利用训练样本对SVM分类器进行训练,并用训练好的SVM多类分类器对测试样本进行故障模式识别。

对表1、表2中的数据集分别进行仿真实验,结果如表3和表4所示,各状态均有较高的分类识别率。仿真实验结果说明,采用本文所提的特征提取方法可以有效地提取出反映轴承运行状态的特征参数,从而实现滚动轴承故障的有效诊断。将本文方法与文献[6]所提方法进行对比实验(均采用SVM作为分类器), 限于篇幅, 此处只列出各数据集整体分类识别率,分别如表3和表4中文献[6]方法一栏所示,可知,采用本文方法所获得的整体分类识别率均高于文献[6]所述方法,结果进一步说明了本文所提方法的有效性。

3.2实验Ⅱ

本实验是在笔者所在实验室的QPZZ-Ⅱ型旋转机械故障平台上进行的,装置如图6所示,轴承型号为N205型滚珠轴承。分别在不同负载、转速下通过安装在轴承座外壳上的加速度传感器采集轴承运行时的振动信号,采样频率为20kHz,模拟的故障分别为外环、内环、滚动体出现细小裂纹的情况。然后采用第2节所述的方法对其进行故障诊断,同样以1024个采样点为一个样本,得到各状态的训练样本和待识别(测试)样本数分别为20和480,诊断结果如表5所示。由表5可知,在不同的故障发生形式下,采用本文所提的方法对各种故障类型进行分析亦具有较高的分类识别率,说明本文所提方法具有一定的实用性。

4 结束语

通过对滚动轴承的振动信号建立时变自回归参数模型,提出了一种用模型的时变参数均值构建表征轴承运行状态特征参数的方法,经支持向量机分类器对所提取特征进行故障诊断与分类,实现了滚动轴承故障的智能诊断。利用不同运行状态下的数据进行仿真实验,结果证明了本文所提方法的有效性。

变参数回归模型 第5篇

在处理测量数据过程中, 通常假设观测值中不含有系统误差, 即只含偶然误差。但是, 这种假定存在着一定程度的缺陷。因此, 在实际问题中我们可以通过建立数据模型来对数据函数近似表达。对于在大地测量中有各种常规的静态模型, 这种假设与实际比较一致。这是因为常规仪器在地面静态测量中通过完整的观测程序和严密的校检方法而把系统误差消除。随着现代测量技术的发展, 经典的参数数据处理方法有时不能满足精度的要求, 尤其是当测量数据中包含不可消除的复杂的系统误差时。因此20世纪80年代发展起来的半参数模型越来越受到关注。目前为止, 对半参数模型的研究, 归纳起来主要的估计方法有:补偿最小二乘法、二阶段估计法、两步估计法和稳健-M估计等。但所用的参数和语言都是纯数学的、相对抽象的, 与具体应用中的实际意义关系不大。例如补偿最小二乘准则中的两个重要量:光滑因子和正则矩阵;偏核光滑估计中的光滑矩阵, 这些量都是一些纯数学含义的量, 没有很直接的测量上的意义。同时这些方法的研究更多的是集中于由于非参数分量的引入, 函数方程如何求解的问题上, 对系统误差仍然有过多的假设。就是通过对半参数回归模型的研究, 把半参数回归模型应用到数据处理的过程中, 得出半参数回归模型可以更好将数据中的系统误差分离出来, 得到比较理想的计算结果。

1 半参数模型及其解算方法

1.1 半参数模型[4]

经典最小二乘平差模型为:[1]

式中, L为n维观测向量, △为n维误差向量, B为列满秩设计矩阵, X为t维参数向量, t为必要观测数。若模型 (1) 中含有系统误差, 则, 为了削弱或消除系统误差对参数估计值的影响, 可以将系统误差与偶然误差分开, 在上述模型中引入非参数分量S, 从而将经典最小二乘平差模型改写为如下形式:

其中△'表示偶然误差, , , 为一个描述模型误差或观测值系统误差的n维未知非随机向量, Si表示第i个观测值中的系统误差 (i=1, 2, …, n) , 通常称之为非参数分量;p为对称正定矩阵, 是观测值L的权。这样在观测方程中既有参数分量又有非参数分量, 因此式 (2) 称为半参数模型。

1.2 半参数模型的解算方法

一般最小二乘配置模型[2]

式中L为观测向量, Y为随机参数向量, X为信号, △为观测噪声, , X1为观测点信号, X2为未测点信号。, 和X的随机特性为:

与 (2.2.1) 对应的误差方程为

按照补偿最小二乘原理估计的准则, 在平差准则中, 引入补偿项, 即

式中P为对称正定方阵, R为正规化矩阵, 为平滑因子;根据半参数的解算方法得到参数的估值如下:

式中。这样当正规化矩阵R正定时, 根据实际问题, 通过选取适当的平滑因子和正规化矩阵R, 就可以得到、的唯一解。

2 正规化矩阵R和平滑因子的选择方法

从上面的半参数模型的计算方法可以看出, 在数据处理中, 要想得到比较理想的计算方法, 就必须选取适当的平滑因子和正规化矩阵R。

2.1 距离法选取R

为了度量非参数分量S, 可选取某种距离法, 即

其中Sij为ti与tj之间距离[4]。如取, 若t表示平面坐标, 则距离:

2.2 时间序列法的特性选取R

把观测值假设成一个时间序列, 若相邻时刻的模型误差si与si+1的差别非常小[3], 则可以用如下模型表达正则化矩阵, 取相邻观测点模型误差之差的平方和, 即

2.3 交叉核实法选取a

交叉核实法其基本思想是:通过选取合适的得到拟合曲线, 以此曲线为标准进行预测, 使所有预测点的均方差误差最小。从交叉核实法的思想知道, 若是一个很小的数值, 则是观测值Li的一个很好的估计。删除ti时刻的观测值Li, 也就是将观测值中相应的Li去掉, 即, 根据给定的平滑因子a, 相应的补偿最小二乘原理为:

上式求得的解为:, 可以将ti时刻的系统误差的估计内插得到观测值的估计, 若a选取合适的话, 则

由式 (10) 和 (11) 可以看出, 平滑因子选取的是否合适对交叉核实法起到很重要的作用, 通常情况下, 如果选取的使式 (10) 达到最小, 就将这个作为平滑因子的估计[3]。从上面也可以看出公式 (10) 和 (11) 的计算量是非常大的, 为了方便计算, 实际应用中有下列公式:

由式 (12) 推算可以得到

3 实例比较

3.1假设观测量, , , , , 偶然误差分布, 观测方程为观测值权P取单位权 (这里采用的方法是用半参数模型的补偿最小二乘法解算) 。

由于观测值是时间的函数, 在这里的R选取采用时间序列法中的选取方法。根据式 (6) , (7) 当a=0.8时, 。 (上面图形中系统误差的比较, 光滑曲线代表系统误差的真值, 折线代表系统误差的估计值。)

数回归模型可以将数据中的系统误差分离出来, 得到的结果也比较接近真实值。但是进行半参数回归模型的计算过程中要考虑正规化矩阵R和平滑因子a的选取。

3.2本实例取自参考文献[1], 设已知四个观测点的重力异常观测值L和它们的坐xiyi (列于表1) , 观测误差 (噪声) 的方差为, 试求和未测点的重力异常估值。

正规化矩阵选用距离法中规定的选取方法, , 利用广义交叉核实法, 选取平滑因子, 计算平滑因子a=1, 由式 (6) 、 (7) 计算得:,

的重力异常估值为:;未测点和的重力异常估值为:;

比较参考文献[1]用最小二乘配置法的计算结果, 可以看出半参数回归模型也可以应用到重力异常计算当中, 并且能有效的分离出数据中的系统误差。

结束语

通过模拟算例和实例可以看出半参数回归模型可以应用到数据处理中, 但是半参数回归模型只有将近30年的发展历程, 它还没有一个完整的模型, 所以在短时间内应用半参数回归模型处理数据的系统误差可以得到理想的结果;随着时间的加长半参数回归模型的参数也要跟着改变。总之半参数回归模型是一个新的发展领域, 很多方面值得我们继续研究。

摘要：在数据处理过程中, 随机因素和确定因素往往是共同存在的, 经典平差往往认为观测值中仅含有偶然误差, 忽略了具有确定性的系统误差。本文的目的就是用半参数回归模型L=BX+S+△将数据处理中的系统误差分离出来, 得出相对比较理想的结果。

关键词：半参数,补偿最小二乘法,数据处理

参考文献

[1]崔希璋, 龄宗铸等, 广义测量平差 (新版) [M].武汉:武汉测绘科技大学出版社, 2001.

[2]潘雄, 半参数模型的估计理论及其应用[D].武汉:武汉大学, 2005

[3]沙跃进, 最小二乘配置法在GPS高程拟合中的应用[J].测绘信息与工程, 2000, 3, 3-5.

[4]胡宏昌, 半参数模型的估计方法及其应用[D].武汉:武汉大学, 2004.

变参数回归模型 第6篇

差分演化 (Differential Evolution, DE) 是一种基于群体差异的演化算法, 该算法是Rainer Storn 和Kenneth Price在1996 年为求解切比雪夫多项式而提出的。差分演化算法在首届IEEE演化计算大赛中表现超群, 已经在数字滤波、化工、阵列天线方向图综合、机械优化设计等领域得到了广泛的应用。

回归分析中的参数估计是指在实际问题中随机变量分布函数的形式己知, 但其中参数未知的情况。如果得到了随机变量的一组样本值后, 希望利用样本值来估计变量分布中的参数值, 这在工程中是一个比较重要的问题。在回归分析中, 最大似然估计法是模型参数估计的基本方法。但在用该方法进行参数估计时, 一般要求解联立的超越方程组, 相当复杂, 用常规迭代算法不易求解, 而且收敛性较差, 甚至有时不能收敛。本文采用DE算法, 以最大似然准则作为适应度函数, 建立回归分析中的参数估算模型。探讨多元线性回归中的参数估计计算。数值仿真分析表明, DE算法可以精确地计算出相关参数。

1 DE算法

差分演化是一种基于实数编码的演化算法, 算法的基本思想及整体构架与遗传算法相类似, 从一代种群到下一代种群都要经过变异、交叉、选择等操作, 也一样有几个至关重要的参数必须事先确定。下面逐一介绍差分演化算法的几个关键性的操作。

1.1 参数的确定

差分演化算法主要涉及以下4个参数:①种群规模大小N;②个体的维数D;③变异因子F;④交叉概率CR。有研究结果表明:群体规模N一般介于5D-10D之间;变异因子F在 (0, 2) 之间取值, 一般F=0.5;交叉概率CR一般在[0, l]之间选择, 一般来说, CR越大, 收敛速度越快, 但易于早熟, 易于陷入局部最优, 算法的稳健性越差, 比较好的选择是CR=0.3。当然这些都只是经验值, 没有严密的理论证明, 对于某些具体的问题也可能取其它的值会得到更好的结果, 需要具体问题具体分析。

1.2 生成初始种群

在D维空间里随机产生满足约束条件的N个染色体, 具体过程如下:

undefined

其中xundefined和xundefined分别是第j个变量的上下界, rand (0, 1) 返回[0, 1]之间的随机数。

1.3 变异操作

一般情况下变异操作有如下两种方式:

Scheme DE1:对于群体中个体xr1 (1≤r1≤N) , 由此产生的新个体x'r1 (1≤r1≤N) 满足下式:

undefined

其中r2, r3∈[1, N], r1≠r2≠r3, F>0为放缩因子。

Scheme DE2:对于群体中个体xr1 (1≤r1≤N) , 由此产生的新个体r'r1 (1≤r1≤N) 满足下式:

undefined

其中r2, r3∈[1, N], r1≠r2≠r3, F>0为放缩因子, xbset为当前产生的最优个体, λ为增加的一个控制变量, 一般取和F相同的值。

1.4 交叉操作

交叉操作是对群体每一个个体xr1 (1≤r1≤N) 以及由此经过变异产生的新个体x'r1 (1≤r1≤N) 之间进行的, xr1和x'r1经过交叉以后产生子代的候选个体v, 经过后面的选择操作, 确定是个体xr1还是个体v保留到下一代。交叉操作主要有两种形式:bin方式和exp方式。

1.5 选择操作

经过交叉和变异以后产生了新的个体v, 根据目标函数值的大小, 从xr1和v中选择一个遗传到下一代。如下式 (以求函数最小值为例) :

undefined

2 多元线性回归模型

线性回归模型在定量分析的实际研究中是最流行的统计分析方法。在许多实际问题中, 某个变量Y往往相关于另外一些变量X1, X2, …, Xp-1, 但是这种相关关系或者由于其机理不甚明确, 或者由于问题的复杂性而不能确切知道, 因此只能说由X1, X2, …, Xp-1的取值部分确定Y的取值。在这些情况下, 可以认为Y的值由两部分构成, 一部分是由X1, X2, …, Xp-1能够决定的部分, 它是X1, X2, …, Xp-1的某个函数, 记为f (X1, X2, …, Xp-1) ;另一部分是众多未加考虑因素 (包括随机因素) 所产生的影响, 被看作是随机误差, 记为ε。于是Y与X1, X2, …, Xp-1的关系可以表示为:

undefined

回归分析即利用Y与X1, X2, …, Xp-1的观测数据, 并在误差项的某些假定下确定f (X1, X2, …, Xp-1) 。利用统计推断方法对所确定的函数的合理性以及由此关系所揭示的Y与X1, X2, …, Xp-1的关系作分析, 进一步应用于预测、控制等问题, 特别是当f (X1, X2, …, Xp-1) 是X1, X2, …, Xp-1的线性函数时, 有:

undefined

此模型称为线性回归模型, 其中β0, β1, …, βp-1是未知常数, 称为回归参数或回归系数;Y称为因变量或响应变量;X1, X2, …, Xp-1称为自变量或回归变量;ε称为随机误差项并假定E (ε) =0。ε是不可观测的随机变量, 而Y与X1, X2, …, Xp-1是可观测的变量。这里只讨论自变量X1, X2, …, Xp-1是非随机变量的情形, 而y与ε有关, 是随机变量, 但它是可观测的。

3 算例分析

3.1 计算模型

多元线性回归模型中的一组参数看作一个个体, 种群中的每一个个体代表模型估计问题中的一个候选解, 于是第i个个体pi表示为:

undefined

另外, 定义适应度函数来评价种群中的每个个体, 根据最大似然估计法, 定义适应度函数如下:

undefined

最优解就是使适应度函数Q (θ) 最小的个体。

3.2 与Weka计算结果的比较

3.2.1 测试1

该算例是Weka3.6.2中自带的CPU性能测试的实例, 多元线性回归示范模型如下:

undefined

其中, 变量X1表示周期 (MYCT, 单位ns) ;变量X2表示最小主存 (MMIN, 单位KB) ;变量X3表示最大主存 (MMAX, 单位KB) ;变量X4表示高速缓存 (CACH, 单位KB) ;变量X5表示最小信道 (CHMIN) ;变量X4表示最大信道 (CHMAX) 。具体数据见表1。

DE算法参数选择如下:个体大小为7;种群大小M=30;迭代次数为800;放缩因子F=0.5;交叉概率CR=0.3。

表2为该算例分别经过DE算法和Weka计算得到的参数估计值的比较。

经过DE算法得到算例1的多元线性回归模型算例的参数估计值为:

undefined

与Weka计算的结果进行对比发现, 使用DE算法得到的参数估计值与Weka软件计算得到的结果一致。

3.2.2 测试2

该算例是一组数据量较小的例子。设已知因变量Y的自变量X1, X2, X3, 共得18组数据, 并已知Y对Xi存在着线性关系, 多元线性回归方程为Y=a0+a1X1+a2X2+a3X3, 求其回归方程。样本数据见表3。

DE算法参数选择如下:个体大小为4;种群大小M=30;迭代次数为800;放缩因子F=0.5;交叉概率CR=0.3。

表4为该算例分别经过DE算法和Weka计算得到的参数估计值的比较。

经过DE算法得到算例2的多元线性回归模型算例的参数估计值为:

undefined

与Weka计算的结果进行对比发现, 使用DE算法得到的参数估计值与Weka软件计算得到的结果一致。

3.3 计算结果分析

使用相对平方根误差 (Root relative squared error) 来评估检验多元线性回归模型的拟合度。计算公式为:

undefined

其中:undefined

p为预测值, a为真实值。

上述两个算例中, 由DE算法得到的相对平方根误差和Weka得到的相对平方根误差的比较如表5所示。

本文通过一些其它数据集的测试, 发现当数据源的数量庞大且有明显的线性关系时, 使用DE算法得到的参数估计值与Weka软件测试得到的结果具有较高的一致性;当数据源的数量比较小且不呈明显线性关系时, 通过DE算法也可以得到参数估计值, 由于Weka本身计算得到的参数估计值就有一定的偏差, 所以不具有比较性。

4 结束语

通过实际算例的仿真实验结果表明, 用DE算法得到的回归结果与Weka统计结果有较好的拟合度, 可以肯定地得知用DE算法来估算回归模型中的参数是可行、可信的。

参考文献

[1]R STORN, K PRICE.Differential evolution-a simple and efficientadaptive scheme for global optimization over continuous spaces[R].Technical report International Computer Science Institute, Berkley, 1995.

[2]ABBASS H A, SARKER R, NEWTON C.PDE:A pareto-frontierdifferential evolution approach for multi-objective optimizationproblems[A].Proceedings of the Congress on Evolutionary Com-putation 2001 (CEC2001) .Volume 2, Piscataway, New Jersey, IEEE Service Center, 2001.

[3]XUE F, SANDERSON A C, GRAVES R J.Pareto-based multi-ob-jective differential evolution[A].Proceedings of the 2003Congresson Evolutionary Computation (CEC2003) .Volume 2, Canberra, Australia, IEEE Press, 2003.

[4]田雨波, 朱人杰, 薛权祥.粒子群优化算法中惯性权重的研究进展[J].计算机工程与应用, 2008 (23) .

[5]王丽, 王晓凯.一种非线性改变惯性权重的粒子群算法[J].计算机工程与应用, 2007 (4) .

[6]王启付, 王战江, 王书亭.一种动态改变惯性权重的粒子群优化算法[J].中国机械工程, 2005 (11) .

变参数回归模型 第7篇

关键词：冷成型钢承重墙,硅酸钙板,数值模拟,热工参数,耐火极限

冷成型钢结构作为木结构的一种替代结构体系是通过自攻自钻螺钉将龙骨立柱、导轨等连接成轻钢骨架,并使用螺钉覆以骨架石膏板、玻镁板、硅酸钙板等建筑板材,以形成楼盖、墙体等承重结构体系。具有建筑材料节约、可回收等一系列优点,是目前符合化解产能过剩、藏钢于民政策理念的理想建筑结构形式。

我国人多地少,引进并将冷成型钢结构体系向多层建筑发展更符合国情,国外已有较成熟的应用范例。如:美国有用薄板轻钢房屋建成的6层公寓,俄亥俄州还建造了6层、9层的轻钢住宅;加拿大建有6层的轻钢住宅,也建有将轻钢墙柱结构体系使用到8层的旅馆。

我国强调建筑的被动防火性能,即通过限定建筑构件的耐火极限,实现良好的防火分隔,把住宅防火问题引入单个住户单元。因此,作为主要承受竖向和侧向荷载的承重墙构件,将面临更为严格的耐火要求。而墙板材料作为墙体的主要防火维护,其自身热工参数等对墙体耐火极限至关重要。国外规范对该类承重墙体耐火极限要求相对较低,导致研究大多针对石膏板覆面墙体,其他类板材覆面墙体的耐火试验则未见相关文献。国内仅有部分厂商进行该类非承重墙体的耐火检测,而鲜有对墙板材料、墙体构造的深入研究。且承重墙体的耐火试验具有耗资大、耗时长等弊端。

针对国内该类墙体常用板材开展研究,给出其热工参数模型并利用有限元模拟进行变参数分析,弥补试验数据量的不足,具有很好的实用性和经济性。

1 承重墙体试件耐火试验简介

结合已开展的冷成型钢承重组合墙体耐火试验建立立柱控制截面传热模型。试验情况介绍如下:

1.1 试件构造及测点布置

足尺墙体耐火试验试件规格为3 000 mm×3 000mm×137mm(长×宽×厚)。试件沿截面厚度方向的墙板、龙骨布置及拼缝位置,见图1所示。

C89龙骨截面规格为:89mm×50mm×13mm×0.9mm(腹板×翼缘×卷边×壁厚)

试验主要考察试件中部高度处A、B两切面(见图1)的温度。热电偶布置情况,见图2所示。

1.2 加载方式及试验结果

试件的力荷载加载情况及耐火极限,见表1所示。

2 硅酸钙板热工参数模型建立

利用文献[7、8]所述方法建立传热模型,并结合对国内某硅酸钙板粉末300 ℃以内热工参数的测定结果以及试件的试验现象与温度数据,给出了硅酸钙板的热工参数模型。测试热工参数所用仪器及方法同文献[7]。

2.1 硅酸钙板热工参数

2.1.1 硅酸钙板相对密度

硅酸钙板的常温密度经批量测试取1 000kg/m３。试验过程中,硅酸钙板处于受力状态且沿切面存在温度梯度,应力及内部水蒸气导致板材发生不同程度的爆裂、分层。由于是内层墙板,爆裂较为轻微,且受石膏板与龙骨夹持,基本没有脱落。板材由灰白色最终变为乳白色,开裂、分层情况见图3所示。

硅酸钙板300 ℃以下密度采用测定值拟合结果,300℃以上取为恒定值,为边长10cm方形试样,高温下不同温度恒温烘烤称重结果,见图4所示。

2.1.2 硅酸钙板比热容

硅酸钙板比热容测定值60 ℃前受仪器精度影响有较大浮动,因此取其合理部分进行拟合。其余部分根据比热容与密度的相关性进行拟合,见图5所示。

2.1.3 硅酸钙板导热系数

硅酸钙板的开裂、分层情况主要通过对导热系数的拟合来考虑。在散点值的基础上,考虑不同温度阶段裂缝、分层的发展过程,并结合数值模拟结果进行拟合,如图6所示。

2.2 数值模拟结果与分析

2.2.1 温度对比与分析

在后处理中,提取模型对应图2测温位置的节点温度,输出节点时间-温度曲线,与试件试验测点对比,见图7、图8所示。

由图7、图8可见:试验开始25min后,试件的a1、b1测点模拟值逐渐高于试验值,原因是单侧墙体耐火试验炉温、炉压精确控制较为困难;试件的a2测点试验值高于模拟值,原因是a2测点对应外层石膏板拼缝(见图1),试验中拼缝发生稍许敞开。而b2测点模拟值与试验值吻合较好,证明a2测点试验值异常。

总体上看,数值模拟结果与试验数据吻合较好,验证了硅酸钙板热工参数模型的合理性。

2.2.2 翼缘温度分布分析

提取试件的A切面龙骨各角点温度,见图9 所示。第92min受火侧翼缘角点(c2-c3)间温差4.6 ℃,背火侧翼缘角点(c4-c5)间温差7.5 ℃,原因是腹板能够很快地传导热量,卷边处则产生热量积聚。总体上看,受火侧、背火侧翼缘宽度范围内温度分布较为均匀,在热-结构耦合模拟中可以简化为等温。

3 墙体截面传热的变参数分析

利用上述试件的传热模型及给出的板材热工参数模型进行变参数分析,主要考察龙骨截面规格、各层墙板厚度、龙骨拼接方式及火源模型对龙骨温度的影响程度,并给出相应构造形式墙体的龙骨温度场数据。

3.1 龙骨截面规格影响分析

仅将原模型(如图1所示)的C89龙骨以C140龙骨代替(N1模型),主要考察龙骨腹板高度对其截面温度的影响。

由图10可见,C140龙骨与C89龙骨模拟点温度基本一致。说明由于龙骨壁薄,其截面高度对温度场分布没有明显影响。

C140龙骨截面规格为140mm×50mm×13mm×1.2mm腹板×翼缘×卷边×壁厚

3.2 内外层墙板厚度影响分析

仅将原模型(见图1)的各层墙板厚度由12mm增加到15mm(N2模型),以便考察墙板厚度对其截面温度的影响。

由图11可见,龙骨两侧各层墙板(外层石膏板、内层硅酸钙板)由12mm改为15 mm,对龙骨温度上升有显著延缓。如仅以原试件受火侧翼缘模拟温度最大值及受火侧、背火侧翼缘模拟温度差值为判断依据,则15 mm板厚组合墙体(N2模型)的耐火极限接近120min。说明墙板厚度对龙骨截面温度值(耐火极限)有很大影响。

3.3 龙骨拼接方式影响分析

仅将原模型(见图1)的单C89龙骨以双拼C89龙骨代替(双拼形式见图12),主要考察龙骨拼接方式对其截面温度的影响并对高温下双拼龙骨的截面温度分布进行简化。其中,N3模型为双C89腹板相拼龙骨,N4模型为双C89卷边相拼龙骨(点焊)。

3.3.1 双拼龙骨翼缘温度分布

N4模型与N3模型相比,受火侧翼缘平均温度高,背火侧翼缘平均温度低,龙骨两侧翼缘温差较大,见图13所示。如仅以受火侧翼缘模拟温度最大值及受火侧、背火侧翼缘模拟温度差值为判断依据,则N3模型的耐火极限大于N4模型。

与原模型相比,N3模型受火侧翼缘温度上升稍缓,背火侧翼缘温度基本一致。判断依据同上,则N3模型的耐火极限约为95min。说明龙骨截面拼接方式对温度场分布没有显著影响。

3.3.2 双拼龙骨截面温度分布简化

对N3、N4模型各模拟点温度进行计算,见图14所示。

对于N3模型,受火侧翼缘模拟点d1低于d2,背火侧翼缘模拟点d5高于d4;对于N4模型,受火侧翼缘模拟点e1高于e2,背火侧翼缘模拟点e5低于e4,同样说明热量在腹板中传导较快。

50min以内,龙骨模拟温度整体较低(见图7、图8),受火侧、背火侧翼缘温度差值与试验相差较大,因此计算结果不稳定,见图14。50 min以后,腹板温度逐渐呈线性变化(第90min,温度梯度误差在1%以内),翼缘温度也逐渐均匀(第90min,N3模型温度梯度误差在1% 以内,N4模型温度梯度误差在3.5%以内),因此可对双拼龙骨截面温度分布进行简化,见图15 所示。Tｈ为受火侧翼缘平均温度,Tｌ为背火侧翼缘平均温度,对腹板温度则进行线性插值。

3.4 火源模型影响分析

3.4.1 EU(0.03)火源模型建立

考虑房间的长×宽×高为3.0m×4.5m×3.0m,窗的高×宽为1.5m×1.2m,墙面及吊顶采用防火石膏板,房间可变火荷载为948 MJ/m２,恒定火荷载为130MJ/m２。房间开口系数约为0.03,室内维护边界热工参数取石膏板材性常温数据,根据文献[9]构造火源模型EU(0.03),见图16所示。

EU(0.03)曲线温度在t１=93 min达到峰值而后下降。t<tｐ时,EU(0.03)曲线温度高于GB/T 9978曲线温度;t>tｐ时,EU(0.03)曲线温度逐渐下降,而GB/T9978曲线温度继续上升。

3.4.2 EU(0.03)火源模拟及结果分析

在N2模型的基础上,以EU(0.03)曲线代替GB/T9978曲线作为火源模型(N5模型),主要考察火源模型对其截面温度的影响。

由图16可见:t≤tｐ时,N5模型的受火侧、背火侧翼缘温度均高于N2模型的相应测点,说明当N2模型的耐火极限tＮ２≤tｐ时,N5模型耐火极限tＮ５<tＮ２。

t>tｐ时,N5模型受火侧墙板表面温度(N5_a1)开始随EU(0.03)曲线逐渐下降,而内层墙板间温度(N5_a2)及龙骨温度(N5_a3)仍保持一段时间的上升。此过程可以定性地看出,在火源温度下降时墙板截面不同位置的温度变化过程,如:N5_a1、N5_a2、N5_a3分别在t１、t２、t３时刻达到最高温度(见图16)。说明墙体截面自受火侧向背火侧不同位置的降温过程具有滞后性。

由于各层板材(石膏板、硅酸钙板)热工参数模型是基于连续升温提出的,降温引起的板材热工参数变化仍有待研究。笔者认为,不适宜在EU(0.03)曲线的降温段定量讨论耐火极限。

4 结论

(1)笔者提出的基于试验的硅酸钙板热工参数模型经数值模拟计算,与试验结果吻合较好,可供同类构造墙体数值计算参考。

(2)在同类构造的冷成型钢承重墙体中:改变龙骨截面规格及龙骨拼接方式对截面温度场分布没有显著影响;改变外覆墙板厚度及火源模型,则有较大影响。

(3)单侧受火双拼龙骨立柱冷成型钢承重组合墙体中,组合龙骨截面温度分布可假设为:受火侧翼缘、背火侧翼缘宽度方向及相连的卷边温度一致,腹板高度方向温度为两侧龙骨温度线性插值。

(4)火源降温时,墙体截面自受火侧向背火侧的不同位置降温过程具有滞后性。

As a substitute structure system of timber structure,cold-formed steel(CFS)structure is built with compositefloor slabs,load-bearing walls and etc.These elementsare lined with plasterboard on both sides of the light steelframe.Gypsum plasterboard,bolivian magnesium board,calcium silicate board and etc.are used as wallboard.Theframe is assembled by C-shaped studs and U-shapedtracks,and connected by self-drilling screws.CFS struc-ture has the advantages of low material consumption andrecoverable.It conforms to the present national policy,which limits the unreasonable over production capacityand encourages people to use more steel.

China has large population,while the land resourcesare limited.It will be more in line with national condi-tions to introduce the CFS structures to China and applythem to multi-story buildings.The application exampleshave been more mature abroad.For example,in UnitedStates,a six-story apartment has been built with CFSstructure.In Ohio,six-story and even nine-story CFSdwellings have been built.Canada has built six-story hou-ses and eight-story hostels using CFS structure.

Unlike foreign countries,we emphasize the passivefire protection performance of buildings so as to ensurethe fire resistance rating of different building elements toachieve effective fire separation and confine fire withinone dwelling unit.Therefore,as the key member to bearvertical and lateral loads,load-bearing wall should meetmore strict fire resistance requirements.As the importantfire prevention enclosure,the thermodynamic parametercurve of wallboard material is an essential factor for fireresistance rating of the wall.The requirements of foreigncodes for fire resistance rating of CFS load-bearing wallsare relatively low,which leads to the fact that most stud-ies about the CFS wall are lined with gypsum plasterboardsheathing.The literature about fire resistance tests of theload-bearing wall with other types of wallboard sheath-ings is rarely seen.The performance of CFS non-load-bearing walls under fire conditions was analyzed by onlyfew domestic manufacturers,and little in-depth studies ofthe wall materials and wall construction was conducted.Fire resistance test of load-bearing walls was costly andtime-consuming.

Therefore,study was done on the commonly usedwallboard in domestic CFS load-bearing wall.The ther-modynamic parameter curve of wallboard material wasproposed and different parameters of the wall by FEMwere analyzed.It would provide more data and would bepractical and economical.

1 Brief introduction to fire test of load-bearing wall

An FEM of crucial cross section was developed basedon the previous fire resistance tests.The overview of testis described as below:

1.1 Cross section geometry and temperature measuringpoints

The specification of the specimen was 3 000 mm×3 000mm×137 mm (width×height×thickness).Thewallboards on each side of the wall,location of the studsand seams are shown in Fig.1.

Note:The specification of C89:89mm×50mm×13mm×0.9mm(web×flange×lip×thickness).

The key target of the test was to observe the temper-ature profile of section A and B,which located at the halfheight of the stud(Fig.1).The arrangement of thermo-couples is shown in Fig.2.

1.2 Loading and test results

Loading conditions and fire resistance rating of speci-men are shown in Table 1.

Note:"load ratio" means the ratio between load val-ue applied to the specimen at fire resistance test and theultimate bearing capacity of the specimen with same con-figuration at static load test.

2 The development of thermodynamic parameter curve ofCalcium silicate board

The thermodynamic parameter curve of Calcium sili-cate board was developed using the methods discussed inliterature [7]and [8],thermodynamic parameters ofpowder of Calcium silicate board in China within 300 ℃and the phenomenon and temperature of full-scale speci-men test under standard fire conditions.The instrumentsand methods used to test the thermodynamic parametersof powder were the same as discussed in literature[7].

2.1 The thermodynamic parameter curve of Calcium silicate board

2.1.1 The relative density of Calcium silicate board

The density of Calcium silicate board at ambient tem-perature is 1 000kg/m３.During the test,Calcium silicateboard was bearing pressure and a temperature gradientcross the board section could be found,so the internalstress and steam caused a varying degree of explosivespalling and layered plates.The explosive spalling wasrelatively minor and did not fall down,as it was the baselayer and located between the gypsum plasterboard andstuds.Color of the layer became milky white from gray e-ventually.The characters of cracking and layered areshown in Fig.3.

Therefore,the density of Calcium silicate board be-low 300 ℃ was obtained from the fitting results of testdata.It was taken as constant equaling to the results of10cm square sample tests above 300 ℃.The test com-pared the weight of sample at ambient and elevated tem-perature conditions(Fig.4).

2.1.2 Specific heat of Calcium silicate board

Specific heat of Calcium silicate board has major fluc-tuations due to the influence of instrument accuracy with-in 60℃.Only the reasonable test data was chosen for fit-ting.The rest part of specific heat was based on the cor-relation of specific heat and density as shown in Fig.5.

2.1.3 Thermal conductivity of Calcium silicate board

The cracking and layered phenomenon of Calcium sil-icate board was considered by the fitting of thermal con-ductivity.The fitting was based on scatter value,the de-velopment process of cracking and layered phenomenonunder different temperature stages and the results of nu-merical simulation(Fig.6).

2.2 Results and analysis of simulation

2.2.1 Comparison and analysis of temperature results

In the post26module of ANSYS,got the node tem-perature from model.The node positions referred tomeasuring points (Fig.2).The comparing results be-tween node time-temperature curve and test result curveare shown in Fig.7and Fig.8.

Fig.7and Fig.8showed that after 25minutes,atpoint a1and b1,the simulation value became graduallyhigher than the test.The reason was that it was difficultto control precisely the furnace temperature and pressurewhen the wall had only one side exposed to fire.At pointa2,test results were higher.The reason was that a2lo-cated at the seam of face layer and the base layer and theseam slightly opened at the later stage of the test.Simu-lation results agreed well with the test at point b2.It alsoproved that the test result was abnormal.

Overall,the results of the numerical simulation werein good agreement with that of the test.It proved that thethermodynamic parameter curve of calcium silicate boardproposed in this paper was quite reasonable.

2.2.2 Flange temperature profile analysis

The temperature of every corner points on section A(Fig.9)was selected.At 92 min,the temperaturedifference between two corner points(c2and c3)of theflange exposed to fire was 4.6 ℃.The temperaturedifference between two corner points(c4and c5)of theflange at the unexposed side was 7.5 ℃.The reason wasthat the web could quickly transfer the heat,while thelips would cause heat accumulation.Overall,it was quiteuniform for the temperature profile of the flanges(alongthe width)both at the fire side and the ambient side.Theflange temperature could be considered as constant whenapplied to thermal-structure coupled simulation.

3 Parameter analysis of heat transfer simulation of crosssection of wall

Parameters have been analyzed based on the FEMand thermodynamic parameter curve of calcium silicateboard proposed above.The effect of different depth ofweb,the joint ways of C-shaped studs,the thickness ofplasterboard and fire scenario on the degree of tempera-ture rise of C-shape section was analyzed.The tempera-ture profile of the cross section of studs has been provid-ed.

3.1 Analyzing the effect of different depth of stud

The stud of original FEM from C89 (Fig.1)waschanged to C140 (FEM N1)and the effect of differentdepth of web on the temperature profile of the stud sec-tion has been analyzed.

Fig.10shows that the temperature of the same simu-late points at the same locations of section C140and C89agreed well.It indicats that for the thin-wall stud,thedepth of web did not affect temperature profile of its crosssection significantly.

Note:The specification of C140:140mm×50mm×13mm×1.2mm(web×flange×lip×thickness).

3.2 Effect analysis of plasterboard thickness

Two layers of plasterboard on both sides of the studwere increased from 12 mm (Fig.1)to 15 mm (FEMN２),and the effect of plasterboard thickness on the tem-perature profile of the section has been analyzed.

Fig.11shows that the thickness of plasterboard in-creased from 12mm to 15mm would delay the tempera-ture rise of the stud significantly.If only the maximumtemperature of simulated point at the flange exposed tofire and the difference value of measuring points of flangesbetween the fire side and ambient side were considered asthe criteria,the fire resistance rating of the compositewall(FEM N2)would reach 120 min.It indicated thatthe thickness of plasterboard affected the temperatureprofile of cross section(fire resistance rating)significant-ly.

3.3 Effect analysis of joint ways of C-shaped studs

Studs of original FEM from single C89(Fig.1)werechanged to double C89(Fig.12),the effect of joint waysof studs on the temperature profile of the cross section ofstud was analyzed.The simplified temperature profile ofdouble studs at elevate temperature has been provided.FEM N3was constructed with double C89studs and wasconnected to web with tapping screw and FEM N4 wasconstructed with double C89studs that were welded tothe lips.

3.3.1 Temperature profile of flanges of double studs

Comparing with FEM N3,the average temperatureof flanges of FEM N4exposed to fire was much higher,while the temperature at the unexposed side was muchlower.The temperature difference of the flanges of studsat exposed and unexposed sides were quite big,as shownin Fig.13.

If only the maximum temperature of simulated pointat the flange exposed to fire and the difference value ofmeasuring points of flanges between the fire side and am-bient side were considered as the criteria,the fire resist-ance rating of FEM N3was longer than that of the FEMN4.

Comparing with original FEM,the temperature ofN3flanges exposed to fire rose slower,while the temper-ature of the unexposed side was basically the same.Ac-cording to the same criteria discussed above,the fire re-sistance rating of FEM N3was about 95min.It indicatedthat the joint ways of studs did not affect the temperaturefield of cross section significantly.

3.3.2 The simplified model of cross section of double-studs

The temperature of measuring points of FEM N3andN4was calculated,as shown in Fig.14.

For FEM N3,the temperature of measuring point d1at exposed side was lower than that of d2and point d5atunexposed side was higher than that of d4.For FEM N4,the temperature of measuring point e1at exposed side washigher than that of e2and point e5at the unexposed sidewas lower than that of e4.These results indicated thatthe heat conduction in the web was quite fast.

Within the first 50 min,the simulated temperatureof the section was quite low (Fig.7and Fig.8)and tem-perature difference of flanges at the exposed side and un-exposed side was quite different from that of the test re-sults.Therefore,the calculation results were not stable,as shown in Fig.14.After 50 min,the temperature ofweb gradually changed linearly and the temperature offlanges also became gradually uniform.At 90 min,thetemperature gradient error of web was within 1%and 1%for flanges of FEM N3and 3.5%for flanges of FEM N4.It indicated that the temperature profile of two types ofsections could be simplified,as shown in Fig.15.Tｈwasthe average temperature of flange exposed to fire.Tｌwasthe average temperature of flange exposed to ambientside.Temperature along the web was calculated by linearinterpolation.

3.4 Effect analysis of design fire scenario

3.4.1 Design fire scenario EU(0.03)

The size of the room was 3.0 m×4.5 m×3.0 m(length×width×height)and the window size was 1.5m×1.2m (height×width).Walls and ceiling of the roomwere covered with gypsum plasterboard.The live load inthe room was 948 MJ/m２and dead load was 130 MJ/m２.The opening factor of the room was 0.03.The parameterof gypsum plasterboard at ambient temperature was usedas the thermodynamic parameter of the enclosure.Thedesign fire scenario EU(0.03)was based on literature[9],as shown in Fig.16.

The temperature curve of EU(0.03)reached thepeak at 93 min (tｐ=93 min)and decreased after that.When t<tｐ,the temperature of EU(0.03)was higherthan that provided in GB/T 9978;When t>tｐ,the tem-perature of EU(0.03)decreased gradually,but the tem-perature based on GB/T 9978continued to rise.

3.4.2 EU(0.03)simulation and results analysis

The fire scenario of FEM N2was changed from GB/T 9978curve to EU(0.03)curve(FEM N5).The effectof fire scenario on the temperature profile of the stud sec-tion was analyzed.

Fig.16shows that when t≤tｐ,the temperature offlanges at fire side and ambient side were all higher thanthat of the measuring points of FEM N2.It indicated thatwhen the fire resistance rating of FEM N2 tＮ２≤tｐ,thefire resistance rating of N5 tＮ５<tＮ２.

When t>tｐ,the surface temperature (N5_a1)ofplasterboard at the fire side started to decrease with EU(0.03)curve gradually,but the temperature of innerplasterboard(N5_a2)and stud (N5_a3)still kept risingfor a period of time.The temperature change across dif-ferent part of the wall could be gotten qualitatively underthe process of fire decay.For example,the measuringpoints of N5_a1,N5_a2,N5_a3reached peak tempera-ture at time t１,t２,t３respectively (Fig.16).It indicatedthat the temperature descending process of different partof the wall section,from the fire side to the ambient side,had a time lag.

However,considering the thermodynamic parametercurve of wallboard(gypsum plasterboard and calcium sili-cone board)were proposed based on continuous heatinghypothesis,the changes of thermodynamic parameters ofwallboard caused by temperature descending still needs tobe further studied.Therefore,it is inappropriate to dis-cuss the fire resistance rating of wallboard at the coolingprocess of EU (0.03)curve.

4 Conclusions

(1)In this paper,the thermodynamic parametercurve of calcium silicate board proposed based on testphenomena is in good agreement with test results accord-ing to simulation results,and can be used to simulate thewall with similar construction.

(2)For CFS load-bearing wall with similar construc-tion,changing of depth of web and the joint ways of studsdid not affect the fire resistance rating of the wall signifi-cantly.However,if the thickness of plasterboard and thedesign fire scenario were changed,the effect was signifi-cant.

(3)For CFS load-bearing wall with double-studs andone side fire exposure,the temperature profile of double-studs section can be simplified.The temperature profileof the flange and lips can be simplified as uniform temper-ature.Temperature profile along the web was calculatedby linear interpolation of temperature between fire sideand ambient side flanges.

变参数回归模型 第8篇

1 文献综述

关于技术对于碳排放的影响, 国内外学者运用不同的方法进行了研究, 取得了不少的成果。Nordhaus[1]分析了外生技术进步下经济增长对环境的影响, 之后又构建了气候变化和经济的动态综合模型 (DICE) , 分析经济增长和环境的相互关系;Grubler和Messner[2]通过建立技术内生化模型, 利用学习曲线, 提出为更好地完成减排目标必须提高研发投入和技术学习能力;Kemfert[3]的研究表明在没有技术进步的情况下减排初期会导致产出量大幅度减少, 从而降低整体的福利, 在有技术进步的情况下产出降低的幅度会减少;Gerlagh[4]认为技术进步一方面可以降低碳价格, 从而降低强制性减排的负担, 另一方面认为存在技术进步情境下的减排可能会产生学习收益, 从而降低减排成本;Smulders和Di Maria C[5]将技术进步作为内生变量分析技术进步对碳排放边际成本的影响;林伯强等[6]对中国1978—2007年的时间序列数据, 采用STIRPAT模型分析了中国人均二氧化碳排放的影响因素;陈劭锋等[7]在揭示科技进步与应对气候变化国家发展趋势的基础上, 探讨了科技进步与碳排放演变之间的关系;张永军等[8]估算了中国1978—2008年碳生产率, 并用拉氏分解法分析了技术进步、产业结构和能源消费结构等影响碳生产率的主要因素, 提出要发挥技术效率在节能减排中的主导作用。

上述研究技术进步与碳减排关系的学者们在对技术进步指标的选择上, 大多选择全要素生产率或专利数量, 并且多是从静态角度去考察当期技术进步对碳排放的影响, 未分析技术进步的累积影响, 而国际上通常采用研究与试验发展 (R&D) 活动的规模和强度指标来反映一国的科技实力和核心竞争力, 认为R&D活动是技术进步的源泉、是促进技术进步最直接的因素, 因此, 我们选择R&D经费内部支出来反映科技进步的水平。另外, 从现实情况来看, 科技投入存在着累积效应, 即每年投入的R&D经费不仅会对当年的产出造成影响, 而且以往研究开发所产生的知识和经验的积累对后续的开发和产出产生积极的影响, 因而本文拟采用胡恩华等[9]提出的以2年作为累积期的做法, 对折算后的科技活动经费内部支出进行累积, 并以累积值作为模型的解释变量来揭示技术进步对碳排放的长期影响及动态关系。

2 模型建立与数据来源

2.1 模型建立

Ehrlich和Holden[10,11]提出IPAT方程用来考察人口对环境压力的影响, 其公式为:I=人口 (P) *富裕程度 (A) *技术 (T) , 指出每一个社会都应该在最有可能改进的地方加以改进, 但IPAT模型存在一些局限性, 即假设各驱动因素和环境因素之间存在线性关系。Richard等[12]在该模型的基础上建立了当前广为应用的STIRPAT (Stochastic Impacts by Regression on Population, Affluence, and Technology) 模型, 其估计方程为:

该模型将二氧化碳排放 (I) 的三个驱动力因素分解为:人口 (P) 、经济规模 (A) 、技术 (T) , 有效避免了IPAT模型无法实现的自变量和因变量之间的非线性关系计量。本文在分析技术进步对二氧化碳排放的影响时, 选择借鉴STIRPAT模型, 在用面板数据进行计量分析时将该模型转化为对数形式, 设计二氧化碳排放与技术进步等变量关系的模型如下:

该模型中, CO2表示各地区二氧化碳排放量, T代表技术进步水平, GDP代表经济发展水平, P表示各地区人口, 下标i指各省份, t指各年份。

2.2 数据来源与处理

(1) 二氧化碳排放量 (CO2) :根据我国各地区能源消费总量及各种能源的碳排放系数进行测算。根据2007年政府间气候变化问题小组 (IPCC) 第4次评估报告, 温室气体增加的主要来源是化石能源燃烧, 而化石能源主要包括煤炭、石油和天然气, 因此, 本文也主要选择这3类能源的消费量作为计算碳排放量的依据。《中国能源统计年鉴》提供了能源消费的具体实物消费量数据, 在计算时需要将各类能源的实物消费量转换为标准煤消费量, 折算系数分别为:1千克原煤=0.7143千克标准煤, 1千克原油=1.4286千克标准煤, 1立方米天然气=1.3300千克标准煤 (参考2011年《中国能源统计年鉴》附录中的折算系数) ;碳排放系数来源于国家发展和改革委员会能源研究所2003年发布的《中国可持续发展能源暨碳排放情景分析》, 具体如表1所示。我国各地区碳排放量的计算公式如下:

此式中, CO2it为i地区第t年的碳排放量, Ejit为i地区第t年第j种能源的标准煤消费量, ηj为第j种能源的碳排放系数。据此可测算出中国各地区1997—2010年的碳排放量, 单位为万吨碳。

注:数据来源于《中国可持续发展能源暨碳排放情景分析》

(2) 技术进步:采用我国各地区R&D经费内部支出总额来表示, 并以2年为累积期进行累积, 体现各地技术进步的水平, 单位为万元。

(3) 经济发展:用我国各地区的人均GDP表示, 反映经济体的经济发展水平, 单位为元。

(4) 人口规模 (P) :用我国各地区的总人口数表示, 单位为万人。

本文以我国1997—2010年的省际面板数据为研究对象, 其中由于西藏地区多年数据缺失, 因而在样本中予以剔除。本文的原始数据来源于《中国统计年鉴》、《中国科技统计年鉴》、《中国能源统计年鉴》、《新中国六十年统计资料汇编》。在对原始数据的处理中, 由于各地经济发展很不平衡, 故各地区的人均GDP和R&D经费内部支出拟采用各地区的国内生产总值指数 (以1997为基期) 来折算, 以减少误差;为了有效消除异方差的同时不改变变量的主要特征, 我们将各个变量进行对数化处理, 分别记为LNCO2、LNK、LNGDP、LNP。

3 实证分析

3.1 面板数据模型的估计

面板数据是包含时间与截面两维性质的数据, 以面板数据为基础建立的模型能够同时反映研究对象在时间和截面单元两个方向上的变化规律及不同时间、不同单元的特性, 因而在计量经济学、社会学等诸多领域已有较广泛的应用[13]。单方程面板数据模型的一般形式为:

其中:N为截面单元个数, T为时序期数, xit为1×K向量, βi为K×1向量, K为解释变量的数目, 误差项μit均值为零, 方差为δu2。根据截距项向量α和系数β中各分量的不同限制要求, 此模型常用的形式有如下三种:

情形1:αi=αj, βi=βj, 即混合模型yit=α+xitβ+μit (1)

情形2:αi≠αj, βi=βj, 即变截距模型yit=αi+xitβ+μit (2)

情形3:αi≠αj, βi≠βj, 即变系数模型yit=αi+xitβi+μit (3)

在现实的经济分析中, 为了避免模型设定的偏差, 改进参数估计的有效性, 需要对样本数据究竟属于哪一种面板数据模型形式进行检验, 通常用协方差分析检验来完成, 即F检验。主要检验如下两个原假设:

假设1:斜率在不同的截面样本点上和时间上都相同, 但截距不相同。

假设2:截距和斜率在不同的截面样本点和时间上都相同。

检验假设2的F统计量为:

其中, S1和S3分别为模型 (3) 和模型 (1) 的残差平方和。如果计算所得到的统计量F2的值不小于给定置信度下的相应临界值, 则拒绝假设H02, 继续检验假设H01;反之, 则认为样本数据符合不变系数模型, 即模型 (1) 。

检验假设1的F统计量为:

其中, S2为模型 (2) 的残差平方和。若计算所得到的统计量F1的值不小于给定置信度下的相应临界值, 则拒绝假设H01, 用变系数模型, 即模型 (3) ;反之, 则认为样本数据符合变截距模型, 即模型 (2) 。

模型 (2) 和 (3) 又可以分为固定影响和随机影响两种情况, 因此, 不论采用模型 (2) 还是模型 (3) , 都要进一步确定影响方式。这一过程常通过Hausman检验来完成。

假设3:个体影响与回归变量无关。

H03:应该建立随机影响回归模型。

检验假设3的统计量为:

其中:表示个体固定效应回归模型参数的估计量, 是离差OLS估计量;表示个体随机效应回归模型参数的估计量, 是可行的GLS估计量。若χ值小于临界值, 则接受假设3, 说明是一致估计量, 两者的差异较小, 应该建立个体随机效应回归模型;若χ值大于临界值, 则拒绝假设3, 说明是一致估计量而是非一致估计量, 两者差异较大, 应该建立个体固定效应回归模型。

为了更好地考察不同地区二氧化碳排放的影响因素, 本文按照全国、东部、中部、西部划分方法来进行实证分析。其中, 东部包括北京、天津、河北、辽宁、上海、江苏、浙江、福建、山东、广东、广西、海南12个地区;中部包括山西、内蒙古、吉林、黑龙江、安徽、江西、河南、湖北、湖南9个地区;西部包括重庆、四川、贵州、云南、陕西、甘肃、青海、宁夏、新疆9个地区, 西藏因缺乏数据未纳入分析。

基于Eviews6.0软件, 利用本文面板数据得出变系数模型、变截距模型、混合模型的残差平方和S1、S2、S3, 计算出F1、F2统计量的值并与临界值进行比较, 各个变量的具体值如表2所示。

结果显示, 在5%的显著性水平下, 无论是全国范围数据、还是分区域数据, F1、F2的统计值都大于临界值, 因此应该选择变系数模型;而且根据Hausman的检验, 支持建立固定效应模型, 这与实际操作中的经验做法一致, 即如果预期以面板数据模型推断样本空间的经济关系, 则模型设定为固定效应模型合理些。

因此, 根据以上F检验及Hausman检验结果, 本文采用截面加权方法进行普通最小二乘估计, 进而建立固定效应变系数模型, 模型估计结果如表3所示。

注:1) 限于篇幅, 省略了t值统计量及相应的概率值;2) ***、**、*分别表示结果在1%、5%、10%水平上显著

3.2 回归结果分析

从回归结果来看, 全国30个省市区中, 辽宁、海南、甘肃、新疆地区的技术进步系数未通过显著性检验, 其余地区都通过了5%水平显著性检验;福建、黑龙江、江西、甘肃地区的人口系数未通过显著性检验, 其余地区都通过了5%水平显著性检验;所有地区的人均GDP系数都通过了5%水平显著性检验。调整的R2系数达到0.994104, F检验值为594.695, 表明模型的解释能力强;同时大部分系数t检验值偏大, p值偏小, 表明多数方程的拟合优度较好。

从系数估计结果来看, 发现技术进步系数通过显著性检验的地区中, 大多数地区的技术进步系数为负值, 表明技术进步可以在一定程度上抑制碳排放强度;就人均GDP而言, 模型中30个地区的系数都为正, 说明经济增长对碳排放具有正向的促进作用;而对于人口弹性而言, 全国有的地区人口系数为正、有的为负, 说明人口增长对环境的影响是双向的。

分地区来看, 东部地区的技术进步对碳排放的抑制作用最为显著, 要远高于中部、西部地区, 说明东部省份的能源利用效率相对较高;而中部地区的技术进步系数又高于西部, 这与目前中国的技术空间分布特征是相符的。但相比于经济增长与人口弹性, 技术进步对碳排放的影响系数最小, 原因在于技术进步对碳排放的影响要考虑其时滞和积累问题, 技术创新后的扩散和应用要经过一段时间的积累才能体现出来, 短期内很难将技术进步快速转变为碳排放的减少。就经济发展而言, 人均总产出对碳排放的影响比其他两个指标明显大很多, 说明经济发展对环境的影响非常大, 尤其中部地区相比东、西部而言, 其影响更为显著一些, 主要原因在于中部地区还处在粗放型的发展阶段, 作为承接产业转移的中间地带, 承接的大部分产业都是重工业, 使得该地区的二氧化碳排放量相对增加, 因而对于发展起步晚于东部的中部地区, 经济发展方式优化是影响现在和以后中部碳排放量大小极其重要的因素。就人口增长对碳排放的影响而言, 东部和中部地区基本上呈现正向关系, 即人口的增长对环境产生了压力, 增加了能源消费, 导致二氧化碳排放增加;而西部大部分地区这二者呈现反向关系, 表明对于西部地区而言, 增加人口反而会降低碳排放的数量, 原因在于对于西部地区, 人口增长有利于人口和经济集聚, 同时会促进技术改革, 从而减轻了对环境的负面影响。

3.3 变量的单位根及协整检验

3.3.1 单位根检验

一般而言, 传统的计量模型有一个基本前提, 即每个变量的时序数据必须是平稳的, 否则容易产生“伪回归”, 其统计推断结论也是值得怀疑的, 因而, 为确保以上面板数据回归结果的准确性, 必须进行数据的平稳性检验。根据所有截面序列是否具有相同的单位根, 可以将检验单位根过程的方法分为同质单位根检验法 (Common Root Test) 与异质单位根检验法 (Individual Root Test) 。同质单位根检验的主要方法有LLC、Breitung, 异质单位根检验的方法主要有IPS、Fisher-ADF、Fisher-PP。本文利用上述方法对中国的碳排放、技术进步、人均GDP、人口规模这4个变量进行单位根检验, 检验结果如表4所示。

根据表4可知, 对于原假设存在单位根的各种检验方法, 变量LNCO2、LNK、LNGDP、LNP在1%的显著水平下不能拒绝有单位根的原假设, 说明这4个变量是不平稳的;其一阶差分系列在1%的显著性水平下, 均拒绝原假设, 故4个变量均具有一阶单位根。因此, 可以判定各个面板的所有变量都是一阶单整的。由于面板数据的不稳定性, 需要继续分析碳排放与各个变量的协整关系, 进而分析模型的长期关系。

3.3.2 变量的协整检验

协整检验是检验变量之间是否存在长期均衡关系的方法。Pedroni以回归残差构造了7个统计量, 其中4个是用联合组内维度 (Within Dimension) 描述, 记为Panel v、Panel rho、Panel PP、Panel AD, 另外3个用组间维度 (Between Dimension) 描述, 记为Group rho、Group PP、Group ADF。基于Pedroni的证明, 在小样本中, Panel ADF和Group ADF检验效果最好, 而Panel v和Group rho效果最差;当结果不一致时, 以Panel ADF和Group ADF为准。考虑到本文数据的小样本性质, 主要采用Panel ADF和Group ADF来判断变量之间是否存在协整关系, 但基于稳健性的考虑, 同时也参考Panel PP、Group PP统计量, 并采用Kao提出的面板协整检验方法。检验结果如表5所示。

从表5分析结果看, Panel ADF、Group ADF、Panel PP、Group PP统计量都通过了1%的显著性水平检验, 并且根据Kao协整检验, 也通过了1%的显著性水平检验, 由此表明碳排放量与技术进步水平、经济发展、人口规模等变量之间存在长期稳定的内生经济关系。

4 结论

(1) 技术进步、经济发展和人口规模与中国的碳排放存在长期稳定的协整关系, 技术进步会在一定程度上抑制碳排放, 人均GDP的增长会带来碳排放的增加, 人口数量对碳排放的影响是双向的, 按影响大小依次为:人均GDP、人口规模和技术进步。