网络机群下多项式预处理EBE-PCG并行算法设计与实现

2024-06-27

网络机群下多项式预处理EBE-PCG并行算法设计与实现(精选6篇)

篇1:网络机群下多项式预处理EBE-PCG并行算法设计与实现

网络机群下多项式预处理EBE-PCG并行算法设计与实现

针对单机上实现困难,计算费用高昂的大规模结构动力学问题,本文采用将总体运算分解到单元上进行的`EBE计算策略和基于区域分裂的SBS存储和任务分配策略,设计了粗粒度EBE-PCG并行算法,并在网络机群环境下得以实现.在PCG迭代法中分别采用Jacobi预处理矩阵和多项式预处理矩阵,比较它们的迭代求解效率.悬臂梁受冲击载荷与吉普车车架振动响应分析问题的数值算例,证明了该算法不但能够显著地提高问题的求解规模,适合大规模结构分析计算;而且还能获得良好的并行效率,是一种适合网络机群并行环境的有效的粗粒度并行算法.

作 者:乐志华 程建钢 姚振汉 作者单位:清华大学,工程力学系,北京,100084刊 名:工程力学 ISTIC EI PKU英文刊名:ENGINEERING MECHANICS年,卷(期):19(5)分类号:O242关键词:SBS策略 EBE策略 PCG法 网络机群并行系统

篇2:网络机群下多项式预处理EBE-PCG并行算法设计与实现

连续属性离散化处理是数据挖掘中一个重要且较活跃的研究领域[1]。由于许多算法要求属性数据取值为离散的, 而现实中经常出现连续属性值的情况, 为了满足这类算法的要求, 常常需要对连续属性的取值进行离散化。同时, 高效地离散化处理能够减少算法的时间和空间开销。

针对连续属性的离散化, 典型的方法主要有基于信息熵的方法[2,3,4,5,6,7]和基于聚类的方法[8,9]。基于信息熵的方法较充分地利用了样本信息, 是最常用的离散化方法之一。Clarke和Barton[2]提出的基于熵和最小描述长度理论的离散化方法是最经典的离散化方法之一, 在计算和确定划分点时利用了信息熵概念[3], 这使得它更有可能将区间边界定义在准确位置, 使得离散化的效果更好, 有助于提高分类的准确性。谢宏等人[4]对每一个候选断点定义了信息熵, 并以此作为断点重要性的量度标准, 提出了基于信息熵的粗糙集连续属性离散化算法。Lee[5]使用信息理论衡量每个区间所包含的信息量, 通过最小化信息损失来离散连续属性值, 提出了一种基于Hellinger的连续数据离散化方法。Wu等人[6]基于二分熵的概念和复合分布索引技术, 提出了一种基于分布-索引的离散化器。赵静娴等人[7]分析了基于熵的离散化标准的切点特性, 提出并证明了一种基于边界点属性值合并和不一致度检验的离散化算法。但是, 这种离散化的方法仅仅局限于小规模的数据集, 对于海量数据的连续属性离散化处理, 将会是一个很大的瓶颈。

面对海量数据, 并行化是提高其效率的有效手段, 而Hadoop是并行计算的一种有效实现方式。因此, 本文在MapReduce编程模型下, 提出一种基于熵的连续属性离散化并行算法, 该算法既保证了离散处理后的分类效果, 又解决了对于海量数据离散化处理效率低的问题。通过大量的实验结果表明, 提出的并行化算法在海量数据集的离散化处理上有很高的效率。

1 相关工作

1.1 MapReduce编程模型

Hadoop是一个开源的框架, 主要由分布式文件系统 (HDFS) 和MapReduce编程模型两大核心组成[10]。HDFS能够实现数据在计算机集群组成的云上高效的存储和管理。MapReduce是一种并行编程模式, 用于海量数据集的并行计算, 它将该并行计算过程高度地抽象为了两个函数:Map函数和Reduce函数。因此, 在MapReduce框架中编写应用程序就是定制化mapper和reducer的过程。这种计算模型要求输入、输出和中间数据的形式都是一组键值对。Map函数接收一个形式的输入, 然后同样产生一个形式的中间输出, Hadoop会负责将所有具有相同中间k2值的v2集合到一起传递给Reduce函数, Reduce函数接收一个如形式的输入, 然后对这个value集合进行处理, 每个Reduce的输出也是的形式。最后, 将结果返回给应用程序。

1.2 基于熵的连续属性离散化算法

在具有n个连续属性的数据集M中, 对于任意一连续属性, 其最大区间数为其样本个数k, 即一个区间至少有一个样本值。将连续属性离散化, 就是将该连续属性的取值范围划分为若干个区间, 将连续的样本值转换成少于k的若干个离散值。因此, 离散化算法的关键在于合理确定划分点的个数和位置[2]。

对于离散型随机变量A的熵定义如下[2]:

文献[2]中已证明H (A) 是一个凹函数。可通过合并区间, 使得区间数k减少从而有效地减小熵值H (A) , 但最终需要达到一个最小熵值和最优区间数k的平衡。假设区间数k的最小取值能够达到1, 即将所有值都合并为一个区间, 依据文献[2]中的平衡条件可得如下公式:

其中:kmax是最大区间数, H (kmax) 是最大区间数时的熵值。当Ck>Ck-1时, 停止区间合并。此时, k为最优区间数。

在此基础上, 本文将基于熵的离散化方法应用到Hadoop平台上, 通过对连续属性值的自身合并, 找到合并前后熵差最小的合并点作为最佳合并点;然后使用式 (2) 作为算法的停止条件。文献[2]中的算法可简单描述为:

输入:有多个连续属性的数据集M

输出:数据集M的各个属性上的合并状态集

步骤:

对于M的每一个属性A:

(1) 计算初始状态熵值H (kmax) ;

(2) 先将M按照属性A进行升序排列, 然后在有序序列中选取相邻两个毗邻区间合并;

(3) 逐一计算每相邻两个区间合并后的熵值H (k-1) , 将其与H (kmax) 作差, 选择熵差最小的合并区间作为最佳合并区间;

(4) 通过式 (2) 计算Ck与Ck-1的值并比较大小, 若Ck

2 基于Hadoop的并行离散化算法

基于MapReduce编程模型的并行算法设计, 最主要的工作是设计和实现Map和Reduce函数。由算法描述可知, 数据集M中的n个属性之间可以并行计算最佳合并区间, 因此, 该算法的MapReduce并行化具有可行性。但在设计离散算法并行时需要计算熵值, 为了提高该离散算法在并行时的效率, 将信息熵公式中求概率p (a) 时所需要的信息进行预处理。因此我们将离散化过程分为两个阶段:离散化预处理阶段和离散化处理阶段。

2.1 离散化预处理阶段

本阶段所要做的工作主要是将数据集M上的各个属性Ai (i=1, 2, …, n) 的所有可能的取值, 及其该值出现次数num统计出来, 为求出连续属性Ai的所有取值的概率做准备。离散化预处理阶段的伪代码可描述为:

2.2 离散化处理阶段

经过离散化预处理后, 将其输出结果作为离散化处理的输入。在该阶段, Map函数将仅仅只是改变其输出形式, 将key中attribute Id的内容即列号取出来作为Map的key输出, 这样一来, 才能在Map与Reduce之间的洗牌阶段将不同的属性列分别传到不同的Reduce节点上, 实现各连续属性之间并行离散化的目的。离散化处理阶段的伪代码可描述为:

将Map输出的不同的key值分别映射到不同的Reduce任务节点上, 实现对数据集中各个连续属性的并行离散化处理。该机制可通过修改mapred_site.xml配置文件中mapred.reduce.tasks的属性值为4来实现。

根据Reduce函数的输出, 可以得到所有区间值, 从而完成离散化处理过程。

3 实验及结果分析

3.1 实验环境

本文进行实验所使用的Hadoop集群中共有37个节点, 每个节点的配置信息及运行环境完全相同。操作系统版本是ubuntu10.12, Hadoop版本是0.20.2, Java版本是1.6.0_21。根据Handoop项目官方网站介绍的方法配置基于hadoop0.20.2版本的集群[11]。其中一台作为Master节点, 其余36台全为Slave节点。

在集群中测试离散效果时使用的基础实验数据来自于UCI dataset中的wine数据集。在测试算法的性能时, 实验中分别构造了256MB、512MB、1GB、2GB和4GB等5个不同大小的数据集。

3.2 实验结果

在保证了离散化效果的前提下, 本文做了大量的实验工作。分别用256MB、512MB、2GB、4GB等不同规模的数据集在2、4、6、8、10、12, 14、16个不同节点的集群上做测试, 得到的实验结果如图1所示。

由图1可以看出, 在分布式集群环境中, 对于小规模的数据集效果甚微, 但对于大型数据集文件, 有明显的优势。该实验结果验证了该算法更适用于规模较大的数据集, 使海量数据的离散化处理成为可能。

同时, 以2GB数据集文件为例, 分别在2、4、6、8、10、12个不同节点的集群上做实验, 并计算了加速比, 加速比可定义为:在单机上的运行时间/集群中的运行时间, 实验数据如表1所示。

将表1中的数据绘制成折线图, 如图2所示。可以明显地看出, 随着节点个数的增加, 加速比呈线性增长, 即对于具有一定规模的同一数据, 随着节点个数的增加, 算法所运行的时间明显减少, 算法的运算速度明显提高。

但在图2中可以看出, 8个节点后, 加速比的增长速度明显减小, 就意味着时间减少的趋势逐渐变缓。这是因为, 节点之间也存在信息的交互, 集群中的节点越多, 节点间的通信量就越大, 使得时间减少的缓慢, 可以通过增加网络带宽来适当解决此问题。对于一定规模的数据集, 随着节点的增加, 运行时间会减少, 但并不是集群中的节点个数越多, 运行时间就会一直减少。不同文件所需要的最佳集群节点个数如图3所示, 可以观察到同一数据文件, 随着节点个数的增加, 运行时间减少的趋势变缓, 并会在某个节点集群之后不再有变动。这是因为Hadoop对数据是分块处理的, 默认数据块的大小是64MB。对于512MB文件, 可划分为8个数据块, 又由于在集群环境配置中设置的每个节点最大运行的map任务数为2, 因此512MB文件在集群节点数为4的时候效果达到最佳。再比如, 1GB的文件在16个节点运行的效果最佳, 在增加节点的个数后, 将趋于平缓。

4 结语

本文对Hadoop环境下连续属性离散化的并行实现做了深入研究, 在MapReduce编程模型下提出了基于熵的连续属性离散化的并行方法, 并予以实现。首先, 简单介绍了Hadoop的体系结构和基于熵的连续属性离散化算法, 然后设计并实现了基于MapReduce编程模型的并行算法。最后, 通过实验结果和实验分析得出, 在保证离散效果的前提下, 该方法可以高效地处理海量数据的连续属性离散化。

本文提到的基于熵的离散化算法虽然能产生较好的离散效果, 但它是一种无监督的, 自下而上的离散化算法。对于有类标记的数据集来说, 将其连续属性离散化应考虑类标记对其连续属性的影响, 因此, 下一步将研究更高效的离散化算法, 提出一种有监督的高效的连续属性离散化算法, 并应用到云计算环境中。

参考文献

[1]Kurgan L A, Cios K J.CAIM discretization algorithm[J].IEEE Transactions on Knowledge and Data Engineering, 2004, 16 (2) :145-153.

[2]Clarke E J, Barton B A.Entropy and MDL discretization of continuous variables for Bayesian belief networks[J].International Journal of Intelligence Systems, 2000, 15 (1) :61-92.

[3]Dougherty J, Kohavi R, Sahami M.Supervised and unsupervised discretization of continuous features[C]//Proceedings of the Twelfth International Conference Machine Learning, 1995:194-202.

[4]谢宏, 程浩忠, 牛东晓.基于信息熵的粗糙集属性离散化方法及应用[J].计算机学报, 2005, 28 (9) :1570-1574.

[5]Lee C.A hellinger-based discretization method for numeric attributes in classification learning[J].Knowledge-Based Systems, 2007, 20 (4) :419-425.

[6]Wu Q, Bell D, Prasad G, et al.A distribution-index-based discretizer for decision-making with symbolic AI approaches[J].IEEE Transactions on Knowledge and Data Engineering, 2007, 19 (1) :17-28.

[7]赵静娴, 倪春鹏, 詹原瑞.一种高效的连续属性离散化算法[J].系统工程与电子技术, 2009, 31 (1) :195-199.

[8]Gupta A, Mehrotra K, Mohan C.A clustering-based discretization for supervised learning[J].Statistics&Probability Letters, 2010, 80 (9-10) :816-824.

[9]Ferreira A J, Figueiredo M A.An unsupervised approach to feature discretization and selection[J].Pattern Recognition, 2012, 45 (9) :3048-3060.

[10]Chuck Lam.Hadoop实战[M].北京:人民邮电出版社, 2011.

篇3:网络机群下多项式预处理EBE-PCG并行算法设计与实现

关键词 BOC 捕获 频谱 仿真搜索 正交

1 引言

由于GPS卫星发射的导航信号比较微弱,而且以固定的频率发射,因此很容易受到敌方的干扰。为了改善这种局面,1977年美国正式提出了“GPS导航战”(Navigation warfare)的概念[1-3]。但由于现行的军用和民用GPS信号的频谱是混叠使用的,所以美方对C/A码施放干扰信号也必然导致军用信号频谱的中央部分不能使用。为了改善GPS这种频谱结构的缺陷,美方经过大量的研究后,确定了最终方案:在现有L1和L2两个频段上附加新的军用信号,即M码。M码采用一种裂谱调制法,它把其大部分功率放在靠近分配给它的频带的边缘外,抗干扰能力主要来自不干涉C/A码或Y码接收机的强大的发射功率,而且M码信号的保密设计基于下一代密码技术和新的密钥结构,具有很高的保密性。实现M码分裂频谱调制方法有很多种,经过大量比较,确定了选用二进制偏移载波(BOC)调制技术,这是因为它不仅使M码信号的频谱能与现有信号的频谱很好的隔离,而且还有更好的性能。本文在介绍BOC调制技术的原理和分析其频谱结构的基础上,提出了基于并行码相位搜索的改进捕获算法,设计了正交IQ捕获算法,对两路并行捕获,捕获性能优于传统的捕获方法,仿真结果显示算法的正确性。

2 BOC调制信号的频谱与自相关特性

二进制偏置载波(BOC)调制,该调制提供两个相互独立的设计参数:

(1)副载波频率fs,单位MHz;

(2)扩频码速率fc,单位Mchip/s。

通过这两个参数,可以自由地将信号能量集中在所分配带宽的指定部分,以减少接收其他信号的干扰。

BOC调制信号通常被记作为BOC(fs,fc)或BOC(a,b),其中fs=af0,fc=bf0(a,b通常为正整数),基准频率f0=1.023MHz。

2.1 BOC调制信号的频谱特性

BOC(fs,fc)信号的归一化基带功率谱4-7:当n=2fs/fc为偶数时[4-7],

(1)

当n=2fs/fc为奇数时,

(2)

采用BPSK调制方式,扩频码码率为fc的GPS C/A码信号的归一化基带功率谱为:

(3)

图1所示为 BOC(1,1)信号和GPS C/A码信号的归一化基带功率谱,可以看出BOC(1,1)相对于 GPS C/A码信号,最大的一个特点就是具有中心分裂的功率谱,这将有助于降低它与GPS C/A码信号的干扰,并且由于它在高频段具有更多的能量,它比GPS C/A码信号拥有更大的Gabor带宽, 这也意味着它具有更佳的捕获性能。

2.2 BOC调制信号的自相关特性

BOC调制信号的自相关函数,可以通过对BOC调制信号的功率谱做逆傅立叶变换得到。理想情况下(接收机前端带宽无限宽)BOC(pn,n)信号的自相关函数为[8]:

(4)

其中,p=1,2,3,…,k=■。

在(4)中令p=1,可得BOC(1,1)自相关函数RBOC (1,1)(τ)为:

(5)

采用BPSK调制方式的GPS C/A码信号的自相关函数RCA(τ)为:

其中k=2τ/Tc,Tc=1/fc。图2所示为BOC(1,1)信号和GPS C/A码信号的自相关函数,可以看到BOC(1,1)的自相关函数相对于GPS C/A码信号而言,具有多个峰值,这也对接收机提出了更高的要求。通常将自相关函数正中间的峰值称为主峰,而其它的峰称为旁峰。

3 捕获算法的设计

经典的并行码相位搜索框图如图3所示[9-10]:

本地产生的副载波与本地产生的伪码相乘之后变换到频域取共轭,输入信号经过傅立叶变换后与经过傅立叶变换的BOC码相乘,输出结果经过傅立叶逆变换 转换为时域信号,傅立叶逆变换输出的模值表示输入信号与本地伪码的相关结果。改进后的并行码相位算法框图如图4所示:

改进捕获算法如下:

(1)将接收到的中频信号取1ms的数据长度,记为x(n);

(2)本地载波取中频信号IF±10KHz范围,步长取500Hz,生成41个本地中频载波,并对生成的本地载波进行采样,记为li(n)(i=1,2,…,21);

(3)将x(n)和li(n)的转置点对点相乘,结果取FFT,记为Fi(k);

(4)产生本地副载波,相位分别为0°和90°,两路信号分别进行FFT变换到频域中,并取共轭。记为f1*(k)和f2*(k);

(5)Fi(k)分别与f1*(k)和f2*(k)点对点相乘,结果记为G1(k)和G2(k);

(6)产生本地的C/A码,并对其进行采样,记为Cs(n)(s为卫星编号);

(7)对Cs(n)进行FFT,变换到频域中,并取共轭,记为Cs(k)*;

(8)将G1(k)和G2(k)分别与Cs(k)*点对点相乘,结果为Rs1(k)和Rs2(k);

(9)对Rs1(k)和Rs2(k)进行IFFT,变换到时域中,得到rs1(n)和rs2(n),得到其模平方rs1(n)2和rs2(n)2;

(10)比较rs1(n)2和rs2(n)2的最大值,选择判决那个大就选择哪个,记为rsi(n)2,最大值对应的坐标即为C/A码的起始点和载波粗频。

4 仿真结果及分析

高斯信道信噪比为0dB时的捕获结果如图5所示,局部放大图如图6所示。由图可以看到有1个主峰和2个次峰,这和BOC(1,1)信号的自相关函数有三个峰值是相符的,由图BOC(14,2)捕获图可以看出,峰值位于码相位45094和频率25.542MHz处,图7为BOC(14,2)信道加高斯干扰0dB-34dB的仿真结果。

由图仿真结果可知,当干扰为-36dB时,捕获算法完全失效,主峰与副峰比值基本相同,捕获结果完全被信号淹没,在干扰为-35dB到-27dB时,捕获算法仍然有效,但效果较差,此算法在大于-27dB高斯信道干扰下完全可以捕获到信号的初始码相位和载波频率,为进一步跟踪提供了粗略值。

5 结束语

改进的捕获算法基于并行码相位搜索算法,针对收发端载波不同的问题,在产生本地副载波时分别产生I路和Q路两路信号,分别基于捕获算法,最后对捕获到的峰值进行比较,这样可以选择最优的捕获结果。

参 考 文 献

[1] 杨东凯,樊江滨,张波,张敏译. GNSS应用与方法. 电子工业出版社. 2011.09 p.16-22

[2] 杨力,基于BOC调制的导航信号同步接收关键技术研究. 南京理工大学.2009.12

[3] 楚横林,李春霞. BOC调制导航信号关键技术研究. 2010 Radio Engineering Vol.40 No.6

[4] 王智,耿相铭. Galileo卫星导航系统中的BOC调制与接收技术.上海航天.2007年第6期

[5] Peter Rinder, Nicolaj Betelsen, Design of a single frequency GPS software receiver. Aalborg University 2004.

[6] D Akopian.Fast FFT based GPS satellite acquisition methods.IEE Proceedings Radar, Sonar and Navigation,2005,1 52(4):277—286.

[7] 杨俊,武奇生. GPS基本原理及其Matlab仿真. 西安:西安电子科技大学出版社. 2006:126-181页

[8] Julien O,Cannon M E,Lachapelle G,et al.A New Unambiguous BOC(n,n) Signal Tracking Technique [A].In:Europe Navigation Conference GNSS 2004[C].2004:17-19P

[9] 刑兆栋. BOC 导航信号处理方法及卫星导航数据的传输. 北京航空航天大学. 2006:58-73页

[10] Kai Borre,Dennis M.AKos,Nicolaj Bertelsen Peter Rinder Soren Holdt Jensen,杨东凯,张飞舟,张波译. 软件定义的GPS和伽利略接收机. 国防工业出版社.2009.3

Based on Parallel Code Phase Search of BOC Signal Acquisition Algorithm of Design and Implementation

Xin Fuguo1,Li Rongfang2

(1. Shaanxi Post and Telecommunication College department of communication,Xianyang 712000,China)

(2. Shaanxi Post and Telecommunication College department of computer,Xianyang 712000,China)

Abstract BOC modulation,as a navigation signal system ,has such better performances as multipath-resistant and anti-jamming ,it can realize spectrum separation under the condition of carrier frequency being shared .This paper first analyzes the characteristic of the frequency spectrum of BOC signal and the autocorrelation function , put forward based on the concurrent code phase search of the algorithm, the design of orthogonal orthogonal IQ acquisition algorithm, the two road parallel acquisition, the simulation results show the correctness of the algorithm and the superiority.

篇4:网络机群下多项式预处理EBE-PCG并行算法设计与实现

现代空间遥感技术的发展已经使传感器的空间分辨率、光谱分辨率和时间分辨率大幅度提高, 使得遥感图像的数据量大增, 有的影像文件大小甚至达到GB级[1], 与此同时, 人们对遥感影像的处理和从影像获取信息的速度却远远落后于影像数据的增长速度, 其原因主要有以下两点: (1) 遥感影像数据量过大, 所需处理图像点的个数过多; (2) 另一个很重要的原因则是遥感影像的处理和信息提取算法越来越复杂, 一次算法的执行往往需要对同一幅图像做多次遍历和处理。对于第一个问题, 可以通过对遥感影像进行多尺度影像分割的方法来降低单幅影像的处理时间[2]针对后一个问题, 在硬件条件和遥感影像大小一定的情况下, 如何提高遥感产品算法的执行速度和并行度, 充分利用计算节点CPU的运算能力, 成为解决问题的关键;在复杂的遥感产品算法执行中利用并行计算技术, 需要从算法本身入手, 对这些算法进行细化和拆分, 以使得对拆分后得到的原子算法集合进行并行执行, 遗憾的是, 现有文献对遥感图像并行处理研究多集中在影像处理的基础算法方面, 比如遥感图像辐射校正、几何校正、图像融合等, 但是这些算法恰恰是遥感图像处理所用到的最基本的算法, 很难拆分出可进行并行的子算法集合;对于复杂的遥感产品算法并行计算研究的比较少, 本文从复杂的遥感产品算法的分解以及遥感产品算法树的构建和执行处着手, 着重分析了基于算法生成树的复杂遥感产品算法的并行计算方法, 以实验数据验证了此方法的可行性, 最后描述了此策略在实际项目中的应用。

2、算法生成树的构建

一个复杂的遥感产品算法的执行, 往往需要很多依赖的遥感产品算法的支持, 我们称这些遥感产品算法为中间算法, 个最终产品的生产往往需要几个低层次的图像处理得到的结果来作为依赖产品, 一般来说越是高级的图像处理功能, 需要越多的的底层图像处理功能的叠加来实现。一幅原始影像数据的处理, 一般要经过以下三个步骤: (1) 图像预处理, 比如几何纠正, 大气纠正, 辐射纠正等均是常用到的图像预处理算法, 这些算法也被称作基础算法。 (2) 图像中间处理, 这里分为两种:第一种我们称之为基础中间算法, 许多中间算法均需要这些算法的计算结果做为参数, 常用到中间基础算法有CCD表观辐亮度 (CCDTOAL) 、CCD表观反射率 (CCDTOAR) 、CCD太阳天顶角 (CCDSOLARZA) IRS像元纬度 (IRSPixelat) 等等;另一种中间算法处于基础中间算法和最终算法之间, 它为最终算法的执行提供所必须的参数。例如归一化植被指数 (NDVI) 、植被覆盖度指数 (VCI) 、比辐射率 (EP) 等 (3) 最终处理在得到最终产品算法执行所需的一系列参数后, 调用最终的遥感产品算法对图像进行最后的处理。一个复杂遥感产品算法, 其实可以看做是一棵由以上四种算法组合得到的算法生成树, 图1展示了一个复杂遥感产品算法执行的各个阶段以及其产出。

下面我们以地表温度产品 (LST) 为例, 来描述算法生成树的构建过程:

地表温度产品的计算公式如下所示:

其中λ、K1为常数, EP和IRSTOAL分别是比辐射率产品算法 (EP) 和基础中间产品算法表观辐亮度 (IRSTOAL) 的执行结果, 因此可知道, LST直接依赖于EP和IRSTOAL, 对于EP来说

其中NDVI、NDVImin、NDVImax均为归一化植被指数产品算法的执行结果 (NDVI) , 则EP依赖于NDVI, 同时EP的计算结果需要根据CCD表观反射率 (CCDTOAR) 的计算结果做校正, 因此EP依赖两个产品算法NDVI和CCDTOAR, CCD表观反射率的计算需要CCD表观辐亮度 (CCDTOAL) 和CCD太阳天顶角 (CCDSolarZA) 的计算结果, NDVI需要CCD表观辐亮度 (CCDTOAL) 的计算结果, 因此可得出LST的算法生成树如图2所示:

3、算法执行策略和并行度分析

在上节中, 以地表温度产品算法为例讲述了算法生成树的构建过程, 在本节中将重点讲述如何以并行方式执行算法生成树

算法生成树的执行遵循如下三个原则

(1) 启动算法调用时, 首先从从最左侧叶子节点开始调用, 同一级别的算法可并行执行然后逐层向上调用, 每一层都是从左到右执行, 启动一个线程调用执行

(2) 只有当一个算法所有依赖算法全部执行完毕, 才执行此算法每个算法调用完毕后生成的与算法同名的结果 (CSV或tif) , 均生成与其算法名对应的图像文件或CSV文件, 供上层算法调用。

(3) 在调用过程中, 如果同名的算法正在被执行或已执行完毕, 则跳过此算法的生产。例如, 图2第四层有两个CCD表观辐亮度产品算法 (CCD-TOAL) , 在执行时, 只需要执行最左边的即可, 这样的执行方式大大节约了计算资源, 缩短了产品算法的执行时间。

仍以地表温度产品为例, 根据上面所述的三个原则, 地表温度产品算法生成树的调用顺序如下图所示:

由于同一级的产品算法可并行执行, 并且同名产品算法的重复执行, 在这里简单的假设每个产品算法的的执行时间均为T, 那么对于算法生成树的执行总时间为TA=T+T+T+T=4T, 而常规的产品算法, 由于未用到并行技术, 执行所用的时间为每一个算法执行的时间总和TC=8×T=8T用时比S=TA/TC=2.63, 当然, 这里, 只是在理论上进行了一个估计, 在算法的实际执行中, 由于不同产品算法执行所用的时间不同, 包括并行执行对计算节点CPU造成的影像, 对比的结果会有一定的差别, 在下节中, 将通过仿真实验进行详细的比较%

4、仿真与分析

为能真实反映本文提出的算法生成树的产品算法执行方法在执行时间和效率方面对比原有产品算法的优势, 本文仿真实验采用进行如下仿真实验:

实验采用8台相同配置的PC机, 其硬件配置如表1所示

仿真实验采用的原始影像为2011年11月1号环境星拍摄的北京区域大小为5898×5987的遥感影像, 每个计算节点对此幅影像计算地表温度 (LST) ,

测试采用两种方式来进行, 第一种采用本文所述的算法生成树并行的方式调用执行此过程, 另一种调用常规的地表温度算法, 8台PC机分作两组, 每四台一组, 采用向计算节点分发任务的方式, 每个任务为一次LST的执行过程, 一次性的向每台PC机分发五个任务, 然后计算出两种执行过程所用的平均时间和CPU平均利用率。

从实验结果可以看到, 由于采用了并行计算技术, 采用基于算法生成树的遥感产品算法的平均CPU利用率约是常规算法的3倍, 由于单位时间内处理的图像点的个数比常规算法的要超出很多, 因此, 算法执行平均所用时间也比常规算法节约了三分之二, 由此可见, 本文提出的基于算法生成树的遥感产品算法调用方法在用时和CPU利用率方面的表现明显比常规算法更为优异。

5、项目应用实例

本文提出的算法生成树的思想已经在遥感环境产品监测与生产分系统中得到了应用, 该系统实现了对200多种遥感产品的并行生产, 其连续处理能力能够满足海量遥感数据处理的需求, 下面, 给出地表温度产品 (LST) 并行生产的监控示意图 (如图4) :

6、结束语

本文针对当前遥感影像处理中存在的算法执行效率不高, 执行时间过长的问题进行了深入的论述和研究, 提出了算法生成树的概念, 并论述了其并行执行策略, 实现了对复杂遥感产品算法的可并行化生产, 并使得算法执行的时间和效率得到一定的提高, 该算法已经在实际项目中得到了应用, 对遥感影像处理的并行化和集群化提供了强力的依据, 尽管如此, 本文提出的遥感产品算法生成树的并行执行策略仍有改进的空间, 在下一步的研究中, 将进一步进行算法优化, 同时还将对遥感数据的并行生产进行研究, 以期实现多级并行, 在计算性能上得到更大的提高。

参考文献

[1]郑友华, 王晓青, 窦爱霞.并行计算在遥感震害中的应用[J].地震, 2011-01-10, 31 (1)

[2]沈占锋, 骆剑承, 陈秋晓.高分辨率遥感影像并行处理数据分配策略研究[J].哈尔滨工业大学学报, 2006-11-30, 38 (11) :1968-1971.

[3]曾志, 刘仁义.一种基于分块的遥感影像并行处理机制[J].浙江大学学报, 2012-03-15, 39 (2) :135-141.

[4]樊兴, 张国平, 夏学知.一种实现动态负载均衡的多机情报并行处理方法[J].解放军理工大学学报, 2009-12-15, 10 (6) :523-527.

[5]桂兵祥, 何健.基于高性能云的分布式数据并行处理机制[J].武汉工业学院学报, 2010-03-15, 29 (1) :60-64.

[6]张树凡, 余涛, 李家国.基于三级并行的遥感业务化处理系统研究[J].计算机工程与设计.2012-10-18, 33 (10) :3828-3832.

篇5:网络机群下多项式预处理EBE-PCG并行算法设计与实现

高光谱遥感影像分类算法根据算法思想可以分为两种,一种是基于地物物性的分类方法,主要是利用反映地物物理光学性质的光谱曲线来识别;另一种思想是基于图像数据的分类方法,主要是利用数据的统计特性来建立分类模型[1]。高光谱遥感图像分类技术在高光谱技术研究领域占据着十分重要的地位,如何快速有效的完成高光谱图像的分类一直是人们研究的重点。目前,国内外对提高高光谱遥感影像分类速度的研究主要集中在三方面,分别是针对分类流程提高预处理或准备工作速度,针对算法和数据结构减少运算时间,通过算法并行化提高处理速度。前两种加速方法虽然可以在一定程度上加快分类总体流程速度,但仍不能满足实际应用的需求。采用并行处理是一种在不降低运算精度的情况下更为有效的手段。通过对分类算法的并行化,可以成倍甚至几十倍的提高分类处理的速度,实现大规模、高吞吐率的对高光谱遥感影像数据分类。

本文基于CUDA和GPU环境,设计并实现了高光谱遥感蚀变填图的SCM并行算法。实验结果表明,并行加速比可达到25,验证了遥感图像分类处理并行化的可行性和有效性。

1 基于GPU的SCM并行算法设计

1.1 GPU原理

计算机图形处理器(Graphics process Unit,GPU)是一种专门用来处理个人计算机,工作站或者游戏机上图形运算的微处理器,其设计的初衷是用于在显示器上渲染计算机图形加速应用程序中的图形绘制运算。由于GPU具有强大的浮点计算能力且其像素渲染程序能够对纹理的多个像素进行并行计算,因此GPU在用于通用并行计算方面有着巨大的潜力。但是在CUDA架构出现之前,开发人员只能通过GPENGL或者DirectX等API来访问GPU,这不仅需要开发人员具备一定的图形编程基础而且需要开发人员将通用计算问题转换为图形计算问题,这大大的限制了GPU在通用并行计算领域的推广[2,3]。CUDA架构的出现很好的解决了上面的问题,在该架构下开发人员可以在不具备图形学知识的前提下通过CUDA C(对标准C的简单扩展)对GPU进行编程[4,5]。本文正是利用CUDA构架完成了并行程序的编写。

1.2 SCM算法简介

SCM(光谱相关系数)是一种常见的高光谱图像监督分类算法,它通过计算像元光谱矢量与波谱库光谱矢量之间的相关系数来确定二者之间的相似度,相关系数越大,说明两者越相似。假设n个波段的两种光谱矢量分别为T=(t1,t2,…tn),R=(r1,r2,…,rn),则相关系数r的求解公式为:

1.3 SCM并行算法设计

本文的主要思想是利用CPU与GPU协同编程模型[6]实现SCM算法的并行化,CPU处理串行任务,包括高光谱影像数据的读取,GPU显存的分配,填图等工作,GPU处理并行任务:波谱库数据并行产生参考光谱,光谱相关系数的计算。

SCM算法具有非常高的数据并行性,原始图像上每个像元与训练样本集的计算过程是完全独立的,即像元间的计算过程能够并行化;本文对SCM算法的并行化分两步进行:1:由波谱库数据并行产生参考光谱2:以像元点为基准并行计算相关系数。下面分别描述。

2.3.1 基于波谱库的参考光谱并行算法

由于采用的波谱库中波段数据范围和实际的高光谱遥感影像波段数据范围并不匹配,因此需要利用波谱库数据进行插值得到和实际高光谱遥感影像波段数据相匹配的数据。本文采用线性插值,设原高光谱遥感影像在波长为λ时进行了采样得到的反射率值为a记为(λ,a),在波谱库中有两点(λ1,a1)(λ2,a2)满足:,λ!(λ1,λ2)则插值出的采样点为(λ,(a2-a1)(λ-λ1)/(λ2-λ1)+a1)。假设原始高光谱影像有N个波段,波谱库中的矿物有M种,则将进行N*M次插值运算,注意到针对每种矿物而言其插值运算是独立的,因此此处采用的并行策略是在GPU中产生M个线程,每个线程负责一种矿物参考光谱的插值运算。

2.3.2 相关系数并行计算

通过公式(1)可以知道各个像元点与波谱库矿物的相关系数计算是相互独立的,本文采用的并行策略是以像元点的数目为基准在GPU中产生相应的线程数,每个线程负责一个像元点相关系数的计算,并根据计算的结果对该像元点分类。

假设遥感影像像元点数为N,则GPU中有N个线程,每个线程计算一个像元点相关系数。

3 实验与结果

3.1 实验环境及实验数据

本文进行实验的计算机的配置是:intel core(TM)i7cpu,2.80GHz处理器,4.0GB内存,NVIDIA GeForce 310M显卡,Win7 32位操作系统。

本文采用的CPU编程环境为Microsoft visual C++2008,GPU编程语言为CUDA,绘图采用的是OpenGL,采用的实验数据是江西某地的高光谱遥感数据,数据格式为envi标准格式。

3.2 实验结果

本文分别实现了SCM的串行和并行算法并统计了各自的运行时间,详细结果见表1,图3。

4 结论与讨论

本文采用的加速比评价指标是指在同一问题规模的前提下,GUP程序和CPU的运行时间比。通过表一的实验数据说明,利用GPU进行高光谱遥感图像相关系数填图具有显著的加速效果,当数据量加大时通常情况下加速效果更为明显。当数据量为158*382*1092时加速效果达到了最好的25.68倍。虽然本文采用并行策略取得了较好的效果但是在有效利用共享内存和数据对齐方面任有可以改进和探讨的地方,有待进一步的完善。

摘要:由于空间和波谱分辨率的不断提高,高光谱遥感影像的海量数据特性导致高光谱遥感影像并行处理成为遥感影像处理技术的发展趋势。本文基于CUDA和GPU环境,设计并实现了高光谱遥感蚀变填图的SCM并行算法。实验结果表明,并行加速比可达到25,SCM并行算法能有效改善算法性能。

关键词:GPU,CUDA,并行计算,SCM

参考文献

[1]童庆喜,张兵,郑兰芬.高光谱遥感原理技术与应用[M].北京:高等教育出版社,2006.6.

[2](美)桑德斯著,聂学军译.GPU高性能编程CUDA实战[M].北京:机械工业出版社,2011.1.

[3]吴恩华.图形处理器用于通用计算的技术、现状及其挑战[J].软件学报,2004,15(10):1-3.

[4]龚敏敏.GPU精粹2—高性能图形芯片和通用计算编程技巧(译)[M].北京:清华大学出版社,2007.

[5]NVIDIA CUDA Programming Guide[M].Version1.1.NVIDIA Corporation.2007.

篇6:网络机群下多项式预处理EBE-PCG并行算法设计与实现

本文针对语种识别系统算法的特点,设计了一种基于TI多核处理器TMS320C6678的语种识别并行实现方法,实现了任务级的并行流水处理和核间的高效通信。

1 平台介绍

TMS320C6678是基于TI公司最新DSP系列器件TMS320C66x、采用8个1.25 GHz DSP内核构建而成的业界首款10 GHz DSP,可在10 W功耗下实现160 GFLOP(Giga-Floating Point Operations per Second)浮点计算性能[3]。不仅能整合多个DSP以缩小板级空间并降低成本,同时还能减少整体的功耗要求,充分满足现代数字信号处理日益增长的需求。

本文语种识别系统的开发在TI公司的最新DSP集成开发环境CCSv5(Code Composer Studio)中基于浮点运算设计完成。

2 基于TMS320C6678的语种识别算法优化

2.1 语种识别算法分解

本文的语种识别系统是基于区分性Model Pushing算法[4]进行构建的,并且对特征参数进行了f DWNAP[5,6]处理,因此系统的测试阶段由特征提取模块、f DWNAP模块及对数似然得分模块3个模块构成,如图1所示。

(1)特征提取模块

特征提取模块的任务包括语音信号预处理、MFCC提取、RASTA滤波、SDC扩展、VAD检测、CMS处理、高斯化等过程,该模块结束即输出56维的特征参数,其需要存储的参数包括汉明窗和梅尔滤波器组总共不到2 KB。

(2)f DWNAP模块

该模块的工作是对所提取的56维特征参数进行处理,以去除与语种无关的各种干扰信息,达到净化语种特征参数的目的。如参考文献[6]介绍,该模块首先将特征参数映射至SVM的高维空间,然后利用训练得到的投影矩阵计算映射后的参数中所包含的干扰信息,再将干扰信息映射至特征空间,从而进行去除。该模块中事先训练得到的投影矩阵P=I-ww T,ww T是对称矩阵,因此存储ww T需要7 MB的存储空间。另外,K-L变换矩阵D是对角矩阵,需要112 KB的存储空间。

(3)对数似然得分模块

如参考文献[4]所述,本模块主要任务是利用训练得到的各语种GMM模型对语音特征参数计算对数似然得分进行输出的判决。

本模块需要存储训练阶段得到的各目标语种的GMM模型及非目标语种的GMM模型,即针对每个语种需要存储2个GMM模型。所有的GMM模型只是均值矢量不同,高斯混元权重及协方差矩阵都是共享UBM模型的。以L个语种为例,需要存储2L个均值矢量,即需要224L KB的存储空间,共享的高斯混元权重需要2 KB的存储空间,协方差矩阵由于是对角化的只需要112 KB的存储空间。

2.2 算法实时性分析

首先对各模块的运算实时性进行分析。以30 s的语音(8 000 Hz采样,帧长25 ms,帧移10 ms)为例,后端模型使用单个语种模型,利用CCSv5的环境进行软件仿真得到各模块处理所花的时钟周期数,然后按照TMS320-C6678芯片的单个内核的工作主频(1.25 GHz)计算得到处理时间,结果如表1所示。

由表1可知,整个语种识别系统测试阶段,在算法代码未经任何优化的情况下,一段30 s的语音在单个TMS320C66x CPU内核上的处理时间约为22.3 s,结果非常不理想,并且特征提取模块和对数似然得分模块耗时较多。

为此,本文从两个方面对代码进行了优化:一是算法本身的约减,二是算法基于TMS320C6678平台的优化。

2.3 算法优化

(1)算法约减

计算过程的优化主要对语种识别系统中对数似然得分模块的算法做约减。对数似然得分过程就是利用已经训练好的各语种GMM模型对输入的语音特征进行似然得分的计算,语种数越多,则该模块的耗时越多。利用Top n的方法,对每个模型选取得分最高的10个高斯用来计算对数似然得分。由于区分性Model Pushing模型是由SVM训练得到的支持向量重构而来,而支持向量由GMM-UBM模型自适应得到,因此,区分性Model Pushing模型与GMM-UBM模型的各高斯分量之间有着很强的对应关系。

上述介绍说明,区分性Model Pushing模型与GMM-UBM模型有着很强的对应关系,可近似认为对同一个特征向量它们得分最大的高斯混元一致[7]。针对拥有512个高斯混元的GMM,似然得分的计算结果必定仅仅集中于很少的几个高斯混元,大部分的高斯混元得分都会非常小以致可以忽略。因此,考虑将得分小的高斯混元结果忽略不计,只计算得分大的高斯混元。鉴于f DW-NAP模块包含特征向量对GMM-UBM计算后验概率的部分,可利用该部分的结果选取Top 10的高斯混元用于后端对数似然得分的计算。

(2)基于TMS320C6678平台的算法优化

基于平台的优化主要是通过选择CCSv5提供的编译优化参数来实现。通过不断的参数选择、搭配,获得最理想的参数优化方式,提高代码中循环运算的性能,使用软件流水调度技术提高代码的并行执行效率。

除此之外,特征提取阶段的FFT和f DWNAP的矩阵运算等算法采用DSPlib中优化的库函数进行替代,利用优化的库函数可以极大地提升代码的运行速度。

(3)算法优化前后识别性能对比

首先检验Top 10算法对系统识别性能的影响。在测试集中模型使用Top 10的区分性Model Pushing,前端特征参数保持不变,在VC++2010的环境下测试系统性能。实验所用语料库为实验室采集的电话信道通话语音,含汉语普通话、日语和英语3个语种,测试集包含汉语1 000段、日语450段及英语750段,共2 200段30 s的语音和3 000段10 s的语音(各语种1 000段)。系统性能用等错误率EER(Equal Error Rate)[2]衡量,实验结果如表2所示。

由表2可以看出,相对于全高斯得分模型,Top 10得分模型系统性能有所下降,主要因为舍弃了其他得分低的高斯成分,而其中必定包含部分语种区分信息,但舍弃掉的这一部分所含的语种信息有限,所以性能下降在可接受范围之内(相对下降小于5%)。该优化方法下模块的运算量下降是显而易见的,同样耗时也会大幅下降。

(%)

(4)算法优化前后系统实时性对比

对经算法优化的系统耗时做如下测试,同样以30 s的语音(8 000 Hz采样,帧长25 ms,帧移10 ms)为例,用CCSv5的环境进行软件仿真得到各模块处理所花的时钟周期数,然后按照TMS320C6678芯片的单个内核的工作主频(1.25 GHz)计算得到处理时间,结果如表3所示。

由以上分析可以看到,算法优化后的系统耗时由22.3 s减少至1.36 s,下降非常明显,其中下降最多的是f DWNAP模块和对数似然得分模块。在整个系统中,经过算法优化,f DWNAP模块耗时所占比例依旧最大,因此在多核任务并行设计时,需要将该模块的任务进行分解。

3 基于TMS320C6678的语种识别算法并行设计

3.1 模块间通信分析

根据语种识别的系统结构,测试过程分为3个模块,各模块的算法都已经进行了相应的优化。这些模块相互配合,通过控制信号完成数据流的交互。任务的控制流程主要是模块的执行次序,任务分配在不同核上的模块之间以传递消息的方式实现同步。模块间数据的传递会造成相应的时间延迟,因此,控制流程的设计准则为最大化系统的处理能力。模块间的数据流程主要是数据的传输方向,描述模块与外部数据间的相互关系。相反,最小化模块间的数据通信量则是数据流程的设计准则。

语种识别系统算法各模块间控制流程和数据流程的通信示意图如图2所示。该图由数据层和控制层两部分构成,控制信号的传输由虚线箭头表示,数据的传输由实线箭头表示。

3.2 模块任务的核映射

为了充分利用所有内核CPU的计算资源以最大限度地提高系统处理速度,根据算法优化前后的系统实时性测试结果及各模块运算量分析,将f DWNAP模块的矩阵乘法任务分配到多个核并行执行。

因本文的语种识别系统适合于数据流模式的任务并行方式,将整个系统的运算任务适当地分配给各个内核,实现任务级的并行流水。由于f DWNAP模块计算复杂度大,制约了整个系统任务级流水的处理速度。为了充分发挥TMS320C6678的性能优势,将该模块任务映射到多个核进行处理。该模块首先需要计算特征矢量对应的自适应GSV;然后通过投影矩阵计算SVM特征域的干扰空间,这一部分的大矩阵乘法占据了整个模块的绝大部分运算量;最后还需要将干扰空间返回映射到特征域,并在特征域去除干扰。整个模块80%以上的运算量都集中在大矩阵的乘法上,故采用将大矩阵拆为小矩阵分配到多个核上并行运算,将其他任务集中在一个核上进行处理。在该模块内还是一个任务级的流水处理方式,矩阵相乘部分是核级相同的并行流水处理方式。

4 基于TMS320C6678的语种识别算法实现

4.1 语种识别算法在TMS320C6678中的实现

根据设计思路,将本文提出的语种识别算法在CCSv5上进行软件仿真。其中,利用SYS/BIOS[8]提供核间任务调度,利用IPC[9]实现核间同步和通信。

启动系统,完成所有核的初始化后,首先调用IPC_start函数让各核进入同步等待状态,然后各核上的程序才能开始执行。从共享存储器划出MSM_IN和MSM_OUT 2块存储区,MSM_IN存储K-L变换矩阵和各语种GMM模型,MSM_OUT存储判决输出结果。投影矩阵数据存储在外接DDR3存储器中的位置信息事先存在Core1中。Core1将投影矩阵数据分成5份,通过Notify_send Event函数将5份数据的地址发送到Core2、Core3、Core4、Core5和Core6。Core2、Core3、Core4、Core5和Core6上的子矩阵乘法任务一直处于悬挂状态,直到Core1发送过来数据地址,矩阵乘法任务才开始并行执行。各核分别根据数据地址从外接DDR3读取数据与Core1传递的数据计算干扰因子向量,计算完毕再利用Message Q_put函数将干扰因子向量数据的Message写入到Core1建立的消息队列上。Core1利用Message Q_get函数从消息队列读取Message,从Message中获取干扰因子向量数据;然后计算补偿后的特征向量;接着Core1利用Message Q_put函数将补偿的特征向量数据的Message写入到Core7建立的消息队列上,Core7上的判决任务开始执行,最后将执行结果的数据写入MSM_OUT。

4.2 实验及结果分析

根据本文语种识别算法的TMS320C6678任务并行设计方案,本节将给出CCSv5平台下浮点算法的软件仿真结果,并进行分析验证。

按照3.2节的描述,将f DWNAP模块设计为并行处理,同样以30 s的语音为例,采用3个语种的模型测试整个系统在TMS320C6678上的实时性能。3个部分的运算处理时间结果如表4所示。

由表4可以看出,三个模块中f DWNAP模块耗时(0.227 s)最多,因此估算该系统的实时倍率至少为132(30/0.227)。

为了验证基于TMS320C6678平台的语种识别系统性能,将采用Top 10优化后的算法与在VC++2010平台中的识别性能进行对比。实验语料保持不变,表5给出了基于两种不同平台的系统EER。

实验结果表明,基于TMS320C6678平台的浮点软件仿真结果和VC++2010平台下的浮点计算结果完全一致,从而验证了TMS320C6678平台实现语种识别系统的正确性。

本文针对语种识别系统的实时性需求,在分析语种识别算法原理和多核DSP任务并行的基础上,分析了系统各模块的运算量,根据各模块的运算量对算法进行了优化。针对优化后算法的特点,设计了基于TMS320C6678平台的语种识别系统。最后从实时性和识别性能两个方面对系统性能进行了测试,结果验证了算法在TMS320C6678中实现的正确性及优化的有效性。

参考文献

[1]WONG K Y E.Automatic spoken language identificationutilizing acoustic and phonetic speech information[D].Queensland:Queensland University of Technology,2004.

[2]徐婷婷.语种识别中的若干问题研究[D].北京:北京邮电大学,2011.

[3]Texas Instrument.TMS320C6678 multicore fixed and floating-point digital signal processor[R].SPRS691C,2012.

[4]刘伟伟,吉立新,李邵梅.基于区分性Model Pushing的语种识别方法[J].电子技术应用,2012,38(4):113-116.

[5]刘伟伟,吉立新,李邵梅.基于区分加权干扰属性投影的语种识别方法[J].中文信息学报,[已录用未发表].

[6]Liu Weiwei,Ji Lixin,Li Shaomei.Robust cepstral featurecompensation for language recognition[C].Guangzhou:Proc.ofBEMI,2012:119-122.

[7]徐颖.语种识别声学建模方法研究[D].合肥:中国科学技术大学,2011.

[8]Texas Instrument.SYS/BIOS inter-processor communication(IPC)and I/O user’s guide[R].SPRUGO6C,2010.

上一篇:走过,才明白中考作文下一篇:中秋节做月饼diy活动方案