三角函数分解

2024-06-26

三角函数分解(共5篇)

三角函数分解 篇1

一、矩阵的三角分解

1. 定义

如果方阵A可分解成一个下三角矩阵L和一个上三角矩阵U的乘积, 则称A可作三角分解或LU分解。如果方阵A可分解成A=LDU (1.1) , 其中L是单位下三角矩阵, D是对角矩阵, U是单位上三角矩阵, 则称A可作LDU分解。

2. A的LDU分解和LU分解求法

(1) A的LDU分解求法:取A的一阶主子式, 作LDU分解A1= (1) (a11) (1) , L1= (1) , D1= (a11) , U1= (1) , 用 (1) 式- (4) 式确定v1, l1, 从而L2=βl1L110β, D2=β0D1 d20β, U2=β0V1 1v1β, A2=L2D2U2。然后重复使用 (1) - (4) 式得到A的顺序主子式的LDU分解。Ak=LkDkUk, k=1, 2, …, n, 当k=n时, 即完成A的LDU分解。

(2) A的LU分解求法:先按照 (1) 的步骤求出A的LDU分解, 再令DU的乘积为U’, 即求出A的LU分解。

3. LU分解的推广和改进

矩阵A的LU与LDU分解都需要假设A的前n-1阶顺序主子式非零。如果这个条件不满足, 可以给A左或右乘以置换矩阵, 就把A的行或列的次序重新排列, 使这个条件满足, 从而就有如下的带行交换的矩阵分解定理。

定理1设A是n阶非奇异矩阵, 则存在置换矩阵P使PA的n个顺序主子式非零。该定理的证明可在矩阵论教科书中查到。

4. LU分解在解方程组中的应用

设A是n阶非奇异矩阵, 则存在置换矩阵P使PA=LDU=LU’ (1.3) 。

其中L是单位下三角矩阵, D是对角矩阵, U是单位上三角矩阵, U’是上三角矩阵。

如果方程组AX=b (1.3) , 系数矩阵A是n阶非奇异矩阵且△k≠0, (k=1, 2, …, n-1) , 则存在三角分解A=LU.于是得到与原方程组同解的以三角矩阵为系数矩阵的方程组βUx=yLy=b (2) (1) , 方程组的 (1) 式解出代入 (2) 式解出, 这就是线性方程组的三角分解法。如果方程组的系数矩阵A的某个顺序主子式△k=0, (k=1, 2, …, n-1) , 可用定理1的推论考虑与其同解的方程组PAX=Pb (1.4) , 即也可用三角分解法求解。

AX=b的解。

解:A的LU分解为

则Ly=b的形式为

于是可求AX=b解向量为

二、矩阵的QR分解

定义如果非奇异矩阵A能够化成正交矩阵Q与非奇异上三角矩阵R的乘积, 即A=QR, 则称A可QR分解。

定理2如果A是n阶非奇异矩阵, 则存在使A=QR成立 (2.1) , 且除去相差一个对角元素的绝对值 (模) 全等于1的对角矩阵因子外, 分解是唯一的。

证记A的列向量α1, α2, …, αn方阵A是非奇异矩阵, 所以列向量α1, α2, …, αn线性无关, 对A的列向量α1, α2, …, αn, 正交化得正交向量β1, β2, …, βn, (α1, α2, …, αn) = (β1, β2, …, βn) K,

于是R=diag (β1, β2, …, βn) K, 则Q是正交矩阵, R是非奇异上三角矩阵, 且有A=QR。

例2试用Schmidt正交化的方法求矩阵A=β121212221β的QR分解。

模式分解保持函数依赖的意义 篇2

规范化需要进行关系的模式分解,评判准则是无损连接性和保持函数依赖,分解必须具有无损连接性而尽量保持函数依赖[2]。一个较大关系被分解为若干个较小关系时,如果分解具有无损连接性,则说明将所有较小关系自然连接可以还原为较大关系,信息不会丢失。如果分解能够保持函数依赖,则说明更新任意一个较小关系时,为了判断更新数据是否遵从原有的函数依赖,只需要考虑该较小关系的函数依赖,而不必连接多个较小关系。一般的国内文献重点研究的是范式定义和模式分解相关算法[3],对于分解保持函数依赖的意义则较少提及[2,3,4]。有些文献虽然讨论了保持函数依赖的意义,但是没有给出严格的证明和易于理解的实例[4,5,6]。本文证明了保持函数依赖的一个定理,说明了模式分解保持函数依赖的意义是在减小数据冗余度的同时不会影响数据更新的性能,并通过实例进一步印证了该结论。

1 相关定义

文中讨论涉及到的规范化概念有模式分解、无损连接性和保持函数依赖,定义介绍如下[2]。

定义1 关系模式R<U,F>的一个分解是指ρ={R1<U1,F1>,R2<U2,F2>,…,Rn<Un,Fn>},其中U=i=1nui,并且不存在UiUj,1≤i,jn,ij,FiFUi上的投影。

定义2 ρ={R1<U1,F1>,R2<U2,F2>,…,Rn<Un,Fn>}是关系模式R<U,F>的一个分解,若对于R<U,F>的任一关系r均有r=i=1nπRi(r)成立,则称分解ρ具有无损连接性,简称ρ为无损分解。

定义3 若F+=(i=1nFi)+,则R<U,F>的分解ρ={R1<U1,F1>,R2<U2,F2>,…,Rn<Un,Fn>}保持函数依赖。

关于模式分解、无损连接性和保持函数依赖的重要结论有[2],若要求分解具有无损连接性,那么一定可以达到4 NF。若要求分解保持函数依赖,那么一定可以达到3 NF,不一定能够达到BCNF。若要求分解既具有无损连接性,又保持函数依赖,那么一定可以到达3 NF,不一定能够达到BCNF。

无损连接性和保持函数依赖是关系模式分解好坏的两个评判准则,分解应该具有无损连接性,同时应该保持函数依赖,原因在于:一个较大关系被分解为若干个较小关系时,如果分解不具有无损连接性,则将所有较小关系自然连接不能够还原为较大关系,丢失了信息。如果分解能够保持函数依赖,则在更新任意一个较小关系时,为了判断更新数据是否遵从原有的函数依赖,需要连接多个较小关系,从而影响了更新性能。

2 保持函数依赖的意义

对于模式分解具有无损连接性的意义,国内外文献均有所介绍,并给出了证明和算法[2,3,4,5,6]。但是对于模式分解保持函数依赖的意义,并没有出现严格的证明和易于理解的实例[4,5,6]。文中证明了保持函数依赖的一个定理,说明了模式分解保持函数依赖的意义是在减小数据冗余度的同时不会影响数据更新的性能。

定理1 ρ={R1<U1,F1>,R2<U2,F2>,…,Rn<Un,Fn>}是关系模式R<U,F>的一个保持函数依赖的分解,关系模式Ri的关系表示为ri,对其中一个关系rm进行更新得到新关系,若满足函数依赖集Fm,则i=1m-1rirmi=m+1nri满足函数依赖集F

证明:

(1)新关系rm满足函数依赖集Fm,其他未更新的关系ri满足相应的函数依赖集Fi,则i=1m-1rirmi=m+1nri满足所有函数依赖集的并i=1nFi,根据Armstrong公理系统[2],它也满足(i=1nFi)+

(2)根据保持函数依赖的定义3,有F+=(i=1nFi)+,由(1)知i=1m-1rirmi=m+1nri也满足F+,而FF+,则有i=1m-1rirmi=m+1nri满足函数依赖集F,定理得证。

通俗的解释,该定理说明,一个较大关系被分解为若干个较小关系时,减小了数据冗余度,如果分解能够保持函数依赖,则在更新任意较小关系时,只要能遵从自己单独的函数依赖集即可,不会影响数据更新的性能。因为所有较小关系自然连接后,产生的较大关系肯定够满足原有的函数依赖,不会丢失任何函数依赖信息。相反的,如果分解不保持函数依赖,则有i=1m-1rirmi=m+1nri满足的函数依赖集(i=1nFi)+F+的真子集,那么就丢失了函数依赖信息,需要重新检测函数依赖集F+-(inFi)+i=1m-1rirmi=m+1nri是否成立,如果成立,则说明rm的更新操作是合法的,否则更新操作是非法的,不能够进行。更为严重的是,该检测需要连接若干个较小关系模式,从而降低了数据更新的性能。

该定理说明了模式分解中保持函数依赖的意义,在实际的数据库系统中也能够得到具体的应用。在关系数据库中,函数依赖有多种表达方式[5]:主码约束,check约束,断言和触发器。在每次数据更新时,检查这些函数依赖的开销很大,因此,将数据库设计得能够高效检测函数依赖是有必要的。如果在检测函数依赖时仅需要考虑一个关系,那么开销就很低。相反,如果必须先连接关系,再检测函数依赖,显然是不切实际的。

3 实例分析

例1 已知关系模式R<U,F>,其中U={Sno,Sdept,Mname},F={SnoSdept,SdeptMname}。R<U,F>的元组语义是学号为Sno的学生正在Sdept系学习,其系主任姓名为Mname,并且一名学生只在一个系学习,一个系只有一名系主任。

(1)关系模式R的主码是Sno,存在传递函数依赖SnoMname,这是非主属性对码的传递函数依赖,因此R属于2 NF,需要模式分解,提高其范式级别。

(2)对R进行分解,只有如下的ρ1和ρ2能够具有无损连接性,同时ρ1丢失了SdeptMname,不能保持函数依赖,但是ρ2保持函数依赖。两个分解ρ1和ρ2中的关系模式都已经达到了4 NF,不需要再进行分解。

ρ1={R1<{Sno,Sdept},{SnoSdept}>,R2<{Sno,Mname},{SnoMname}}

ρ2={R1<{Sno,Sdept},{SnoSdept}>,R2<{Sdept,Mname},{SdeptMname}}

(3)R的关系实例如下表1所示,ρ1和ρ2中的R1相同,它们的关系实例如下表2所示,ρ1和ρ2中的R2关系实例分别如下表3和表4所示,所有表中未加下划线的记录是原始内容,加下划线的记录是更新需要增加的记录。

考虑向数据库中添加(3,计算机系,王五)记录信息,模式分解之前的R需要检查函数依赖SdeptMname,发现“计算机系”的系主任姓名不惟一,因此拒绝该更新操作。

(4)模式分解后,先考虑ρ1,需要向R1和R2对应的关系实例表2和表3中分别添加(3,计算机系)和(3,王五),这两个操作没有违反各自的函数依赖集F1和F2,表面上看是可以允许的,但是考虑到分解丢失了函数依赖SdeptMname,不能够保持函数依赖,需要将表2和表3自然连接起来,形成表1,检查SdeptMname是否成立,结果发现是有问题的,从而拒绝该更新操作。如果不考虑后面的自然连接再判断的过程,则会将非法的数据加入到数据库中,这显然是不允许的。

(5)采用同样的方法,对分解ρ2进行分析,因为ρ2保持函数依赖,根据定理1,只要更新操作分别满足ρ2中的函数依赖集F1和F2,则更新操作一定是合法的。向R1和R2对应的关系实例表2和表4中分别添加(3,计算机系)和(计算机系,王五),检查后发现F2中的SdeptMname不成立,因此立即判断更新操作是非法的,不需要执行(4)中的自然连接再判断的过程,从而提高了数据更新的性能。同时,与(3)中原关系模式R的操作相比,数据更新性能并未受到影响。

经过以上分析可知,分解ρ1和ρ2能够具有无损连接性,但是ρ2还保持函数依赖,在数据更新时具有更好的性能,因此应该选择ρ2为最终的分解。在文献[2]中存在类似的实例分析,它指出在不保持函数依赖的分解ρ1中,当计算机系所有学生毕业后,相应的系主任信息也被删除,因此存在删除异常。当计算机系还没有学生时,相应的系主任信息无法加入,因此存在插入异常。这说明了在某种极端情况下不保持函数依赖带来的问题,并未体现在一般情况下不保持函数依赖会降低数据更新性能的后果。

4 结束语

关系数据库中的规范化理论在数据库设计中具重要作用。一个属于较低级别范式的关系模式可以通过模式分解得到若干较小的关系,减小了数据冗余度,使其达到更高一级的范式。关系模式分解的评判准则是无损连接性和保持函数依赖,分解必须具有无损连接性而尽量保持函数依赖。无损连接性能够保证不丢失数据库的信息,保持函数依赖能够数据更新性能不受到影响。规范化的过程中如果能够保持函数依赖,则在减小数据冗余度的同时不会影响数据更新的性能。如果不能保持函数依赖,则在减小数据冗余度的同时会降低数据更新的性能,此时数据库设计人员需要根据实际情况进行权衡。只有真正理解模式分解保持函数依赖的意义,才能设计出一个好的关系数据库。

摘要:规范化过程通过模式分解提高关系的范式级别,分解的评判准则是无损连接性和保持函数依赖。文中证明了保持函数依赖的一个定理,说明了模式分解保持函数依赖的意义是在减小数据冗余度的同时不会影响数据更新的性能,通过具体分析实例进一步说明了该结论。

关键词:规范化,模式分解,保持函数依赖

参考文献

[1]Date C J.深度探索关系数据库:实践者的关系理论[M].熊建国,译.北京:电子工业出版社,2007.

[2]王珊,萨师煊.数据库系统概论[M].4版.北京:高等教育出版社,2006.

[3]于思江,王小兵.极小函数依赖集的正确算法及其证明[J].西安邮电学院学报,2009,14(3):110-112.

[4]Neil P O,Neil E O.数据库原理、编程与性能[M].2版.周傲英,俞荣华,季文赟,等,译.北京:机械工业出版社,2002.

[5]Silberschatz A,Korth H F,Sudarshan S.数据库系统概念[M].5版.杨冬青,唐世渭,译.北京:机械工业出版社,2006.

三角函数分解 篇3

以求解线性方程组的矩阵三角分解 (LU分解) 为代表的数值计算在现代科学研究和工程技术中得到广泛应用。当今的计算机在高速的CPU和低速的存储器之间设置了二级或更多级的高速缓冲存储器 (简称高缓, cache) , 高缓的存取速度很快, 但容量较小, 将计算机当前正在运算的数据或经常需要访问的数据置于高缓中, 让数据存取的速度适应CPU的处理速度, 从而提高计算机整体的运算效率。对于大型方程组, 其系数矩阵很大, 标准LU分解直接对大矩阵进行计算, 一次参与运算的数据量很大, 不能很好地利用当今计算机的高缓, 运算效率不高。本文对LU分解运用分块算法, 将大矩阵的运算分解成分块子矩阵的运算, 使各部分小的分块矩阵的数据运算能够更多地在高缓中进行, 同时采取矩阵转置、循环展开等优化技术, 对分块算法中各个模块的计算设计了较高效的算法, 提高了算法的总体运算速度。

1 LU分解分块算法

LU分解用于求解线形方程组Ax=b。从高斯消去法出发, 其系数矩阵A可分解为单位下三角阵L和上三角阵U之积:

或:

A=LU (2)

利用分解后的LU, 通过回代可方便地求出方程组的解。由矩阵乘法可推导出:

式 (3) 即为标准LU分解的计算公式, 整个数据计算在大矩阵A的范围内进行。

下面介绍LU分解分块算法。将矩阵A分解成2×2块, 方程A=LU可写成:

由矩阵乘法有:

A11=L11U11 (5)

A12=L11U12 (6)

A21=L21U11 (7)

A22=L21U12+L22U22 (8)

LU分解分块算法的步骤如下:

① 对式 (5) 中的分块矩阵A11进行LU分解, 可得L11、U11。

② 将①中得到的L11求逆, 并代入式 (6) , 有U12 = L11 -1A12 。

③ 将①中得到的U11求逆, 并代入式 (7) , 有L21 = A21 U11 -1。

④ 改写式 (8) , 设A′22=A22-L21U12, 则式 (8) 成为A′22=L22U22, 即为A′22的LU分解。

对A′22的LU分解可继续采用上述分块过程进行计算。这样, 只要每次选取适当的分块大小, 分块后的各子矩阵间的运算就有可能在高缓中进行, 从而提高运算速度。

2 LU分解分块算法的设计与实现

综上所述, LU分解分块算法的实现涉及六个模块的计算, 分别是分块矩阵的LU分解、单位下三角阵求逆、上三角阵求逆、单位下三角阵乘矩形阵、矩形阵乘上三角阵和矩形阵乘矩形阵。

2.1 分块矩阵LU分解

分块矩阵 (A11) 的LU分解计算采用式 (3) 所示的标准LU分解公式。与标准LU分解算法不同之处在于: (1) 形参的个数不同。对于分块分解, 需提供分块矩阵入口地址p, 阶数d以及在原矩阵内的偏移量displace和原矩阵的阶数n。 (2) 选主元的范围不同。标准LU分解是在大矩阵A的范围内选主元, 而块内LU分解是在分块矩阵的范围内选主元。

2.2 单位下三角阵求逆

可证明单位下三角阵的逆矩阵仍为单位下三角阵, 上三角阵的逆矩阵仍为上三角阵且对角线上的元素分别变为原来的倒数。

因为L11为单位下三角矩阵, 求逆后对角线上元素仍为1, 所以在计算时, 只需对对角线以下部分求逆。设单位下三角阵为L= (lij) d×d, 其逆矩阵为B= (bij) d×d。由矩阵乘法k=ji (lik×bij) =0 (j<i) , 有bij=-k=ji-1 (lik×bkj) (i=2, 3, , d; j=2, 3, …, i-1) , 对角线上元素为bii=lii=1。在该算法中可直接将算出的bij保存在原lij的位置上。

2.3 上三角阵求逆

设上三角矩阵为U= (uij) d×d, 其逆矩阵为C= (cij) d×d。同样有cii=1uiik=ij (uik×ckj) =0, 得到cij=-1uiik=i+1j (uik×ckj) (i=d, d-1, , 1; j=d, d-1, …, i+2) 。计算出的cij直接保存在原uij的位置上。

上面分别设计了针对单位下三角阵和上三角阵的求逆算法, 而没有采用常用的原地全选主元高斯约旦消元法, 这样做可以充分利用三角矩阵的特点, 避免无效的计算。且高斯约旦消元法对矩阵求逆时, 要多次全选主元, 进行很多次比较, 新设计的两个算法中不需要进行比较。

2.4 矩阵相乘

矩阵相乘是LU分解分块算法中的核心运算, 提高矩阵相乘的运算速度是提高整个分块算法效率的关键。本文采用了矩阵分块、矩阵转置和循环展开等方法来提高矩阵相乘的运算速度。

矩阵分块 将大矩阵逐次一分为二, 直到划分到分块矩阵小于指定的大小。把这些分块矩阵当作元素来处理, 分块矩阵内部用普通三重循环来计算。设两相乘矩阵为AB, 则过程如下所示:

(A11A12A21A22) (B11B12B21B22) = (A11B11+A12B21A11B12+A12B22A21B11+A22B21A21B12+A22B22)

测试表明, 采用矩阵分块相乘技术可使运算速度得到较大提高。但对矩阵相乘仅仅采用分块还不够, 有必要引入矩阵转置和循环展开技术。

矩阵转置 对于矩阵Am×pBp×n=Cm×n, 无论分不分块, 最终都要用三重循环实现。对C语言来说, 是采取按行存放元素, 最内层循环内AC是顺序存取, 而B是不连续的。当矩阵的阶数较大时, 访问B中上下相邻两个元素跨度很大, cache 内可能放不下, 访问时不能命中, 仍要访问内存。这样的访问会进行m次, 当m很大时, 会降低效率。所以将B转置, 在最内层循环中用A的一行乘B的一行, 这样ABC都是顺序存取, 提高了数据访问速度。

循环展开 循环展开可减少测试循环终止条件的次数, 它不仅可以减少循环的开销, 还可以在两个方面提高程序的性能。一是使程序利用多个执行部件并行计算, 充分发挥现代计算机多执行部件的特性。二是循环展开有助于重用装入寄存器中的数据, 降低了寄存器与cache以及cache 与内存之间对通信带宽的需求。展开多层循环的外层循环可以使一些变量得到重用。

下面介绍三处矩阵相乘的具体实现。

2.4.1 单位下三角阵乘矩形阵

计算U12 = L11-1A12 , 先将A12转置, 然后相乘, 结果保存于一临时矩阵T中。计算时需使用一大小为m的额外向量空间, 来保存A12转置后矩阵当前循环到的一行。乘运算从右边矩阵中的第一行开始, 逐行向下进行。采用循环展开技术, 第一层和第二层循环各展开1次。此时需要考虑两矩阵的行数为奇数的情况。当左边矩阵行数为奇数时, 计算结果存入右边矩阵的最后一列中;当右边矩阵行数为奇数时, 计算结果存入右边矩阵的最后一行中。

由于L11-1的对角线上元素全为1, 所以它们不需和T中的对应元素相乘, 可直接将这些元素作为最内层循环乘积项累加的变量初始值, 这样可减少m×n次乘操作。又L11-1为下三角阵, 只需将其对角线以下元素与T中对应元素相乘即可。

2.4.2 矩形阵乘上三角阵

在计算L21 = A21 U11-1时, 将结果保存在左边矩阵 (A21) 中, 使用一大小为n的额外向量空间来保存左边矩阵当前循环到的一行, 从左边矩阵的第一行开始, 逐行向下运算。计算时, 第一层和第二层各循环展开1次。同样要考虑两矩阵的行数或列数为奇数的情况。当左边矩阵行数为奇数时, 计算结果存入左边矩阵最后一行;右边矩阵列数为奇数时, 结果存入左边矩阵最后一列。

在该模块中, 矩形阵A21的列数大于等于上三角阵U1l的列数, 且U11对角线以下元素全为 0, 所以没有对其转置, 只需将三角阵对角线以上元素与左边矩形阵中对应元素相乘即可。

2.4.3 矩形阵乘矩形阵

A′22=A22-L21U12中, L21U12为两矩形阵相乘。计算时, 对L21进行分块, 将U12转置并分块。由2.4.1节中得到的结果T可直接作为L21所乘的对象, 并在做乘法时将所得的积取反, 相当于完成减运算, 这些操作有利于提高整个算法的效率。

L21、T分块时, 取m=n=d, 即分块矩阵为d阶方阵。这样一直分块下去, 直至L21最后所剩一块矩阵的行数为N%d+d, 列数为d;T最后所剩一块矩阵的行数为d, 列数为N%d+d, 其中N为原矩阵A的阶数。

分块子矩阵相乘时采用循环展开技术, 第一层和第二层循环各展开1次。此时由于参与计算的分块矩阵均为稠密型矩阵, 可直接展开不必象有三角阵参与那样考虑边界处的处理, 计算结果最终保存在L21分块对应的位置上。

3 测试与分析

3.1 测试结果与比较

两种算法在PC机上测试, CPU为Pentium 4, 主频2.8GHz, 内存512MB, 编译系统采用VC++6.0。在验证算法计算正确的基础上, 依次取矩阵的阶数为300~6000, 对两种算法进行运行和比较。对于分块算法, 取每次分块大小d=150。测得的算法运行时间如表1所示, 表中第四行为分块算法比标准算法效率提高的百分比。

表1 (续)

由表1可见, LU分解分块算法明显优于标准算法, 在大矩阵情况下, 运算速度提高50%以上。

3.2 测试结果分析

标准LU分解算法和LU分解分块算法中的运算时间主要用在乘操作和除操作上, 可分别计算出两种算法中乘操作和除操作的运算次数, 限于篇幅在此只给出结果, 作简要分析。

分块算法 (取每次分块大小为d) 的乘运算量约为n3/3+ (n%d-0.5) n2, 除运算量约为 (d+1) n/2。标准算法的乘运算量为n (2n-1) (n-1) /6, 除运算量为n (n-1) /2。当矩阵阶数n很大时, 两算法的乘运算量相差不大, 而分块算法除运算量比标准算法降低了一阶。计算机做除法最耗时, 且标准算法没能充分利用高缓。分块算法在较好地利用高缓的同时, 通过矩阵分块、矩阵转置和循环展开等技术加速了矩阵乘法的运算, 所以运算效率显著提高。

LU分解分块算法比标准LU分解算法多用约d×n的存储空间。

4 结束语

本文将分块算法引入LU分解, 并加以有效实现, 使算法的运算速度得到较大提高。由算法分析和测试可知, 在分块算法中, 两矩形阵相乘的用时最多, 约占总时间的70%。因此, 要想进一步提高运算速度, 必须提高两矩形阵相乘的效率。可以对矩阵相乘作并行处理, 事实上在LU分解分块算法中, 有些模块是可以并行计算的。再者就是研究出达到最佳运算效率时总阶数n与分块大小d之间、cache大小与d之间的关系。以上这些问题将是以后进一步研究的目标。

摘要:对稠密型线性方程组的系数矩阵进行分块LU分解, 更充分地利用高速缓存, 提高运算效率。对LU分解分块算法进行了研究, 用VC++6.0对分块算法进行实现, 并与标准的LU分解算法进行比较。在大矩阵情况下, 分块算法比标准算法运算速度提高50%以上。

关键词:LU分解,矩阵分块,矩阵快速相乘,VC++6.0

参考文献

[1]陈建平.LU分解递归算法的研究[J].计算机科学, 2004, 31 (6) :141-142.

[2]李玉成, 朱鹏.BLAS的加速方法与实现技术[J].数值计算与计算机应用, 1998, 19 (3) :227-240.

[3]李玉成.LAPACK中的分块算法及其效果[J].数值计算与计算机应用, 2001, 22 (3) :172-180.

[4]李忠泽, 陈瑾, 龙翔, 等.基于PentiumPro的高性能BLAS的设计与实现[J].北京航空航天大学学报, 1998, 24 (4) :454-457.

[5]王小牛, 冯百明.基于存储的矩阵乘积优化算法[J].西北师范大学学报, 2005, 41 (1) :22-24.

[6]刘萍.数值计算方法[M].2版.北京:人民邮电出版社, 2007:85-88.

[7]Elmroth E, Gustavson F, Jonsson I, et al.Recursive Blocked Algorithms and Hybrid Data Structures.SIAMReview, 2004, 46 (1) :3-45.

[8]Toledo S.Locality of Reference in LUDecomposition with partial pivo-ting[J].SIAM Journal of Matrix Analysis&Application, 1997, 18 (4) :1065-1081.

三角函数分解 篇4

随着电力电子技术的日益发展,大量非线性负荷波动,各种变频调速装置、电力电子装置在电力系统中广泛应用,使得电力系统中电压、电流波形发生畸变,导致电力系统的谐波问题日益严重,同时非整数次谐波——间谐波和次谐波,也引起了国内外学者的广泛关注[1]。谐波和间谐波在电力系统中的大量存在导致了用户侧电流的畸变和电流波动频繁,因此需要可以同时准确分析谐波和间谐波频谱的方法[2]。

目前常见的分析方法有加窗插值FFT算法[3]、小波算法[4,5]、支持向量机算法[6]、Prony方法[7]、多信号分类法(MUSIC)[8,9]等,各种方法都存在优点和不足,特别在低信噪比时检测精度受到较大影响。

本文利用信号自相关矩阵的特征值分解理论,将信号的自相关矩阵分解为信号子空间和噪声子空间,利用2个子空间的正交性,进一步分解噪声子空间,对其进行变换构造出基于噪声子空间分解的特征多项式(DNS-MUSIC函数),求解此多项式后得到信号基波和谐波频率预估计,结合消去法[10]逐步得到频率的精确估计,并利用扩展Prony方法估计得到信号的幅值和相位。

实例仿真实验与其他方法比较表明,所提方法具有较高的频率分辨率,而且对数据长度没有要求,还具有一定的抗噪声能力。

1 数学模型

电力系统的电压(或电流)信号为

其中,x(t)为电压或电流信号;p为所含的谐波和间谐波的个数;u(t)为白噪声项;Ai为第i次谐波的幅值;fi为第i次谐波的频率;φi为第i次谐波的初始相位。

特征多项式是基于复信号计算的,可以借助Hilbert变换将原信号相移90°,将式(1)转换成复频率表达式:

2 空间谱频率估计的理论基础

本文假设噪声向量的每个元素都是零均值的复白噪声,相互独立,并且具有相同的方差σ2。由自相关矩阵的定义知:

利用特征值分解的理论(详细理论可参见文献[11]),自相关矩阵可分解为

自相关矩阵R的特征值为

这表明,当存在加性观测白噪声时,很容易将自相关矩阵R的前2p个特征值与后面的M-2p个特征值区分开。因此称前2p个特征值为信号的主特征值,其余的M-2p个特征值为噪声特征值。根据信号特征值和噪声特征值,又可以将特征矩阵U的列向量分成2个部分,分别由信号特征向量和噪声特征向量组成。

3 构造DNS-MUSIC函数

经典MUSIC法要在频率轴上进行全域搜索,计算量较大,而且频率分辨率较低,耗时较长,在信噪比较低的情况下可能产生虚假频谱[11]。为了改善经典MUSIC法的不足,有学者提出了求根MUSIC法,通过研究噪声空间和信号空间的正交性,构造MUSIC型函数[12]:

该多项式的阶数为2(M-1),通过求解式(6),可得到(M-1)对根,然后选出位于接近单位圆上的2p个根,估计出信号的频率。求根MUSIC法只采用单位圆上的根的信息来估计频率,造成了非单位圆上根的信息的丢失,在强噪声环境中对单位圆上根的误判可能产生虚假频谱,在修正过程中用P T(z-1)代替PH(z),而实际过程中PT(z-1)=PH(z)未必成立也造成了估计偏差。上述因素使得求根MUSIC法在估计电力系统信号频谱时检测精度受到了影响。鉴于此,本文提出了一种新的方法,通过噪声子空间分解DNS-MUSIC函数来求解频谱信息。

由文献[12]知定义多项式:

其中,a(z)=[1 z…z M-1]T;ei是自相关矩阵R的第i个特征向量。

由于信号子空间与噪声子空间正交,当zi=exp(jωi)时,即多项式的根正好位于单位圆上时,a(exp(jω))是一个空间频率ω的导向矢量。由特征结构类算法可知,a(exp(jω))就是信号的导向矢量,则zi=exp(jωi)(i=1,…,2 p)应该是式(7)中M-2 p个多项式的2p个根,即

从式(8)可知,这些多项式应该有一个能够用F(z)表示的2 p次的最大公因式,所有信号的频率可以通过求解F(z)来获取。F(z)的求解方法如下。

假设存在(M-2p)×1维向量:

通过式(9),对噪声子空间的特征向量进行变换,即

其中,VM为M×(M-2p)阶矩阵;VM1为2p×(M-2p)阶子矩阵;VM 2为(M-2 p)×(M-2 p)阶子矩阵;M为阵元个数,p为信号个数。

显然,式(10)是M-2 p个噪声空间的噪声向量的线性组合,噪声子空间与信号子空间正交,所以它的线性组合与信号空间也是正交的。因此,这种变换不会改变噪声特征向量的参数信息,而是把噪声子空间的参数信息最大限度地集中到VM1子矩阵上,这样便于利用VM1中的噪声特征向量来提取空间的参数信息。

由式(10)知:

由式(12)知,通过对噪声空间的特征向量进行变换,将噪声空间的参数信息集中到向量a中,防止空间参数信息的泄漏和丢失。据此可以构造出DNS-MUSIC函数:

该多项式的阶数为2p,即有p对共轭根,并且正好分布在单位圆上,与式(6)相比降低了特征多项式的阶数,简化了特征多项式的求解运算过程,减少了有效根解的影响因素,求解的结果将更真实。

求式(13)的根得:

所以有

DNS-MUSIC算法得到2p阶的特征多项式,只有p对共轭根(p个频率成分),即使信号中混有噪声,数据矩阵存在误差时,也不会产生虚假频率成分,即不会产生伪频谱。

DNS-MUSIC算法对噪声子空间的特征向量进行线性变换,可直接得到特征多项式的系数,不需要对其修正处理,防止了求根MUSIC法在修正过程中用PT(z-1)代替PH(z)造成的误差,也可以在一定程度上提高检测的精度。

DNS-MUSIC函数充分利用了噪声特征向量提取空间参数信息,有效克服了求根MUSIC法从很多个多项式中求解少量单位圆上根时带来的较大误差,也防止了在强噪声环境中求根MUSIC法会对单位圆上根的误判而产生伪谱的情况。因此该方法能大幅提高频率分辨率,并且使检测结果更加稳定准确。

综上分析,通过构造DNS-MUSIC函数,可以比求根MUSIC法更准确地估计信号频率,而且不会丢失信号信息或产生虚假信号信息。

4 谐波和间谐波频率检测过程

电力系统信号中谐波、间谐波分量的幅值一般仅为基波分量幅值的百分之几或更小。当对其进行非同步采样时,大幅值频率分量的频谱泄漏有可能淹没小幅值的频率分量[3]。

本文采用文献[12]中提到的递减处理思想作为本文检测过程的处理方法,其主要思想是:先检测原信号中的大信号,为了消除大信号的频谱泄漏对小信号的影响,在原信号的基础上减去大信号然后再依次检测小信号。这样不会丢失原信号的信息,并可以提高检测的精度。

设电力系统的信号为x(t),则

其中,νh(t)为大幅值信号;νi(t)为小幅值信号;

对电力系统信号按照奈奎斯特抽样定理抽样并把式(16)写成向量形式:

其中,W=[a b]T为待估系数;e为误差向量;N为采样点数。

系数W=[a b]T采用最小二乘法进行预估计处理,则有

由于式(8)的频率已由DNS-MUSIC算法确定,则矩阵R已知,向量X为采样信号,则通过式(18)可方便求出W,即可确定νh(t)的信息,然后由式(16)便可知:

由于在时域中减掉大幅值信号,从而减弱或消除了大幅值信号所带来的频谱泄漏,再继续采用DNS-MUSIC算法对式(19)定义的信号重采样进行频谱分析,即可准确地检测出小幅值信号的频率。

通过这种方法在电力系统信号中顺次减掉大幅值的基波、谐波和大幅值间谐波成分,从而检测出更小幅值的间谐波分量,有效抑制了基波、谐波和间谐波之间的干扰,进一步提高频率检测的精度和频率分辨率,检测过程如图1所示。

5 谐波和间谐波幅值和相位检测

通过以上过程得到谐波和间谐波的频率和个数,下面采取扩展Prony法来获取信号的幅值和相位[13]。扩展Prony法采用2p个具有任意幅值、相位、频率与衰减因子的指数函数为数学模型,式(2)可表示为

令阻尼因子α=0,并且用(采样近似值)代替x(n)(实际值),可知DNS-MUSIC函数求得的根即为Prony的极点。利用DNS-MUSIC算法求出信号中谐波和间谐波的频率fi之后,p和zi就成为已知量(对于实正弦信号,fi还应包括与正频率对应的负频率,频率个数2p应为正频率个数的2倍)。

于是式(16)的最小二乘解(详细原理可参见文献[13])为

求出bi后根据式(16)可以得到:

对于实信号,信号的实际幅值是式(22)计算幅值的2倍。

6 仿真实例和性能分析

根据电力系统信号(未加噪声)中谐波/间谐波的特点,设待检测信号为

信号基波频率50 Hz,含有100 Hz和200 Hz的谐波,根据间谐波特性分别设置40 Hz、125 Hz和201 Hz的间谐波成分,采样频率为1 000 Hz,阵元个数为40,采样点数为1000。

为了验证本文算法的性能,分别在无噪声和信噪比为10 d B和15 d B的情况下,用MATLAB7.4进行了插值FFT算法、求根MUSIC法与本文所提出的DNS-MUSIC算法进行仿真比较,结果如图2—4所示。

在信噪比为10 d B情况下,插值FFT算法无法检测到小幅值的间谐波成分,如图2(a)所示,只有当信噪比达到15d B才能检测到间谐波频率,如图2(b)所示,而且存在较大的误差,详见表1。这说明插值FFT算法中小信号对噪声比较敏感,其频率分辨率为fs/N,在数据量较少的情况下其频率分辨率也比较低[3]。

图3为求根MUSIC法的功率谱估计图。可以看出,求根MUSIC法在一定程度上提高了频率分辨率,并具有一定的抗噪性,但是在信噪比较低的情况下,基波的频谱泄漏较为严重,致使40 Hz的间谐波被基波频谱噪声完全淹没,而且出现了270 Hz的伪谱,对电力系统信号的频率造成误判,如图3(a)所示。

随着信噪比的提高,当信噪比达到15 d B时,40 Hz的间谐波成分被检测出来,而且抑制了伪谱成分,如图3(b)所示,但是在频率检测精度上有一定的偏差,详见表1。可见,插值FFT算法和求根MUSIC法在低信噪比的情况下不能准确检测电力系统信号的频率成分。

图4为DNS-MUSIC算法的功率谱估计图。可以看出,在低信噪比情况下,DNS-MUSIC算法能够检测出电力系统信号的频率个数和频率值,如图4(a)所示,克服了低信噪比情况下分辨率低、抗噪性差的不足。随着信噪比的提高,DNS-MUSIC算法在精度上大幅提高。在相同的仿真条件下,DNS-MUSIC算法在分辨率和精度上都远高于前2种算法,详见表1。

由表1可以看出,在无噪声的情况下,插值FFT算法、求根MUSIC法和DNS-MUSIC算法都可以检测出电力系统信号的频率、幅值和相位。但插值FFT算法检测出的电力系统信号成分存在较大误差,不能达到电力系统信号的高精度的检测要求。求根MUSIC法可以准确检测出基波和谐波的成分的相关参数,但是对间谐波的检测存在一定的误差。DNS-MUSIC算法在无噪声的仿真条件下可以准确检测出电力系统信号的各个成分。此外还可以看出,插值FFT算法和求根MUSIC法在低信噪比的情况下有可能淹没小幅值的间谐波成分,即便信噪比足够高,谐波和间谐波的估计精度也不高,而DNS-MUSIC算法都能准确分辨出谐波和间谐波成分,克服了求根MUSIC法在低信噪比情况下产生虚假频谱的不足,而且当信噪比相同的情况下DNS-MUSIC算法估计谐波间谐波频率的精度远高于插值FFT算法和求根MUSIC法,在频率达到较高精度的前提下Prony法可以准确估计出电力系统信号的相位和幅值。表1中还给出了通过Prony算法得到的幅值和相位的检测结果。

下面采用均方根误差(RMSE)来验证本文算法检测信号频率的准确性。

其中,Ns为仿真实验次数;fi为第i次实验获得的频率估计;f为频率的真实值。

利用3种检测方法,在不同信噪比下,对频率为200 Hz的信号分量分别做20次仿真实验,得到均方根误差曲线如图5所示,表明DNS-MUSIC算法在低信噪比情况下均方误差明显低于其他2种算法,说明该算法在低信噪比情况下具有良好的稳定性。随着信噪比的增加,均方根误差改善明显,在相同的仿真条件下,本文提出的算法大幅提高了频率分辨率和抗噪声的能力,而且具有较好的稳定性和准确性。

7 结论

本文在传统MUSIC的基础上,通过分析噪声空间的特性,构造出DNS-MUSIC函数,结合递减消噪法,应用到电力系统信号的检测过程中,得到电力系统信号的频率成分的估计值,通过Prony法,检测出电力系统信号的相位和幅值,取得了良好的频谱效果,尤其在间谐波成分的检测中取得了更高的精度。该方法克服了插值FFT算法在短数据情况下分辨率低的不足,弥补了求根MUSIC法在低信噪比情况下出现虚假频谱的不足。该方法比求根MUSIC法具有更高的精度,而且计算量小。但是由于Prony法对噪声比较敏感,采用Prony法估计电力系统信号的相位和幅值有一定的误差,有待进一步改进。综合而言,本文提出的新算法具有稳定性好、分辨率高、抗噪性强、运算时间短的特点,对电力系统信号的实时精确检测具有一定意义。

摘要:针对电力系统中存在的谐波和间谐波问题,提出了基于噪声子空间分解MUSIC(DNS-MUSIC)函数的谐波/间谐波检测方法。利用信号自相关矩阵的特征值分解理论,将信号的自相关矩阵分解为信号子空间和噪声子空间,利用2个子空间的正交性进一步分解噪声子空间,对其进行变换,构造出基于噪声子空间分解的特征多项式(DNS-MUSIC函数),求解该多项式得到信号基波和谐波频率预估计,结合消噪思想检测电力系统信号频率成分,然后利用扩展Prony法检测信号的幅值和相位。通过仿真实验与其他经典算法比较,结果证明了所提算法的可行性、高效性和稳定性。

三角函数分解 篇5

关键词:环境全要素生产率,SBM,方向性距离函数,GML指数

改革开放以来,我国经济增长取得了令人瞩目的成绩,1978~2009年的年均增长率在9%以上。在盘点所取得的辉煌业绩时,我们发现高速的经济增长是以能源的高投入、高消耗为特征,经济对能源的过度依赖导致大气中二氧化碳等温室气体浓度增加,诱发全球气候变暖及环境污染,从而降低了经济增长质量,使得经济增长速度大打折扣。我国正处于工业化、城市化、现代化快速发展阶段,重化工业发展迅速,大规模基础设施建设不可能停止,能源需求的快速增长一时难以改变,碳排放量在短时间内也不能得到遏制。能源耗竭和二氧化碳排放日益成为制约经济可持续发展的约束条件,为此,节能减排被提到前所未有的战略高度,我国政府在“十一五”发展规划中提出:到2020年单位GDP二氧化碳排放比2005年下降40%~50%,达到这一目标的关键是运用技术提高生产效率。因此,我们有必要深入分析资源环境约束下我国全要素生产率处于一个什么样的状态,通过分析环境全要素生产率增长的源泉,为缓解经济高速增长与能源耗竭、二氧化碳排放之间的矛盾提供一些政策启示。

1 文献综述

国内外学者已经注意到能源消费对产出的影响。Hu和Wang(2006)最早采用规模报酬不变(CRS)的DEA模型,将能源消费引入生产函数中,研究了中国省际间的全要素生产率[1]。由于经济增长是一个伴随着非期望产出(环境污染)不断产生的期望产出增加过程,我国经济增长对能源消费的路径依赖,以及以煤为主的能源结构致使我国能源消费与碳排放量随着GDP的提高也在不断增加,而忽视二氧化碳等非期望产出的传统生产率测度思路会导致有偏差的生产率增长(Chung,1997)[2]。袁晓玲、张宝山等(2009)利用Malmquist指数测算出了包含非期望产出环境污染的中国省际全要素能源效率,结果表明,我国总体能源效率偏低,在维持一定水平后呈现缓慢下降的趋势,区域差异较为显著[3]。但以上文献的测度思路是把污染排放指数的倒数作为产出变量带入模型求解,这与将坏产出作为投入变量一样,是将好产出和坏产出进行非对称处理,扭曲了对经济绩效和社会福利水平的评价,从而误导政策建议(Hailu,Veeman,2001)[4],并且是在采用保持坏的产出不变的情况下利用基于产出距离函数的Malmquist指数计算的,得出的结果存在一定的偏差。Chambers,Chung和Fare(1996)归纳了Luenberger(1992)利润函数,提出了非径向的新思路——方向性距离函数[5]。该函数是Shephard产出距离函数的一般化,可以处理投入产出同时变化的情况,并在此基础上构建了可用来处理非期望产出的Malmquist-Luenberger(ML)生产率指数。方向性距离函数和ML指数在测算及分解全要素生产率时能够从经济增长与二氧化碳等非期望产出减少2个方面考虑,从而分析节能减排效果。面对保增长和减排放的双重约束,沈能(2012)测算了加入环境约束因素后我国工业各行业的环境效率,结果表明,若考虑非期望产出的影响,中国工业环境效率有明显的下降[6]。Fare(2001)、Kumar(2006)以及胡鞍钢(2008)利用传统的M生产率指数和ML生产率指数测度发达国家和发展中国家的全要素生产率,研究结果表明,考虑环境因素的全要素生产率年均增长率与忽略环境因素的全要素生产率的年均增长率之间存在显著的差异[7,8]。陈诗一(2009)同样采用方向性距离函数和ML生产率指数模拟了中国工业节能减排的损失和收益[9]。孙传旺等(2010)基于ML生产率指数,测度并分解了低碳经济发展中我国的全要素生产率,发现技术进步是推动碳强度约束下全要素生产率提高的主要因素[10]。很多文献都是运用径向的、角度的DEA来计算方向性距离函数,当存在投入或产出过度,即存在投入或产出的非零松弛量时,径向的DEA方法会高估效率;而角度的DEA方法忽视了投入或产出的某一个方面,计算的结果并不准确。非径向、非角度的基于松弛变量SBM(Slack-Based Measure)方向性距离函数可以克服上述2个缺陷。王兵等(2010)运用SBM方向性距离函数和ML生产率指标测度了环境约束下中国区域全要素生产率的增长情况,认为能源的过多使用以及污染物的过度排放是环境无效率的主要来源[11]。庞瑞芝等(2011)同样采用SBM方向性距离函数测算了我国各省份工业部门的“新型工业化”生产力[12]。

本文将能源消费量作为生产投入之一,同时将产生的二氧化碳作为非期望产出,测度及分解了碳排放约束下中国区域全要素生产率。以上文献中测算环境全要素生产率时所用的ML生产率指数一般采用2个当期ML指数几何平均的形式,导致ML指数在测度跨期方向性距离函数时,如果t+1期的投入产出值在t期的技术下是不可行的,则面临着一个潜在的线性规划无解的问题。此外,以几何平均形式表示的ML指数不具有循环性或传递性(Circularity or Transitivity)Oh(2010)将Global Malmquist生产率概念和方向性距离函数相结合,构建了ML指数的替代方法——Global Malmquist- Luenberger(GML)指数[13]。该方法不仅可以应付多输出和多输入,以及生产带来的环境污染非期望产出问题,而且避免了传统ML指数的线性规划无解的问题,同时是可循环累加的。本文尝试将SBM方向性距离函数和GML指数运用到碳排放约束下中国全要素生产率累积值的测度中,并将全要素生产率进一步分解为纯技术效率变动、纯技术进步、规模效率变动和技术规模变动四部分,以便从能源节约、二氧化碳减少以及GDP增长3个方面分析中国经济增长质量的长期变动趋势,考察增长节能减排三目标的实现情况。

2 环境全要素生产率的测算方法

本文首先构建一个既包括期望产出(国内生产总值)又包括非期望产出(二氧化碳)的生产可能性集,即环境技术。假设每一个决策单元使用N种投入x=(x1,x2,Λ,xN)∈RN+,得到M种期望产出y=(y1,y2,Λ,yM)∈RM+,以及J种非期望产出b=(b1,b2,Λ,bJ)∈RJ+。利用第K个省份第t期的投入和产出值(xk,t,yk,t,bk,t),可以构造生产可行性集Pt(xt)∶Pt(xt)={(yt,bt):xt能生产(yt,bt)},xt∈RN+,t=1,2,Λ,T。Pt(xt)满足闭集和有界集、投入与期望产出强可处置性、非期望产出弱处置公理、期望产出和非期望产出零结合公理等条件。

Pt(xt)是依赖t时刻当期的生产技术构造的参考生产集,集合中的每一个数据仅仅是那个时刻得到的观测值。运用当期数据来确定当期的生产前沿有可能导致技术倒退结论的出现,运用全局生产技术集(Global Production Technology Set),即每一年的参考技术由当期及其以前所有可能得到的投入产出数据确定。Oh(2010)将全局生产技术集定义为所有当期生产技术集的并集,即PG(x)=P1(x1)∪P2(x2)∪Λ∪PT(xT)。运用整个时间段内贯穿整个生产集的观测数据设置的全局生产技术构建了一个单一的参考生产前沿,所有的决策单元都以此为“标杆”,增强了决策单元技术效率之间的可比性。

2.1 方向性距离函数

很多学者已经证实,方向性距离函数能较好地解决包含“坏”产出的效率评价问题,本文同样利用方向性距离函数计算生产可行性集的最优解,以此反映经济活动中对环境技术的使用情况。给定投入组合后,构造的方向性距离函数如下所示:

Do(x,y,b;g)=max{β∶(y,b)+βg∈P(x)} (1)

其中,g=(y,-b)为产出水平扩张的方向向量,β为方向性距离函数值。生产可行性集与方向性距离函数的关系如图1所示:

根据Fukuyama & Weber(2009),资源环境约束下的SBM方向性距离函数的表示形式如下所示[14]:

Di(x,y,b;g)=max1Νn=1Νsnxgnx+1Ν(m=1Μsmygmy+j=1Jsjbgjb)2

s.t. t=1Τk=1Κzksykms-smy=yimtt=1Τk=1Κzksbkjs+sjb=bijtt=1Τk=1ΚzksXkns+snx=xint,

zkt≥0,snx≥0,smy≥0,sjb≥0

其中,g=(-gx,gy,-gb)为代表投入压缩、好产出扩张和坏产出压缩的方向性向量,而(snx,smy,sjb)则代表投入、好产出和坏产出的松弛变量,zkt表示每一个横截面观测值的权重。以设定的方向性向量为权数,同时寻求投入(x)最小化、“好”产出(y)最大化以及“坏”产出(b)最小化,同时兼顾环境规制和经济增长,实现经济“又快又好”地发展。

2.2 Global Malmquist-Luenberger指数模型

由于采用几何平均形式的ML指数测算出的全要素生产率不具有循环累乘性,只能进行相邻期间生产效率的短期变动分析,无法观察生产效率的长期增长趋势。而且其中混合方向性距离函数容易导致线性规划无可行解,图1中的A点就是全要素生产率无解的最好证明。而GML指数是基于全局生产技术集构造的,能有效避免线性规划无解的缺陷。同时,这种连续生产前沿面避免了生产前沿向内偏移的可能性,也就是说,可以避免“技术倒退”现象出现的可能性,从而也就避免了生产效率的“被动”提高。

全要素生产率的分解一般遵循Fare等(1994)和Ray & Desli(1997)提出的方法,虽然都将生产率的增长分解为:技术进步、纯效率变化和规模效率变化,但Fare等是在规模报酬不变的假设下分解,Ray & Desli是在规模报酬可变的假设下分解[15,16]。Grifell & Lovell(1999)认为虽然Fare的分解并不准确,但对生产率增长的测度是准确的,而Ray & Desli的分解思路准确,但其生产率增长的测度并不准确[17]。本文将两者的分解方法结合,在不变规模报酬下将 分解为纯效率变化(GPEC)、纯技术进步(GPTC)、技术规模变化(GSTC)和规模效率变化(GSEC)4个因子,如下所示:

GMLtt+1=GPECtt+1×GPTCtt+1×GSTCtt+1×GSECtt+1

其中,D0G(xτ,yτ,bτ;gτ)=max{β∶(yτ,bτ)+βgτ∈PG(xτ)}为全局方向性距离函数,记S(x,y,b;g)=1+D(xyb;g)。BΡGGτ=SVG(xτyτbτ;yτ,-bτ)/SVt(xτ,yτ,bτ;yτ,-bτ)是τ时刻沿着方向向量(yτ,-bτ)标示的射线方向,当期技术前沿和全局技术前沿之间的最佳实践差距,于是,表示t和t+1时期之间最佳实践差距BPGG,τ变动的指数GPTCtt+1,可用来测度在更多期望产出和更少非期望产出的方向上,当期技术前沿向全局技术前沿移动的紧密程度,为不同时期的技术变动情况提供了一种新的测度形式。而GPECtt+1、GSECtt+1和TSTCtt+1分别表示纯技术效率变动、规模效率变动和技术规模变动。如果GPECtt+1>(<)1,则表明在t和t+1时期之间,决策单元向着前沿面移动(远离前沿面),从而使生产活动的效率改善(效率恶化)。GPTCtt+1>(<)1表示相对于t期的生产技术,t+1期的生产技术更加接近(远离)全局生产技术,反映了技术进步(技术倒退)。GSECtt+1表示规模效率变化,是由效率值变化引发的规模效应,GSECtt+1>(<)1表示规模效率提高(下降)。GSTCtt+1测度技术进步的规模效应,GSTCtt+1>(<)1表示技术偏离不变规模报酬(向不变规模报酬移动)。全要素生产率变动指数表示决策单位向最佳生产实践前沿面的移动。若碳排放约束下全要素生产率得到提升则GMLtt+1>1,反之,若全要素生产率下降则GMLtt+1<1。技术进步(技术倒退)、效率改善(效率恶化)技术规模提高(下降)和规模效率提高(下降)将对全要素生产率的提升起到促进作用(阻碍作用)。

为了测算与分解当期和全局Malmquist-Luenberger指数,需要借助线性规划方法在不变规模报酬和可变规模报酬下分别计算4个方向性距离函数。以t期为例,其当期方向性距离函数可通过求解如下线性规划得到:

D0t(xt,yt,bt;yt,-bt)=maxβ

s.t. k=1Κzktykmt(1+β)ymtm=1ΚΜk=1Κzktbkjt=(1-β)bjt, j=1,K,J (2)

k=1Κzktxknt(1-β)xnt, n=1,K,N

zkt≥0, k=1,K,K

全局方向性距离函数在构建生产可能性集时,需检测整个时间段内的生产技术,因此,t期的全局方向性距离函数求解模型如下所示:

D0G(xt,yt,bt;yt,-bt)=maxβ

s.t. t=1Τk=1Κzktykmt(1+β)ymtm=1ΚΜt=1Τk=1Κzktbkjt=(1-β)bjt, j=1,K,J (3)

t=1Τk=1Κzktxknt(1-β)xnt, n=1,K,N

zkt≥0, k=1,K,K

模型(2)或(3)中的最优值β*反映的是决策单元投入削减、期望产出增加和非期望产出减少的最大提升空间。本文将模型(2)或(3)测算得到的省际全要素生产率得分界定为σ=1-β1+β,该指标同时反映了节能和减排实施效果。

3 实证分析结果

本文以2001~2009年中国内地的29个省市(重庆市归于四川省,西藏因数据不全不在考察范围内)为决策单元,测算并分解各省市碳排放约束下的全要素生产率(本文将其称之为环境全要素生产率)。数据来源于历年的《中国统计年鉴》、《中国能源统计年鉴》和《中国环境统计年鉴》。

3.1 变量选取与数据说明

假定生产过程中需要3种投入要素:劳动力、物质资本存量、能源消费,生产出期望产出GDP和非期望产出。

(1)劳动力。

本文使用年末社会从业人员总量作为生产中劳动力的投入。

(2)物质资本存量。

本文采用永续盘存法:Ki,t=Ii,t+(1-δi,t)Ki,t-1估计每年的资本存量,其中Ki,t是地区i第t年的资本存量,Ii,t是地区i第t年的投资,δi,t是地区i第t年的固定资产折旧。本文采用张军(2004)的相关数据,并根据其方法将时间延长到2009年[18]。为了研究的可比性,本文将各省市历年的资本存量按照2000年的可比价格进行了折算。

(3)能源消费。

用各省每年的能源消耗量表示能源投入,由于各省的能源消费种类不一,所以统计上把煤炭、石油、天然气和水电等4种主要一次性能源的消费量转换成统一单位加总而得。

期望产出用各省每年的GDP表示期望产出,并采用GDP平减指数将其换算为2000年为基期。由于二氧化碳排放是诱发全球变暖的主要因素,是环境规制中的主要控制对象,且二氧化碳排放主要是由能源消费产生的。因此,以二氧化碳排放量作为非期望产出指标。能源消费CO2排放的计算公式为:CΟ2=i=17CΟ2i=i=17Ei×ΝCVi×CCi×COFi*(44/12)。

其中,CO2,i代表第i种能源消费的CO2排放量,本文考虑了煤炭、焦炭、汽油、煤油、柴油、燃料油和天然气共7种能源消费种类(i=7);Ei表示第i种能源的消费量(t,m3);NCVi为第i种能源的平均低位发热值(KJ/Kg,m3);CCi为第i种能源的碳含量(tc/TJ),COFi是第i种能源的碳氧化率。各种能源的CO2排放系数见表1。由于生产主要利用电能作为动力,本文还根据电力消耗与二氧化碳的折算系数估算我国生产过程中产生的二氧化碳。根据国家发改委2008年公布的中国区域电网的基准线排放因子,华东电网的电量边际排放因子为0.954tCO2/MWh,本文以此作为二氧化碳排放因子,对消耗电力产生的间接二氧化碳排放也进行估计。将化石能源燃烧产生的直接二氧化碳和电力消耗产生的间接二氧化碳加总,得到总的二氧化碳排放量。

注:数据来源于IPCC(2006)及国家发改委能源研究所(2007)

3.2 考虑非期望产出的中国各省市环境效率分析

运用中国各省市2001~2009年间的投入产出数据,利用上述线性规划,在可变规模报酬以及投入产出双导向的情况下,得出历年各省市的全局环境效率,结果如图2所示。

从图2可以看出:当考虑二氧化碳排放对经济增长的约束作用时,我国各省市的经济增长普遍存在着环境无效率,而且环境效率在省际间的分布不均,起伏波动较大。2000年位于技术前沿面上的省份为内蒙古、福建、广东、海南、甘肃和青海,而北京、天津、上海、江苏、浙江等经济发达省市的环境效率仅在0.6~0.7左右,四川、贵州、云南、陕西、新疆等经济不发达的省市环境效率也较低,这表明环境效率似乎与现有经济发展水平没有多大的必然联系。环境效率的省际差异大,但观察可知环境效率高和环境效率低的省份分别在地理位置上相互邻近,表现出高——高、低——低的分布格局,即环境效率具有空间相关性和集聚特征。环境效率普遍低下的原因主要在于:在生产过程中使用的全局生产技术与自身经济发展条件(包括经济增长水平、资本存量、资源环境等)不相匹配,不能充分利用最佳实践技术;经济增长主要是靠能源、资本等要素投入而不是技术效率推动的,在创造大量期望产出的同时也不可避免造成二氧化碳排放等非期望产出,只有通过减少能源消耗投入或二氧化碳排放才能提高环境效率,否则碳排放约束将给经济的可持续发展带来较大的掣肘。

从全局环境效率的平均增长速率来看,大部分省市的技术效率是恶化的,原本位于生产前沿的省市环境效率变得无效,另有一部分技术无效的省市环境效率继续无效。平均而言,全局环境效率改善的省市有11个,内蒙古、广西、甘肃、宁夏、青海等省市的恶化程度较大。最初环境效率有效的甘肃、内蒙古、青海出现不同程度的环境效率恶化,而北京、天津、上海、福建、广东、江苏、山东等省市则由环境效率无效转变为有效或持续保持有效。这表明经济发达省市已意识到资源环境对经济可持续的重要性,减少能源和物质资本的过量投入,引进或创新已有的环境技术,以更为有效地方式组织生产活动,进而实现了环境效率改善。而环境效率恶化的省市多是经济不发达省市,对它们而言经济的高速增长是位于首位的,在生产低水平扩张过程中导致技术——技能匹配度下降。这从侧面印证了“环境库兹涅茨曲线”所描绘的状况,只有积极发展低碳经济发展模式,在促进效率改善上做好文章,才能扭转环境效率恶化趋势和实现经济可持续地健康发展。

3.3 基于GML指数和ML指数的全要素生产率变动分析

为了表现碳减排的治理需要一定的代价,从而将环境规制引入研究框架中,本文利用SBM方向性距离函数和GML指数计算并分解了中国29个省市2001~2009年间的全要素生产率变动、纯效率变动、技术进步、规模效率变动、技术规模变动,分析各省市生产效率的增长特征与差异,并对两种测度结果进行了比较,具体结果如表2所示。

由于基于全局生产技术集的GML指数具有循环可加性,表2中的GML、GPEC、GPTC、GSEC和GSTC数据实际上为各个省市在2001~2009年间的累积变化值,反映了近10年来中国环境全要素生产率、技术效率、技术进步、规模效率和技术规模的总体变动程度。表2显示,在低碳经济视角下,2001~2009年间我国环境全要素生产率总体上呈上升趋势。全国平均环境全要素生产率增长指数为1.096,10年间共增长了9.6%。根据《中国环境经济核算研究报告2008》统计,与2004年相比,2008年我国环境退化成本增长了74.8%,虚拟治理成本增长了75.4%,远远高于同期GDP的增长速度,这意味着伴随着经济的高速增长,中国在节能减排方面所付出的努力日益增大。环境规制是保持我国现有经济发展可持续性的关键,探讨建立发展的同时追求环境保护的长效机制势在必行。

续 表

注:各个指数的变动为期末值/期初值,变动指数的全国均值是各省市变动指标的算术平均。

全要素生产率的增长指数可进一步分解成纯技术效率变化、技术进步、规模效率和技术规模效率变化,从表2可知,10年间节能减排的约束条件使得纯技术效率不仅没有上升,反而恶化了11.9%,而技术水平、规模效率以及技术变动的规模效应都呈现增长趋势,其中技术水平提高了8%,而且技术进步也导致了生产呈现规模报酬递增状态,两者在推动环境全要素生产率增长中占据主导地位。

从环境全要素生产率增长及其分解因子的地区差异来看,技术变动差异最小,规模效率的省际差异程度最大,要素在地区间集聚规模的差异可能是规模效率差异大的原因。规模效率也是造成环境全要素生产率差异的主要原因,技术变动引发的规模效应对环境全要素生产率增长差异的相对贡献也较大,表明规模因素将成为缩小地区间环境全要素生产率差异的关键。东中部地区的环境全要素生产率呈增长趋势,其中东部地区环境全要素生产率10年间增长了22.8%,中部地区增长了4.9%,西部地区的环境全要素生产率呈下降趋势,10年间下降了1.3%。东部地区已经充分认识到经济发展与环境协调的重要性,并不断地通过经济转型降低经济增长的代价,在经济增长质量的提高中达成合作、达到人与自然的和谐,率先体现了可持续增长的内涵。相对而言,中西部地区经济增长效率仍存在相当大的改善空间。其中,中部地区因邻近东部地区,和东部地区的经济合作较为密切和频繁,其发展速度和质量要好于西部地区,表明中部地区是较好的“追赶者”,或已经意识到要注重当地环境的保护,而西部地区经济基础较为薄弱,经济增长使其首要任务,导致西部地区在发展过程中更多地依靠资源的投入,是粗放型增长模式,对当地环境的破坏性较大,阻碍了全要素生产率的增长。东部地区环境全要素生产率的增长主要依靠技术进步推动,中部地区是由技术进步的规模效应起到主导作用,西部地区环境全要素生产率的下降主要是由于纯技术效率的恶化。东部地区的技术进步率最高,中部其次,西部最低。就技术进步的规模效应而言,3个地区的技术边界都偏离不变规模报酬技术,表现出规模经济,推动全要素生产率的增长。其中,西部地区的技术规模效应最高,东部其次,中部最小。这表明西部大开发战略的实施吸引部分要素从东部地区流向中西部地区,西部地区因技术水平低下,要素的流动和初始集聚极大地推动了其技术水平的提升,东中部地区要素集聚规模和结构已经相对成熟,对技术水平的影响相对减少。东、中、西部的纯技术效率在10年间都有所恶化,恶化程度依次为东部地区(4.9%),中部地区(9.4%)和西部地区(21.5%),这与环境全要素生产率增长地区差异的表现形式是一致的,三大地区技术效率的改善有很大空间,将成为提高全要素生产率的关键。进入“十五”规划后,中西部地区规模效率都有所提高,而东部地区的规模效率则下降了4.3%。中西部地区的规模效率虽然都有所增长,但其环境全要素生产率的增长却不是最大,这表明规模效率并不是环境全要素生产率增长的主导,中西部地区较低的技术水平和严重恶化的技术效率是造成其环境全要素生产率增长缓慢的主要原因。东部地区的规模效率下降得最为严重,这是由于东部地区的要素集聚已经出现了集聚过度的特征,进而导致规模经济效率的下降。如东部地区的江苏、浙江、山东、海南等省市都是典型的要素集聚程度极高的地区,但是这些地区的规模效率出现不同程度地下降。中西部地区的规模效率呈现增长趋势,原因主要在于经济落后的省市集聚的要素规模较小,还不能满足现有生产技术对其的要求,因而对已有要素的利用效率较高,呈现规模报酬递增。这说明要素集聚程度过高的地区会造成规模效率的损失,只有不断利用技术水平要素质量和集聚结构,才能实现规模报酬递增,真正意义上提高经济发展质量。

4 结 语

在促进我国经济高速增长的同时也要关注经济增长与能源消耗、环境污染之间的矛盾,这对提高经济运行质量及转变经济增长模式提供了良好的理论基础和政策导向。考虑到Malmquist-Luenberger指数和方向性距离函数不具有循环累加性,不能表现出全要素生产率的长期变动趋势,而且容易导致潜在的线性规划无解问题。且运用径向的、角度的DEA计算方向性距离函数时,经常忽略投入或产出的某一个方面。为了避免出现这些问题,本文在SBM方向性距离函数和GML指数的基础上,测度了碳排放约束下2001~2009年中国29个省份的环境效率、环境全要素生产率增长,其中环境全要素生产率增长分解为技术进步、纯技术效率变动、规模效率变动和技术规模变动,得出的主要结论和启示如下。

能源的过度使用以及二氧化碳的过度排放使我国各省市普遍存在环境无效率,且环境效率在省际间的分布差异较大;低碳经济下我国各省市环境全要素生产率的增长呈现地域不均衡特征,技术水平和规模效率呈现改善趋势,且技术进步也导致了生产呈现规模报酬递增状态,技术进步和技术规模效应在推动环境全要素生产率增长中占据主导地位,而纯技术效率则出现不同程度地恶化;规模效率的省际差异程度最大,是造成环境全要素生产率差异的主要原因,技术变动引发的规模效应对环境全要素生产率增长差异的相对贡献也较大;东中部地区的环境全要素生产率呈增长趋势,西部地区经济基础较为薄弱,经济增长是其首要任务,粗放型的增长模式对当地环境的破坏性较大,导致其环境全要素生产率呈下降趋势;中西部地区较低的技术水平和严重恶化的技术效率是造成其环境全要素生产率增长缓慢的主要原因,而东部地区的要素集聚己经出现了集聚过度的特征,导致其规模效率则有所下降;当期生产技术下的环境全要素生产率普遍低于全局生产技术下的测算结果,这种差异是由全局生产技术下生产前沿持续向外扩展造成的。

上一篇:化学中的兴趣教学下一篇:染色问题论文