1.清华大学 自动化系, 北京 100084;
2.华能新能源股份有限公司, 北京 100036
收稿日期: 2016-01-07
基金项目: 中国华能集团公司重大专题科技项目(HNKJ13-H20-01)
作者简介: 封宇(1992-), 男, 博士研究生。
通讯作者: 黄必清(1966-), 教授, E-mail: hbq@tsinghua.edu.cn
摘要：海洋蕴含的风能资源比陆地更加丰富。由于海上与陆地气象要素的差异, 在进行近海及海上风资源评估时, 通常需要考虑海上风资源时间特性和空间特性的影响。该文根据某公司海上测风塔和近海测风塔的实际测风资料, 研究了海上风资源在时间上的不同分布, 重点使用了EM(expectation-maximization)算法, 研究传统风资源评估方法几乎不考虑的风速在昼夜表现出的不同分布。同时, 基于实测数据, 应用机器学习算法, 结合Monin-Obukhov相似性理论、参数离散替代方案以及海面空气动力学粗糙度参数化方案, 对海上风资源的空间特性进行了分析研究, 有效地反映出海洋表面动力学粗糙度变化对于风速垂直变化的影响, 弥补了传统切变公式的不足。结果揭示了海上风资源的时空变化规律, 为近海及海上风资源评估和海上风电场的前期规划提供了依据。
关键词： 风能 海上风资源 时空特性 风资源评估 Monin-Obukhov相似性理论 空气动力学粗糙度
Temporal and spatial characteristics of offshore wind resources
FENG Yu1, HE Yan2, ZHU Qihao1, GUO Chen2, FENG Xiaodan2, HUANG Biqing1
1.Department of Automation, Tsinghua University, Beijing 100084, China;
2.Huaneng Renewables Corporation Limited, Beijing 100036, China
Abstract:Offshore wind energy resources are more abundant than on land. However, sea and land winds have different meteorological elements so offshore wind resource assessments need to take the impact of the temporal and spatial characteristics of the offshore wind resources into consideration. Data from an offshore wind measurement mast was used to study offshore wind distributions at different time scales. The EM (expectation-maximization) algorithm was used to study the differences in offshore wind distributions between day and night, which is normally not considered in traditional wind assessment methods. The spatial characteristics of offshore winds will than analyzing using a machine learning algorithm, Monin-Obukhov similarity theory, and a parameter replacement scheme in discrete calculations and in an ocean surface aerodynamic roughness model. This method efficiently reflects the impact of ocean surface aerodynamic roughness changes on the vertical variations of the wind speed that is not considered in traditional wind shear formula. The results show the temporal and spatial characteristics of the offshore wind resources as a basis for better offshore wind assessments for planning offshore wind farms.
Key words: wind energyoffshore wind resourcetemporal and spatial characteristicswind resource assessmentMonin-Obukhov similarity theorysurface aerodynamic roughness
由于海上风资源优于陆上，加上风机技术的快速发展，近年来海上风电的发展十分迅速。进行海上风电场的建设需要对拟建风电场位置处的风资源进行评估，评估不仅为风电场建设提供客观依据，更是风电场能否取得良好经济效益，达到预期收益的关键，是整个风电场开发建设的重要环节。
目前各国虽然已经存在相关的风资源评估标准，但大多是根据陆地上的环境特点制订的。海上风资源与陆上有一定的相似性，但是海上风资源还有其自身的特点，进行海上风资源评估还需要考虑气温、水温、潮位、昼夜等因素的影响。传统风资源评估方法几乎不考虑风速在昼夜表现出的不同分布，但是这种差别在风切变普遍较小的海上显得十分明显，考虑昼夜因素能够更加准确地评估风电机组的发电量。同时，由于海洋表面每时每刻都在发生变化，动力学粗糙度不再像陆地上一样保持恒定，传统的切变公式未能反映海洋环境的变化对于动力学粗糙度的影响，因此需要结合Monin-Obukhov相似性理论和海面空气动力学粗糙度参数化方案来描述海上风资源的垂直变化规律。本文根据某公司海上测风塔的实测数据，对海上风资源的时间分布特点和垂直空间分布规律进行了详细的分析研究，为海上风资源评估和海上风电开发提供了有益的参考。
1 海上风资源的时空特性研究方法1.1 风资源的时间分布研究方法1.1.1 传统方法传统的风资源评估主要研究风速在单日的变化趋势，在一年中不同月份的变化规律，以及一年中不同季节的变化规律等。风速的单日变化是指每日同一时刻风速特征值的变化。在气压梯度不大的情况下，风速在一天内将表现出规律性的变化趋势。通常来讲，风速的日变化主要与大气层下垫面的性质有关，也与海陆昼夜稳定度变化情况有关。风速的月变化是指一年中每月平均风速特征值的变化，一般与季节的更迭以及当地的气候状况有关。风速的季节变化是指一年中四季平均风速特征值的变化。赵燕华等[1]基于统计方法对风资源在不同季节中单日的变化进行了研究，并比较了内陆和沿海变化规律的差异。符平等[2]根据测风塔的长期数据，研究了海上风资源的月变化和年际变化特征。李常春[3]则使用了瑞利分布结合NASA(National Aeronautics and Space Administration)数据进行了风资源的评估及发电量的估算。传统方法几乎不考虑风速在昼夜表现出的不同分布，但是昼夜分布差异在风切变普遍较小的海上显得十分明显。
1.1.2 方法改进风速的昼夜变化是指白昼和黑夜风速特征值的变化，在传统方法中经常被忽略，但是当风能评估与电力负荷相联系时，忽略昼夜的变化很可能导致评估出现较大误差。特定地区的风资源在白昼与黑夜通常表现出不同的分布，2个分布之间也可能存在着内在的联系。本文使用双参数混合Weibull分布描述风速，并结合测风数据使用EM算法进行了参数估计，得到昼夜2个不同分布的模型。
在概率统计中，如果某个随机变量X的概率密度函数可以表示为多个概率密度函数的和，即：
| $\eqalign{ & f\left( {X = x} \right) = \sum\limits_{i = 1}^k {{p_i}} {f_i}\left( {x;{\theta _i}} \right) = \cr & {p_1}{f_1}\left( {x;{\theta _1}} \right) + \cdots + {p_k}{f_k}\left( {x;{\theta _k}} \right) \cr} $ | (1) |
本文假设风速V服从由2个成分组成的混合分布，其中每个成分均服从双参数Weibull分布。即假定：k=2，θ
| $\eqalign{ & f\left( {V = v;\phi } \right) = \cr & {p_1}{f_1}\left( {\upsilon ;{k_1},{c_1}} \right) + {p_2}{f_2}\left( {\upsilon ;{k_2},{c_2}} \right) \cr} $ | (2) |
下面的问题就是如何根据测风数据对该混合分布进行参数估计。由于本文具有一系列从风速混合分布模型中采样得到的风速样本，因此可以使用最大似然法进行模型参数的估计。风速混合分布的似然函数为
| $\eqalign{ & L\left( {\upsilon ;{p_1},{p_2},{k_1},{c_1},{k_2},{c_2}} \right) = \cr & \mathop \prod \limits_{i = 1}^n \left\{ {{p_1}\left[ {\frac{{{k_1}}}{{{c_1}}}{{\left( {\frac{{{\upsilon _i}}}{{{c_1}}}} \right)}^{{k_1} - 1}}\exp \left[ { - {{\left( {\frac{{{\upsilon _i}}}{{{c_1}}}} \right)}^{{k_1}}}} \right]} \right]} \right. + \cr & \left. {{p_2}\left[ {\frac{{{k_2}}}{{{c_2}}}{{\left( {\frac{{{\upsilon _i}}}{{{c_2}}}} \right)}^{{k_2} - 1}}\exp \left[ { - {{\left( {\frac{{{\upsilon _i}}}{{{c_2}}}} \right)}^{{k_2}}}} \right]} \right]} \right\}. \cr} $ | (3) |
1) 人工设定混合分布参数的初值p1、p2、k1、k2、c1、c2.
2) E-step：
假设隐变量为w1=p1，w2=p2分别代表2类成分的权重。对于每个数据vi，估计该数据由第j(j=1,2)个成分生成的概率
| ${\gamma _{ij}} = \frac{{{\omega _j}f\left( {{\upsilon _i};{k_j},{c_j}} \right)}}{{\sum\limits_j {{\omega _j}f\left( {{\upsilon _i};{k_j},{c_j}} \right)} }}$ | (4) |
优化每个分量的参数
| $\mathop {\arg \max }\limits_{{k_j}{c_j}} \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^2 {{\gamma _{ij}}} } \log \frac{{{\omega _j}f\left( {{\upsilon _i};{k_j},{c_j}} \right)}}{{{\gamma _{ij}}}}$ | (5) |
| $P\left( {\upsilon \leqslant {\upsilon _g}} \right) = 1 - \exp \left[ { - {{\left( {\frac{{{\upsilon _g}}}{c}} \right)}^k}} \right].$ | (6) |
| $ln\left\{ { - \ln \left[ {1 - P\left( {\upsilon \leqslant {\upsilon _g}} \right)} \right]} \right\} = k\ln {\upsilon _g} - k\ln c.$ |
| $\eqalign{ & y = ln\left\{ { - \ln \left[ {1 - P\left( {\upsilon \leqslant {\upsilon _g}} \right)} \right]} \right\}, \cr & x = ln{\upsilon _g}, \cr & b = - k\ln c, \cr & a = k \cr} $ |
| ${\text{y = ax + b}}$ |
| $\eqalign{ & k = a, \cr & c = \exp \left( { - \frac{b}{a}} \right) \cr} $ |
| $\eqalign{ & W = \cr & \int_{\upsilon = {\upsilon _{in}}}^{{\upsilon _{out}}} {\frac{k}{c}} {\left( {\frac{\upsilon }{c}} \right)^{k - 1}}\exp \left[ { - {{\left( {\frac{\upsilon }{c}} \right)}^k}} \right]f\left( \upsilon \right)d\upsilon \bullet t. \cr} $ | (7) |
1.2.2 方法改进根据Monin-Obukhov相似性理论，大气边界层平均风廓线可表示为对数—线性律[6]：
| $u\left( z \right) = \frac{{{u^*}}}{\kappa }\left[ {ln\left( {\frac{z}{{{z_0}}}} \right) - \varphi \left( {\frac{z}{L}} \right)} \right].$ | (8) |
为了提高Monin-Obukhov在大气处于稳定状态时的预测精度，Gryning等[7]在公式中加入了修正项：
| $u\left( z \right) = \frac{{{u^*}}}{\kappa }\left[ {ln\left( {\frac{z}{{{z_0}}}} \right) - \varphi \left( {\frac{z}{L}} \right) + \frac{z}{{{L_{MBL}}}}} \right].$ | (9) |
| ${L_{MBL}}{\text{ = }}\frac{{{u^*}}}{f}/\left[ { - 2\ln \left( {\frac{{{u^*}}}{{f{z_0}}}} \right) + 55} \right].$ | (10) |
| $L = A{\left( {{u^*}} \right)^2},A = 1.1 \times {10^3}{{\text{s}}^2}{m^{ - 1}}$ |
当大气处于中性稳定状态时，公式中的线性项(Monin-Obukhov普适函数)将为0，可以直接从公式中剔除；当大气处于稳定状态时，Monin-Obukhov普适函数将与高度和Monin-Obukhov长度的比值成正比；当大气处于不稳定状态时，则该项的形式更加复杂。φz/L在学术界普遍认为有如下形式
| $\varphi \left( {{z \over L}} \right) = \left\{ \matrix{ 2\ln \left( {{{1 + {\Phi _m}} \over 2}} \right) + \ln \left( {{{1 + {\Phi _m}^2} \over 2}} \right) - \hfill \cr 2\arctan \left( {{\Phi _m}} \right) + {\pi \over 2},\left( {{z \over L}} \right) ＜ 0; \hfill \cr - \beta \left( {{z \over L}} \right),\left( {{z \over L}} \right) > 0. \hfill \cr} \right.$ | (11) |
在陆地上，通常人们将表面动力学粗糙度长度看作常数，但是由于海面状况时刻都在发生变化，因而海洋表面动力学粗糙度是一个难以确定的量。为了方便处理，许多****都对海洋表面动力学粗糙度的参数化方案进行了研究。Charnock模型[11]认为海洋表面动力学粗糙度只与摩擦风速有关：
| ${z_0} = {z_{ch}}\frac{{{{\left( {{u^*}} \right)}^2}}}{g}.$ |
Smith[12]则在Charnock模型基础上加入了海面光滑流的影响：
| ${z_0} = {z_{ch}}\frac{{{{\left( {{u^*}} \right)}^2}}}{g} + 0.11\frac{\upsilon }{{{u^*}}}.$ |
| $\eqalign{ & Z = {[{z_1}{z_2} \cdots {z_n}]^T}, \cr & {\text{V}} = {{\text{[}}{\upsilon _1}{\upsilon _2} \cdots {\upsilon _n}]^T}, \cr} $ |
J(Z,V;u*)=$\sum\limits_{i = 1}^m {} $[uzi;u*－vi]2.
根据实际情况调整学习速度α，那么学习过程可以表示为
u*:=u*－α·$\frac{{\partial J\left( {Z,V;{u^*}} \right)}}{{\partial {u^*}}}$.
直到算法达到规定的迭代步数或者u*的变化很小。
使用机器学习方法能够有效地避免脏数据的出现对于风廓线拟合产生的影响，增强模型预测不同高度风速时对于异常输入的适应性。
1.3 海上风资源评估方法本文结合了海上风资源的时间特性与空间特性，提出具有海上特色的风资源评估方法，如图1所示。
 |
| 图 1 海上风资源评估方法 |
| 图选项 |
首先，基于实际测风资料得到风资源在不同粒度的时间序列，并根据序列通过拟合双参数Weibull分布得到风速的概率分布模型。其次，获得测风位置的海洋参数，并通过海洋粗糙度参数化方案结合相似性理论得到垂直的风速廓线，并推算出轮毂高度的风速大小。结合风电机组在白昼和黑夜不同时间的可利用率数据，最终推算出风电机组的发电量。并结合白昼和黑夜的不同电价，计算风电场的收益。
2 海上测风数据分析2.1 海上测风塔信息测风塔的地理位置在北纬39°13′(N 39° 13′)，东经118°20′(E 118° 20′)，位于我国渤海西部渤海湾处，西部与天津的陆岸相邻，南部与渤海相通，有海岸线23.5 km。数据持续时间为2014年1月1日—2014年12月31日，采样周期为10 min。数据包括每个通道内每10 min的平均值、最大值、最小值和方差。
2.2 数据的验证在处理测风数据时，必须从数据的完整性、连贯性和合理性着手，对提取的测风数据进行检查，挑选出不合理的、可疑的以及漏测的数据，对其进行适当的修补处理，从而整理出较实际合理的完整数据以供进一步分析处理。
根据国家标准GB/T 18710—2002对风电场测风原始数据制定的合理性标准，对不合理的数据的处理方法如下：
1) 首先使用其他测风塔同期的相同高度实测数据进行补充；
2) 同一层高度所测数据均不合理或缺失时，由相邻高度的本测风塔测风数据根据风切变指数进行插值；
3) 同期本测风塔所有数据均不合理或缺失时，对与本测风塔最近的测风塔进行相关性分析，选取相关性好的测风数据进行补充。
经过对不合理数据的订正后，测风塔数据完整率的计算结果如表1所示。
表 1 测风数据的完整率计算
| 变量名称 | 应测数据量 | 缺测数据量 | 无效数据量 | 数据完整率 |
| 计算结果 | 210 240 | 58 | 806 | 99.6% |
表选项
2.3 海上风资源的时间分布1) 单日变化。
风速的单日变化如图2所示。图2中各条曲线由上到下分别代表80、50、30和10 m高度的测风数据一年内同一时刻的平均值，可以看出，近海及海上区域风速相对于陆地区域较大，且呈现明显的峰谷特性，其峰值出现在夜间，谷值出现在清晨，白天风速却相对较小。
 |
| 图 2 风速的单日变化曲线 |
| 图选项 |
2) 月变化。
对2014年1月—2014年12月一年的数据进行整理后，测风塔70 m的逐月平均风速变化如图3所示。从风速的月际分布来看，风能资源在不同月份呈现不同的分布，风能资源最好的月份是5月和12月，风能资源最差的月份是8月和9月。从风速的季节分布来看，春季风能资源最好，冬季次之，夏季风能资源最差。
 |
| 图 3 风速的月变化曲线 |
| 图选项 |
3) 昼夜变化。
针对风速的昼夜变化做研究，对2014年1月—2014年12月一年的数据进行整理后，以6：00—18：00为白昼，18：00—次日6：00为黑夜，各测风塔80 m高度的逐月平均昼夜风速变化如图4所示。
 |
| 图 4 风速的昼夜变化曲线 |
| 图选项 |
从图4中可以看出，除了6月份的夜晚平均风速与白天平均风速相差不大之外，其余的月份均呈现出明显的夜晚平均风速高于白天平均风速的情况。在80 m高度上，一年中每个月黑夜的平均风速高于白昼风速，黑夜中的风资源高于白昼。在其余的高度上亦可获得相似的结论。由于海水的热容量大，因此夜间易形成海水温度高于气温的现象，使得海面大气呈不稳定层结，有利于空气对流发展和高空动量下传，从而造成海上夜间风速较大的特点；白天则恰恰相反，常会有海水温度低于气温的情况，使得海面大气层结稳定，造成风速略小或变化不大的现象。
根据昼夜划分，通过最大化似然函数，使用EM方法得到的Weibull分布参数如表2所示。
表 2 海上风资源昼夜分布参数
| 高度/m | 白昼 | 黑夜 | ||
| 尺度参数c | 形状参数k | 尺度参数c | 形状参数k | |
| 80 | 6.625 | 1.894 | 7.068 | 1.998 |
| 50 | 6.526 | 1.897 | 6.792 | 2.009 |
| 30 | 6.301 | 1.949 | 6.441 | 2.032 |
| 10 | 5.654 | 1.927 | 5.353 | 1.900 |
表选项
从得到的Weibull分布c参数可以看出，夜晚风速比白天风速普遍要高，符合之前统计得出的结论。通过白昼和黑夜计算得到的不同Weibull分布参数值进行发电量的预测，能够获得更加精确的结果。
2.4 海上风资源的空间分布本文使用测风塔的实测数据进行训练，针对替代与简化后Monin-Obukhov方程进行参数估计，对海上风速廓线进行拟合。使用传统Charnock模型得到的结果如图5所示，使用Smith模型得到的结果如图6所示。
 |
| 图 5 Charnock模型对于风速垂直变化的拟合结果 |
| 图选项 |
 |
| 图 6 Smith模型对于风速垂直变化的拟合结果 |
| 图选项 |
从图5(a)，5(b)和图6(a)，6(b)可以看出，近海及海上风速廓线在高度较低的时候近似为一直线，近海及海上的风切变系数较小，风速随高度的变化不显著，并且Charnock模型对于风速垂直变化的估计普遍偏高，Smith模型对于风速廓线的拟合较好，更适合于描述海面风速的垂直变化。使用Monin-Obukhov模型结合海平面粗糙度模型进行风速廓线拟合得到的结果达到了不错的水平。从图5(c)，5(d)和图6(c)，6(d)可以看出，当测风塔80 m高度处出现异常数据时，Monin-Obukhov模型则能够有效地去除由于测风数据异常对于模拟海面风速廓线的影响，准确地保留风速变化趋势，更加准确地模拟稳定状态下其余高度的风速值。
针对全体测风数据，使用不同的方法对海上风速廓线进行拟合，得到的均方根误差如表3所示。
表 3 不同风廓线拟合方法的均方根误差
| 拟合方法 | 均方根误差RMSE/(m·s-1) |
| 风切变拟合 | 1.7 |
| Charnock-Monin-Obukhov | 2.2 |
| Smith-Monin-Obukhov | 1.5 |
表选项
结合表3和图6可以看出，传统的风切变拟合方法对于海上数据的拟合效果虽然误差较小，但是对于出现逆向切变的情况不能进行有效地拟合，受脏数据的影响较大，对于异常数据的敏感性强。而Monin-Obukhov理论能够有效地减小异常数据对于理论预测的影响，取得不错的拟合结果。
3 结论根据海上测风塔的实测风资料，对海上风资源的时空特性进行了研究与综合分析，可为海上风资源评估和海上风电场建设提供参考：1) 海上风资源在昼夜上明显体现出不同的分布，进行发电量预估时应考虑昼夜差异。2) 海上风速沿垂直高度方向变化较小，切变系数较小，使用Monin-Obukhov理论能更好地解释风速在垂直高度的变化，并且对于异常数据的适应性更好。3) 使用Monin-Obukhov 长度在离散高度上的近似替代方案能够有效地描述海上垂直风速的变化。4) 海面粗糙度的Smith模型更加适用于描述我国海上风资源的垂直变化规律。
参考文献
| [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->赵燕华, 蒋翠花, 吴杰. 宿迁70m高度测风塔风速和风能的变化特征分析[C]//S8大气探测与仪器新技术、新方法. 北京:中国气象学会, 2012:349-353. ZHAO Yanhua, JIANG Cuihua, WU Jie. Analysis of characteristics of wind resource based on 70-meter wind measurement mast in Suqian[C]//S8 New Technology and Method of Atmospheric Detection and Instrument. Beijing:Chinese Meteorological Society, 2012:349-353. (in Chinese) |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.-->符平, 秦鹏飞, 张金接. 海上风资源时空特性研究[J].中国水利水电科学研究院学报,2014, 12(2): 155–161.FU Ping, QIN Pengfei, ZHANG Jinjie. Research of spatio-temporal characteristics of offshore wind resources[J].Journal of China Institute of Water Resources and Hydropower Research,2014, 12(2): 155–161.(in Chinese) |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->李常春. 风资源评估方法研究[D]. 呼和浩特:内蒙古工业大学, 2006. LI Changchun. Study on the Method of Wind Resource Assessment[D]. Huhhot:Inner Mongolia University of Technology, 2006. (in Chinese) |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->杜燕军, 冯长青. 风切变指数在风电场风资源评估中的应用[J].电网与清洁能源,2010(5): 62–66.DU Yanjun, FENG Changqing. Application of wind shear index in the assessment of wind resources of wind farm[J].Power System and Clean Energy,2010(5): 62–66.(in Chinese) |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.-->王有禄, 李淑华, 宋飞. 风电场测风数据的验证和处理方法[J].电力勘测设计,2009, 16(1): 60–66.WANG Youlu, LI Shuhua, SONG Fei. Treatment method and inspection of anemometric dates in wind power plant[J].Electric Power Survey and Design,2009, 16(1): 60–66.(in Chinese) |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Businger J A, Wyngaard J C, Izumi Y, et al. Flux-profile relationships in the atmospheric surface layer[J].Journal of the Atmospheric Sciences,1971, 28(2): 181–189. |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Gryning S E, Batchvarova E, Brümmer B, et al. On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer[J].Boundary-Layer Meteorology,2007, 124(2): 251–268. |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Venkatram A. Estimating the Monin-Obukhov length in the stable boundary layer for dispersion calculations[J].Boundary-Layer Meteorology,1980, 19(4): 481–485. |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lange B, Larsen S, H?jstrup J, et al. Importance of thermal effects and sea surface roughness for offshore wind resource assessment[J].Journal of Wind Engineering & Industrial Aerodynamics,2004, 92(11): 959–988. |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Högström U L F. Non-dimensional wind and temperature profiles in the atmospheric surface layer:A re-evaluation[J].Boundary-Layer Meteorology,1988, 42(1-2): 55–78. |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Charnock H. Wind stress on a water surface[J].Quarterly Journal of the Royal Meteorological Society,1955, 81(350): 639–640. |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Smith S D. Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature[J].Journal of Geophysical Research:Oceans,1988, 93(C12): 15467–15472. |
