郑崇伟
Chongwei ZHENG
中图分类号:
P743
收稿日期:
2016-04-5
修回日期:
2016-06-5
网络出版日期:
2017-07-13
版权声明:
2017 《北京大学学报(自然科学版)》编辑部 《北京大学学报(自然科学版)》编辑部 所有
基金资助:
展开
摘要 利用欧洲中期天气预报中心的ERA-Interim海
关键词:
Abstract Based on the 36-year (1979-2014) ERA-Interim 10 m sea surface wind data from European Centre for Medium-Range Weather Forecasts (ECMWF), the climatic long term trend of the wind energy resource of the Gwadar Port of Pakistan was analyzed. Using two linear regression and artificial neural network (ANN) techniques, the wind energy resource in the long term was predicted. The results show that wind energy in summer is richer than that in winter. The stability in summer is better than that in winter. For the past 36 years, the wind power density, occurrence of effective wind speed and energy level occurrence have significant annual decreasing trends, of -0.78 W/(m2 · a), -0.21%/a, -0.22%/a separately. These trends mainly exhibit in summer, while no significant variation in winter. The stability (coefficient of variation, monthly variability index and seasonal variability index) does not have a significant long term trend for the past 36 years. From the prediction value, the wind energy resource in 2015 is similar to the multi-year average value, while the wind energy resource in 2016 is richer than the multi-year average value. For the year 2015-2016, the prediction wind energy will be more unstable than the multi-year average status. The results can provide scientific reference for the 21st Maritime Silk Road construction, development of remote islands and ports in the China seas.
Keywords:
PDF (52842KB)
元数据
多维度评价
相关文章
收藏文章
作为“21 世纪海上丝绸之路”的关键节点之一, 巴基斯坦瓜达尔港(Gwadar Port)的风能资源开发利用将有利于提高其可持续发展能力, 同时可以造福沿线人民, 成为我国的海外示范工程。资源开发应评估先行。郑崇伟等[1]利用来自 ECMWF 的 ERA-interim的风场资料, 率先评估瓜达尔港风能资源的季节特征, 结果表明: 瓜达尔港的风能资源常年可用, 年平均风能密度为 121 W/m2, 年平均有效风速频率为 43%; 风能常年稳定地由偏西南向的风贡献。Chang 等[2]对台湾地区的风速特征和风能资源进行分析, 发现该地区蕴藏着较为丰富的风能资源。陈飞等[3]利用30多年的气象资料, 分析了连云港地区及其近海的风能资源特征, 发现该地区风能储量丰富, 风速变化稳定, 有广阔的风能开发前景。毛慧琴等[4]利用 100 多个观测站数据, 对广东省沿海的风能进行分析, 发现广东省风能的潜在开发区分布于沿海和河口一带, 可开发量理论值为7.99×l03 MW。周荣卫等[5]指出我国沿海风能资源非常丰富, 沿海多年平均风功率密度约在 300~800 W/m2之间。郑崇伟等[6,7]对我国南海重点岛礁的海浪发电和海上风电可行性展开论证, 为南海岛礁建设提供科学依据。前人对我国近海的风能资源评估做了大量研究, 但目前关于瓜达尔港风能的研究还较为稀少, 针对其历史变化趋势和预测的研究更是凤毛麟角。
在风力发电和风力提水等风能的实际开发过程中, 往往关注资源的历史变化趋势以及未来预测, 以便为长期计划提供科学依据。本研究利用来自ECMWF 的近 36 年ERA-Interim海表 10 m 风场资料, 计算瓜达尔港的风能资源在近 36 年期间的历史变化趋势, 并利用线性回归和 BP 神经网络两种方法, 对该港的风能资源进行长期年度预测(30 分钟内的预测为超短期预测, 72 小时内的预测为短期预测, 10 天以上的预测为长期预测[8]), 为风能开发的长期规划、战略支撑点建设等提供科学依据, 为“21世纪海上丝绸之路”建设提供决策辅助。
利用欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts, ECMWF)近36年(1979—2014年)的ERA-Interim海表 10 m 风场资料。ERA-Interim 再分析资料是继其早期产品ERA-40 之后的新产品, 使用分辨率更高的气象模式, 在观测资料的应用及同化方法方面也有很大改进[9,10]。ERA-Interim 使用最新的 12 小时窗口四维变分同化技术, 同化的资料包括卫星资料、常规观测资料以及模式数据, 其作用是提供 ECMWF 早期产品与新一代产品之间的衔接, 目的是对 ERA-40 以及更早的数据进行完善。ERA-Interim 风场的空 间范围为 90°S—90°N, 180°W—180°E; 空间分辨率有0.125°×0.125°, 0.25°×0.25°, 0.5°×0.5°, 0.75°×0.75°, 1.0°×1.0°, ..., 2.5°×2.5° (本文选择的空间分辨率为0.125°×0.125°); 时间序列为1979年1月1日00:00至今(不断更新中)(本文选择的时间序列为 1979—2014 年); 时间分辨率为逐 6 小时(每 6 小时一个数据)[11]。Dee等[9]、Song等[10]、Bao 等[12]及马永锋 等[13]曾将 NCEP-CFSR, NCEP-NCAR, ERA-Interim 和 ERA-40 Reanalysis这几种数据与观测资料进行比较, 发现 ERA-Interim 在均方根误差和偏差方面更优。
利用 ERA-Interim 海表 10 m 风场资料, 根据风能密度的计算方法[14,15], 首先计算得到近 36 年逐 6 小时的瓜达尔港海表 10 m 的风能密度数据。基于该数据, 采用一元线性回归方法, 计算分析瓜达尔港风能资源的历史变化趋势, 主要包括风能密度、有效风速频率、能级频率、大风频率、稳定性在近 36 年的历史变化趋势(含 1 月、7 月和逐年的变化趋势)。此外, 还利用BP神经网络和线性回归预测方法, 对瓜达尔港的风能资源进行长期年度预测, 主要包括风能密度、有效风速频率、能级频率和稳定性的预测。目前线性回归和 BP 神经网络预测方法比较成熟[16,17], 广泛运用于气候分析与预 测[18]、城市可持续发展[19]、海平面预测[20]、生态系统风险研究[21]、金融预报[22]等诸多领域。
线性回归模型: 假设预报量为 y, 预测因子向量为 x, 则方程形式为
BP 神经网络拓扑结构如图 1 所示(以 3 层为例)。图 2 是相应于图 1 的网络运行流程, 其中 b1和 b2 为相应神经元的阈值(可看做一个特殊的权值, 故在图 1 中未标识), 黑色实线及箭头表示数据间的映射关系和流动方向。输入层及隐含层的输入“1”可以叠加到输入向量中。输入向量为 X = (x1, x2, …, xi, …, xn)T, 隐层输出向量为Y = (y1, y2, …, yj, …, ym)T, L为激发函数, 输出层向量为 O = (o1, o2, …, ok, …, ol)T, 目标输出为 D = (d1, d2, …, dk, …, dl)T。输入层到隐层之间的权值向量为 vi (
将瓜达尔港 1979 年 1 月逐 6 小时的风能密度做平均, 得到该月的平均值, 同理得到1979年1月至2014 年 12 月逐月瓜达尔港的风能密度。利用逐月的数据计算瓜达尔港在 1 月(代表冬季)、7 月(代表夏季)以及逐年的风能密度变化趋势, 结果见图3。
图3 近36年1月、7月和逐年的风能密度 Fig. 3 Monthly wind power density in January, July and annual value for the past 36 years
1 月, 瓜达尔港的风能密度在 50~170 W/m2 之间较为剧烈地波动, 线性相关度|R|=0.01, 未通过显著性检验, 表明瓜达尔港1月的风能密度在近36年期间无显著性变化趋势。
7 月, 在强劲西南季风的影响下, 风能密度在90~210 W/m2 之间波动, 且数值明显高于 1 月。显然, 这是由于该区域冬季冷空气的强度没有夏季西南季风强劲。值得注意的是, 瓜达尔港的情况与我国南中国海有着很大差异, 南中国海的冷空气强度明显大于西南季风。风能密度的线性相关度 R= 0.33, 通过 95%的信度检验, 线性趋势为-0.89, 即在近 36 年瓜达尔港 7 月的的风能密度以-0.89 W/(m2 · a)的速度显著性递减。从整体上看, 瓜达尔港 7 月的风能密度特征分为两个时段: 1979—1990年在 150 W/m2 上下波动, 1991—2014 年在 120 W/m2上下浮动。
逐年, 风能密度的线性相关度 R=0.57, 通过99.9%的信度检验, 线性趋势为-0.59, 即在近36年期间, 瓜达尔港的风能密度以-0.59 W/(m2·a)的速度显著性逐年递减。从整体上看, 年平均风能密度分为两个阶段: 1979—2000年在120 W/m2上下波动, 2001—2014年在110 W/m2上下波动。
通常风速在 5~25 m/s 之间有利于风能资源的采集与转换, 称为有效风速[23,24,25]。利用 1979 年 1月逐 6 小时的风速, 统计该月瓜达尔港有效风速出现的频率(简称有效风速频率), 用同样的方法得到1979 年 1 月至 2014 年 12 月逐月的有效风速频率。利用逐月的数据计算瓜达尔港在 1月、7 月以及逐年的有效风速频率变化趋势, 结果见图4。
图4 近36年1月、7月和逐年的有效风速频率 Fig. 4 Monthly occurrence of effective wind speed in January, July and annual value for the past 36 years
1 月, 有效风速频率在 30%上下浮动。线性相关度 R=0.07, 未通过显著性检验, 表明瓜达尔港 1月的有效风速频率在近 36 年期间无显著性变化趋势。
7 月, 有效风速频率的线性相关度 R=0.37, 通过 95%的信度检验, 线性趋势为-0.25, 即近 36 年瓜达尔港 7 月的有效风速频率以-0.25%/a (这里的%是有效风速出现的频率, 不是变率)的速度显著性递减。从整体上看, 瓜达尔港7月的有效风速频率分为两个时段: 1979—1990 年在 50%上下波动, 1991—2014年在40%上下浮动。
逐年, 有效风速频率的线性相关度 R=0.60, 通过 99.9%的信度检验, 线性趋势为-0.17, 即在近36 年期间, 瓜达尔港的有效风速频率以-0.17%/a的速度显著性逐年递减。从整体上看, 年平均有效风速频率分为两个阶段: 1979—2000 年在 45%上下波动, 2001—2014年在40%上下波动。
研究表明, 虽然瓜达尔港 200 W/m2 以上丰富能流密度出现频率较低, 但 100 W/m2 以上可用能流密度出现频率较为乐观[1]。在此, 统计分析瓜达尔港 100 W/m2以上能流密度出现的频率(简称 100 W/m2以上能级频率)的历史变化趋势。利用 1979年 1 月逐 6 小时的风能密度, 统计该月瓜达尔港100 W/m2以上能级频率, 用同样的方法得到 1979年 1 月至 2014 年 12 月逐月的 100 W/m2以上能级频率。利用逐月的数据, 计算瓜达尔港在 1 月、7月以及逐年的 100 W/m2以上能级频率的变化趋势, 结果见图5。
图5 近36年1月、7月和逐年的100 W/m2以上能级频率 Fig. 5 Monthly occurrence of wind power density greater than 100 W/m2 in January, July and annual value for the past 36 years
1 月, 100 W/m2 以上能级频率在 30%上下浮动。线性相关度 R=0.08, 未通过显著性检验, 表明瓜达尔港的100 W/m2以上能级频率在近36年期间的1月无显著性变化趋势。
7月, 线性相关度R=0.37, 通过95%的信度检验, 线性趋势为-0.22, 即在近 36 年的 7 月, 瓜达尔港的 100 W/m2以上能级频率以-0.22%/a (这里的%是 100 W/m2以上能级频率, 不是变率)的速度显著性递减。从整体上看, 瓜达尔港 7 月的 100 W/m2以上能级频率分为两个时段: 1979—1990 年在45%上下波动, 1991—2014年在40%上下浮动。
逐年, 线性相关度 R=0.64, 通过99.9%的信度检验, 线性趋势为-0.18, 即在近 36 年期间, 瓜达尔港 100 W/m2以上能级频率以-0.18%/a 的速度显著性逐年递减。从整体上看, 年平均 100 W/m2以上能级频率在30%~45%之间波动。
风力等级频率是风能开发和防灾减灾中较受关注的要素[26]。研究表明, 瓜达尔港出现频率最高风力为 3 级(34.29%), 其次是 4 级(28.32%)和 2 级(23.11%), 6 级以上大风频率很低, 常年在 1.5%以内[1]。我们统计得到 1979 年 1 月至 2014 年 12 月逐月的风力等级频率, 通过计算, 发现瓜达尔港 1月、7 月以及逐年的 6 级以上大风频率无特别显著的变化趋势(图略)。
风能资源的稳定性直接决定着风机对资源的采集和转换能力以及风机的寿命。Cornett[27]曾计算全球海域波浪能的变异系数(coefficient of variation, Cv)、月变化指数(monthly variability index, Mv)和季节变化指数(seasonal variability index, Sv)来判断其稳定性。本文参考他的计算方法, 利用近36年逐6 小时的风能密度, 计算瓜达尔港风能密度的Cv, Mv和Sv来判断其稳定性, 结果见图6和7。
图6 近36年1月、7月和逐年的变异系数 Fig. 6 Monthly coefficient of variation in January,July and annual value for the past 36 years
图7 风能密度在近36年逐年的月变化指数(a)和季节变化指数(b) Fig. 7 Monthly variability index (a) and seasonalvariability index (b) for the past 36 years
变异系数: 1月在1.4上下浮动, 且无显著的长期历史变化趋势; 7月在1.0上下浮动, 也不存在显著的变化趋势, 在1995年存在一个突变点。对比1月和 7 月的变异系数不难发现, 7 月的数值整体上偏小, 即 7 月的稳定性比1月好。变异系数的年平均值在1.15上下波动, 无显著的逐年变化趋势。
月变化指数: 1979—1994 年期间在 14 上下浮动。近36年期间, 瓜达尔港风能的月变化指数无显著的变化趋势, 但在 1996—2014 年期间却以0.0394的速度逐年显著性递增。
季节变化指数: 线性相关度 R=0.04, 未通过显著性检验, 表明瓜达尔港风能的季节变化指数在近 36 年期间无显著性变化趋势, 数值在 0.7 上下 浮动。
为了更真实地反映风能各要素的变化趋势, 我们对时间序列进行 5 点滑动, 分析风能密度、有效风速频率、100 W/m2以上能级频率、变异系数、月变化指数和季节变化指数的历史变化趋势, 结果见表 1。可以看出, 瓜达尔港风能各要素在 1 月都不存在显著的变化趋势。7 月, 风能密度、有效风速频率以及 100 W/m2以上能级频率都表现出显著的递减趋势, 通过 95%的信度检验; 滑动平均后, 上述 3 个要素的递减趋势更为显著, 均通过 99.9%的信度检验。5 点滑动后, 风能密度、有效风速频率和100 W/m2以上能级频率分别以-0.78 W/(m2 · a), -0.21%/a 和-0.22%/a 的速度逐年递减, 均通过99.9%的信度检验; 变异系数、月变化指数和季节变化指数在近36年期间无显著的变化趋势。
表1 瓜达尔港风能各要素的历史变化趋势
Table 1 Climatic long-term trend of the wind energy parameters of the Gwadar Port
历史变化趋势分析可以较为直观地展现风能资源在过去几十年的变化。在海洋开发过程中, 往往更关注未来半年、一年或两年的风能资源情况, 以便为长期计划提供依据。短期预报和历史变化趋势分析都不能满足需求, 这就需要长期预测。本文利用线性回归和 BP 神经网络两种方法, 对瓜达尔港的风能资源展开年度预测。
利用 1979—2012 年逐月的风能密度、有效风速频率、能级频率和变异系数做预测训练, 对2013年 1 月至 2014 年 12 月逐月的数值进行预测。本文主要利用线性回归和 BP 神经网络两种方法进行预测, 结果见表 2。预测值与实际值之间的相关系数(correlation coefficient, CC)都通过 99.9%的信度检验(CC>0.60), 均方根误差(root mean square error, RMSE)、偏差(Bias)、平均绝对误差(mean absolute error, MAE)、离散系数(scatter index, SI)以及标 准均方根误差(normalized root mean square error, NRMSE)也在可以接受的范围内。从整体上看, BP神经网络的预测效果要优于线性回归预测, 基于人工神经网络的风能各要素的预测值与实际值之间的相关系数基本上都在 0.7 以上, 明显高于线性回归预测值。此外, 从RMSE, Bias, MAE, NRMSE和SI来看, 人工神经网络的误差均小于线性回归预测, 这在风能密度、有效风速频率、能级频率和变异系数等几个风能要素上都得到较好的体现。我们对比 3, 5, 7层BP神经网络, 发现5层效果更优(表略)。因此, 下面利用BP神经网络(5 层网络)对风能各要素展开预测。
表2 两种预测方法对瓜达尔港风能资源的预测训练效果
Table 2 Prediction precision of the wind energy of the Gwadar Port
为了更直观地展现预测效果, 我们绘制风能密度、变异系数的真实值和预测值曲线(图 8)。从整体上看, 线性回归和方法 BP 神经网络方法对风能密度和变异系数都有一定的预测能力, 但 BP 神经网络的预测值曲线在走势上与真实值保持更好的一致性。BP 神经网络对风能密度在 2013 年 5 月和2014 年5 月的两个波峰、变异系数在 2014 年 1 月的波峰体现得更明显。值得注意的是, 两种预测方法都未能很好地刻画风能密度在 2013 年 9 月的相对峰值和变异系数在2014 年 11 月的相对峰值。
图8 线性回归和人工神经网络对风能密度和变异系数的预测训练效果 Fig. 8 Prediction experiment of wind power density and coefficient of variation of the wind energy
利用线性回归和人工神经网络两种方法, 基于1979—2014 年逐月的风能密度、有效风速频率、能级频率和变异系数, 对上述各要素做 24 个月的长期年度预测, 并将预测值与多年平均值进行比较, 结果见图9。
图9 风能密度、有效风速频率、100 W/m2以上能级频率和变异系数未来24个月的预测效果((a)~(d))以及预测值与多年平均值的差值((e)~(h)) Fig. 9 Prediction values of wind power density, occurrence of effective wind speed, occurrence of wind power density greater than 100 W/m2 and coefficient of variation for the future 24 months ((a)-(d)), and the abnormal values between prediction and multi-year average value ((e)‒(h))
由图 9(a)~(c)可见, 对于风能密度、有效风速频率和100 W/m2以上能级频率, 2015年的预测值与多年平均值基本上接近, 2016年的预测值则明显高于多年平均状态, 尤其是风能密度。这意味着瓜达尔港的风能资源在 2015 年与多年平均状态持平, 2016 年则趋于更丰富。由图 9(d)可见, 变异系数的
预测值在 2015 年 5—11 月与多年平均值接近, 其余月份略高于多年平均值, 表明瓜达尔港 2016 年风能的稳定性比多年平均状态稍差。
为了更清晰地展现预测值与多年平均值的差异, 本文计算了预测值与多年平均值的差值, 结果见图 9(e)~(h)。从图 9(e)可见, 2015 年的预测值与多年平均值接近, 仅 2015 年 6 月明显低于多年平均值; 2016年的风能密度比多年平均状态丰富(高出 0~40 W/m2)。从图 9(f)可见, 有效风速频率在2015 年的预测值与多年平均值接近, 2016 年的预测值则明显高于多年平均值(高出 0~10%)。从图 9(g)可见, 100 W/m2 以上能级频率的预测值 2015 年与多年平均值基本上接近, 2016 年的预测值比多年平均值高出 0~9%。从图 9(h)可见, 变异系数 2015 年和 2016 年的预测值都略高于多年平均值(尤其是2016 年冬季), 表明预测的风能密度稳定性比多年平均的风能密度稳定性差。
本文利用来自 ECMWF 的 ERA-Interim 海表10 m 风场资料, 计算了瓜达尔港的风能资源在近36 年期间的历史变化趋势, 并利用线性回归和 BP神经网络两种预测方法, 对该港的风能资源展开长期年度预测, 结论如下。
1) 瓜达尔港夏季的风能密度为 90~210 W/m2, 明显比冬季(50~170 W/m2)丰富; 年平均风能密度为 100~140 W/m2。夏季的有效风速频率(45%左右)和 100 W/m2以上能级频率(40%左右)明显高于冬季(有效风速频率、100 W/m2以上能级频率都在 30%左右); 年平均有效风速频率为 35%~50%, 100 W/m2以上能级频率的年平均值为 30%~45%。风能资源的稳定性夏季好于冬季, 年平均变异系数、月变化指数和季节变化指数分别为 1.0~1.4, 0.6~2.4和0.3~1.3。
2) 冬季, 瓜达尔港的风能密度、有效风速频率、能级频率和变异系数无显著变化趋势。夏季, 风能密度、有效风速频率和能级频率分别以-0.93 W/(m2 · a)、-0.25%/a和-0.22 %/a的速度显著递减; 变异系数无显著变化趋势。近 36 年期间, 风能密度、有效风速频率和能级频率分别以-0.78 W/ (m2 · a)、-0.21%/a 和-0.22%/a 的速度逐年显著递减, 稳定性(变异系数、月变化指数、季节变化指数)无显著性逐年变化趋势。
3) 从预测值来看, 风能密度、有效风速频率和 100 W/m2以上能级频率的预测值 2015 与多年平均值基本上接近, 2016年的预测值则明显高于多年平均状态, 尤其是风能密度; 2015—2016 年, 预测风能资源的变异系数略高于多年平均值, 表明其稳定性比多年平均状态略差。
The authors have declared that no competing interests exist.
瓜达尔港的风能开发可行性论证 Assessment of wind characteristics and wind turbine characteristics in Taiwan 连云港沿海地区及近海风能资源评估 广东省风能资源区划研究 中国近海风能资源开发潜力数值模拟 中国南海岛礁建设: 重点岛礁的风候、波候特征分析 中国南海岛礁建设: 风力发电、海浪发电 风电功率预测技术综述 The ERA-Interim reanalysis: configuration and performance of the data assimilation system Comparison of wind data from ERA-Interim and buoys in the Yellow and East China Seas ERA-Interim wind data Evaluation of NCEP-CFSR,NCEP-NCAR,ERA-Interim,ERA-40 reanalysis datasets against independent sounding observations over the Tibetan Plateau ERA-Interim再分析和NCEP FNL分析资料在东南极中山站至 Dome A 断面的适用性研究 江苏环港风电场风能资源特征分析 江苏省风能资源重新估算与分布研究 一种基于人工神经网络的基本块重排方法 人工神经网络在星载散射计海面风场反演建模中的应用 近年来CCA在气候分析与气候预测中的应用 基于人工神经网络的城市拓展区可持续发展指数序列研究 基于长期验潮数据的海平面预测方法与案例应用 基于人工神经网络的青藏公路铁路沿线生态系统风险研究 人工神经网络及其在金融预报中的应用 Active power regulation of wind power systems through demand response Assessment of the global ocean wind energy resource 经略 21 世纪海上丝路: 海洋资源、相关国家的海洋开发状况 全国低风速气象特征分析 A Global wave energy resource assessment // Proceedings of the Eighteenth(
/
1 资料与方法
1.1 资料
1.2 方法





2 风能资源的历史变化趋势
2.1 风能密度
2.2 有效风速频率
2.3 能级频率
2.4 大风频率
2.5 稳定性
2.6 滑动平均后的趋势
风能要素 1月 7月 逐年 原始值 滑动平均 原始值 滑动平均 原始值 滑动平均 风能密度/(W · m-2 · a-1) -0.03 -0.05 -0.89* -0.93** -0.59** -0.78** 有效风速频率/(% · a-1) -0.07 -0.08 -0.25* -0.25** -0.17** -0.21** 能级频率/(% · a-1) -0.07 -0.09 -0.22* -0.22** -0.18** -0.22**
3 风能资源的中长期预测
3.1 预测训练
项目 风能密度 有效风速频率 能级频率 变异系数 线性回归 神经网络 线性回归 神经网络 线性回归 神经网络 线性回归 神经网络 CC 0.6593* 0.8608* 0.6724* 0.7331* 0.6634* 0.7002* 0.6352* 0.8484* Bias 5.7897 3.0560 1.7890 -0.2516 1.9375 -0.5617 0.1101 0.0955 MAE 33.4049 18.2077 9.9247 7.3124 9.6059 9.0033 0.2073 0.1494 NRMSE 0.3347 0.1811 0.2795 0.2318 0.3174 0.2856 0.2028 0.1458 RMSE 42.8393 23.1813 12.4186 10.0006 12.5472 11.7034 0.2452 0.1762 SI 0.3316 0.1811 0.2766 0.2617 0.3136 0.2768 0.1812 0.1225
3.2 各要素预测
4 结论
参考文献
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
技术支持:北京玛格泰克科技发展有限公司
〈
〉