地球物理学报  2010, Vol. 53 Issue (2): 457-465   PDF    
煤矿微震监测台网优化布设研究
巩思园1,2 , 窦林名1,2 , 曹安业1,2 , 贺虎1,2 , 杜涛涛1,2 , 江衡1,2     
1. 中国矿业大学煤炭资源与安全开采国家重点实验室, 江苏徐州 221008;
2. 中国矿业大学矿业工程学院, 江苏徐州 221006
摘要: 煤矿冲击矿压现象日益严峻,作为矿山动力灾害的主要监测手段,微震监测系统已在许多矿井广泛使用,为保证矿震定位和能量计算的准确,提高预测预报煤矿冲击矿压的可行性,应建立一套台网布设优化及评价系统.应用微震定位和D值优化设计理论,结合煤矿实际条件研究了影响矿震定位精度的主要因素和不利条件,并提出了采用综合指数法确定煤矿高微震活动区域和区域内矿震发生的概率,制定了台站候选点和监测区域确定的一般原则.通过理论分析震中和震源标准差反映台网定位能力的不足,建立基于数值仿真实验方法的震中与震源误差期望值模型,最终形成台网布设优化及评价系统.实验和现场应用结果表明,该系统能够快速确定台网最优布设方案,准确评价台网定位能力,满足煤矿微震监测的需要.
关键词: 冲击矿压      微震监测系统      台网布设      D值优化设计      定位      综合指数法      数值仿真实验     
Study on optimal configuration of seismological observation network for coal mine
GONG Si-Yuan1,2, DOU Lin-Ming1,2, CAO An-Ye1,2, HE Hu1,2, DU Tao-Tao1,2, JIANG Heng1,2     
1. State Key Laboratory of Coal Resource and Mie Safety, China University of Mining & Technology, Xuzhou, Jiangsu 221008, China;
2. School of Mines, China University of Mining & Technology, Xuzhou, Jiangsu 221006, China
Abstract: The rockburst phenomenon caused mostly by mining-induced tremors becomes progressively more severe as the mining depth and extent increase. Microseismic monitoring technique, as the major investigation tool for the monitoring of rockburst and other mining dynamic hazards, has been widely used in many coal mines around the world. To ensure the accuracy of seismic source locating and energy calculation, and further improve the feasibility of rockburst forecasting, an optimal design and assessment system of network configuration needs to be built first. In this paper, according to practical conditions and requirements of coal mines, especially the accuracy of picking seismic wave arrivals at stations and building velocity model, the main influencing factors and bad conditions which influence the accuracy of source location were analyzed using the theory of D-optimal design and localization. A view of using the comprehensive index method to determine the high microseismic activity zones and the occurrence probability of mining-induced tremors were also proposed. Then, the general principle of choosing candidate points of geophone stations and monitoring area were established. In addition, since the location capacity of network can't be reflected accurately using the source standard error calculated by the traditional theoretical analysis, the model of expected epicenter and hypocenter error was created by numerical emulation experiment, and then the optimal design and assessment system of network was built up finally. The results from experiment and field practice show that, the optimal design and assessment system can be used to solve the optimal plan of network quickly, evaluate the location ability of network exactly, and can satisfy the requirements of microseismic monitoring in coal mines quite well..
Key words: Rockburst      Microseismic observation system      Microseismic network configuration      D-optimal design      Localization      Comprehensive index method      Numerical emulation experiment     
1 引言

近年来,国内许多煤矿都已进入深部开采阶段,以煤矿矿震为代表的深部开采诱发的灾害更具突发性,表现出明显的非线性动力特征,严重时导致冲击矿压的发生,造成煤岩体破坏、支架和设备损坏和人员伤亡[1, 2].为监测煤矿矿震,中国矿业大学与波兰矿山研究总院合作并引进了先进的“SOS”(Seismological Observation System)微震监测系统,并在现场取得了成功的应用.这很大程度上建立在震源定位和能量计算的可靠性和准确性上,而煤矿矿震定位的准确度依赖于以下几种因素[3, 4]:微震台网布设、台站P波到时读入的准确性、背景噪音的特点和仪器的采样频率、求解震源算法、速度模型和区域异常(例如采空区)所导致的传播路径的变化.其中,速度模型和区域异常等因素可通过联合震中测定技术来消除[4~6],而P波到时读入准确性和震源到台站几何特征等随机因素则无法根本消除,只能通过优化台网布设和降低随机因素大小等手段降低求解震源的非线性方程组的条件数,提高台网的容差能力,从而提高求解系统的稳健性.

台网布设优化问题最初来源于地震台网优化布设.Sato和Skoko[7]提出了蒙特卡洛算法进行地震台网监测能力的数值计算研究,并绘制了监测区域震源参数的误差等值线,但是对煤矿领域却没有涉及,参考价值较小.随后,Kijko[8, 9]和Mendecki[3, 10]分别基于D值和C值最优设计理论提出了设计微震台网的方法.D值理论认为震源参数协方差矩阵行列式大小正比于误差椭球体体积,行列式越小,椭球体体积越小,震源参数分布越集中,参数估计就越准确.该方法成功的应用到了多个领域[11, 12].C值理论则从影响非线性系统稳健性的角度分析了台网布设与所形成方程组条件数的关系.以上两种方法虽能够用于评价台网的优劣,但是仍不能准确反应台网的定位能力,也没有提出确定煤矿冲击矿压危险监测区域和台站位置候选点的一般原则.

基于上述存在的问题,本文则根据煤矿矿山的实际条件,研究微震台网优化布设中选取台站候选点和确定冲击矿压危险监测区域的一般原则,利用D值优化准则形成微震监测系统最终方案,并基于统计学提出检验监测方案的数值分析仿真方法,并在现场和人工数据中得到验证,最终形成适于煤矿的微震台网布设优化系统.

2 微震定位理论

在煤矿井田范围尺度下安装微震探头,要保证震源定位具有较高的准确性,通常选择比较容易辨认的P波进行定位,与其他波相比,P波初次进入时间的确定误差较小,定位精度较高.由于初期成本的原因,实际应用中,多假设煤岩体为均质、各向同性介质,即P波在各个传播方向上保持速度不变,为一定值.图 1从震源传播到台站的最短时间可由(1)式描述:

图 1 矿震定位示意图 Fig. 1 Schematic representation of the minet remor location process

(1)

式中,H=(x0y0z0)和Xi=(xiyizi)分别为震源和第i个台站的坐标,V为P波波速,t0为矿震的发震时刻,ti为读入的P波初至到达时刻,i=1,2,… nn是矿井中安装的台站数目.

对于均匀和各向同性速度模型,自震源H到第i个台站的走时为[4]

(2)

式(2)中P波波速VP为一常数.

式(1)有θ=(x0y0z0t0)4个未知数,要解这个方程至少需要4个观测站的数据.目前在各个矿区投入使用的SOS微震系统都采用16个台站的布置形式,所以最多可以列出16个类似于(1)式的方程.为进行震源定位,目标函数可写成如下形式:

(3)

式中,ri为残差,即观测值ti与P波计算到时值TiHVXi)+t0之差,本文中q值定为2.通过求解(3)式的最小值,所求的参数值为参数θ的最小二乘估计.为了估计,通常先提供尝试矢量θ (0),然后以校正矢量δθ n来更新尝试矢量θ n,并减少目标Θ的值.对走时TiHVXi)应用一阶泰勒式线性化后,在每次迭代过程中:

(4)

式中,δrn为在空间内点θ n上的时间残差矢量,A为在θ n上计算的(1)式对参数θ的(n×4)偏微分矩阵.

(5)

3 台网设计的D值最优设计理论

受随机误差ξi影响,(1)式有如下普遍形式:

(6)

通常在台网安装前,对台网内所有台站总是假设随机误差满足相同的正态分布ξ~N(0,σ2I),I为单位矩阵,σ为随机误差的方差.

由以上假设,Gallant[13]确定最小二乘估计参数近似满足正态分布:

(7)

CθX)=(ATA-1σ2为求解参数θ的协方差矩阵,根据正态分布的特点,满足自由度为n-4的χ 2分布为:

(8)

(8)式的几何意义非常重要,描述了在某一置信区间下估计参数的分布特征.对不同的台网布设方案X,根据最小二乘法可获得不同的参数估计,参数估计的好坏则由此式表示的椭球体积大小来评价.体积越小,估计参数的分布就越集中,定位就越准确,布设方案X就越有利.由D值优化准则知,椭球体体积与成正比.协方差矩阵CθX)的行列式越小,椭球体体积越小.满足det[CθX)]最小的台网设计方案X*称为D值最优台网布设方案[4].由估计参数的方差表达式可以看出,并不需要计算协方差矩阵,而可由偏微分矩阵A表示.求(ATA-1的行列式最小值同样满足D值优化准则.在考虑随机误差中P波波速的影响后,协方差可写成(ATWA-1的形式,W为对角矩阵,对角元素分别为:

(9)

式中,σVPσt分别为P波波速和P波到时读入的方差,并且对于所有台站取值相同.

以上D值优化准则仅适用于矿震集中在相对较小区域时,煤矿中实际情况更加复杂.由于许多矿井开采和掘进工作面不止一个,矿震活动危险区域较多,某点上的最优布设方案X*由多个区域组成的整体区域上的最优方案ΩX*所替代.假设在某点Hj上发生矿震的概率为pHj),同样可表述为该点的重要性,则mindet[CθX)]可以由整个区域ΩH中目标函数所替代:

(10)

离散形式为:

(11)

式中,ne为冲击危险区域内需要计算的震源点数量.在计算偏微分矩阵A时,需要代入震源位置Hj和台网布设方案X.

因为监测区域中的大部分震动都不能激发所有台站,况且小的矿震是分析大震动的基础,是重要的信息源,所以偏微分矩阵A的行数必须是变化的.不同能量级的震动有不同的可探测距离,只有处于监测范围的探头才能被激发并包含于矩阵A.利用能量E和可探测距离r的经验公式E=λrq[3]q接近2,λ为某一常值,可以确定某一能量下,某一震源点的接收探头数目,且规定在可探测距离内至少要有5个通道能够接收到波形.可以看出,q值越大,可探测距离就越短,那么在同样能量下,能够触发的探头个数就越少.此时,最优台网布置就越紧密,所能覆盖的区域也越有限.

4 监测区域和台站候选点确定原则

根据D值优化理论,微震系统安装前,为确定台网布设方案,应首先根据影响冲击矿压危险状态的地质因素和开采技术因素确定矿区内须重点监测的高微震活动区域.在分析已发生的各种冲击矿压灾害的基础上,利用综合指数法[1, 2],确定各种因素的影响权重,然后将其综合起来评价各区域内的冲击危险程度,最后由冲击矿压危险状态等级评定综合指数确定区域内发生矿震的概率,如下表 1.

表 1 高微震活动区域内矿震发生概率 Table 1 The probability of mine tremors in high seismic activity zones

确定危险区域后,由于煤矿中受巷道布置、开采、施工和现场条件等因素限制,并不是所有的地点都可以安装微震探头,所以初期必须根据一定的原则,选入一些可行的监测点作为台站位置的候选点,再通过D值优化准则进行优化组合选择,最终确定台网的布设方案.为尽可能避免随机因素中P波波速和P波到时读入误差的影响,减少震源定位的误差,候选点的选择还要考虑所处的环境因素和开采活动的影响.由此确定选择候选点的一般原则为:

(1)危险区域周边应尽量在空间上被候选点均匀包围,候选点数不能少于5个,并避免近似形成一条直线或一个平面,并具有足够和适当的空间密度;

(2)一部分候选点应尽可能接近待测区域,避免较大断层及破碎带的影响.但是受巷道布置的客观条件影响,常见情况如下图 2,接近监测区域的探头只能安装在A和B两条直线巷道中,与原则(1)相悖.为尽量挺高定位精度,一方面增加在A和B中备选探头的数目,但距超前支护段的距离不应小于50m;另一方面结合客观条件考虑在监测区域其他方位的地面上选择合适的候选点.

图 2 候选点选择的不利条件 Fig. 2 The disadvantage for selecting candidate points

(3)候选点应远离大型电器和机械设备的干扰,例如皮带机头,转载机等,尽量利用现有巷道内的躲避硐室,远离行人和矿车影响.为减少波的衰减,探头尽量安装在底板为岩石的巷道内.

(4)既要照顾当前开采区域,又要考虑未来一定时期内的开采活动.

5 台网D值优化方案的数值仿真

A.Kijko定义震中位置标准差为平面圆的半径,该圆的面积等于在(x0y0)处标准误差椭圆的面积[9].由此,确定:

(12)

式中,CθXij为协方差矩阵CθX)的元素.由于椭圆的两个轴对应协方差矩阵的特征值(λx0,λy0),所以(12)式标准误差又可写成(13)式的形式:

(13)

两式的区别在于,(12)式由协方差矩阵CθX)对应X轴和Y轴的子矩阵计算而得,即(13)式中的(λx0λy0)不是来自于协方差矩阵CθX),导致由(11)式计算的震中位置标准差可能会出现大于震源误差的情况.但是Kijko[8, 9]并没有定义相应的震源位置的标准差,与震中标准差定义类似,通过建立球体和椭球体等体积的关系式,利用X轴、Y轴、Z轴方向的协方差矩阵CθX)的特征值(λx0λy0λz0)同样可以计算出σxyz

(14)

时,震中标准差大于震源标准差,与震源定位误差应大于震中误差不符.多出现在台网立体分布较好,且平面分布较差时,但煤矿实际条件下,出现的可能性较小,反而多是出现由于台网空间分布不均匀,密度不高,λz0偏高的情况.而且由于标准误差椭圆或椭球置信区间为39.4%,只能包含近40%的点,所以利用特征值计算的σxyσxyz并不能恰当的反应台网的定位能力.当特征值差异比较大时,λz0将被低估,计算的震源标准差偏低.

为体现台网的定位能力,一种非常简便的方法就是对监测区域进行大量的仿真实验.实验过程主要考虑随机因素P波波速和P波到时读入误差的影响,假设它们在所有台站分别服从相同的正态分布,即.受随机误差污染后,监测区域内Hj点到台站Xi的P波传播时间为

(15)

式中,Dist(HjXi)为Hj到台站Xi的直线距离,〈Vp〉和〈ξ〉为随机产生的样本值.根据微震定位理论,当n≥4时,即可利用污染后的tij计算新的震源位置H'jH'jHj的震中距离和震源距离即可作为污染后的定位误差.在多次重复实验(大样本)下,定位误差的期望值(16)即可作为对台网在Hj上定位能力的评价:

(16)

式中,NmHj点上重复实验次数,一般大于1000.为保证与SOS系统算法的一致性,求解震源位置算法选为鲍威尔[14]算法.

由于仿真实验也可得到定位误差的期望值,所以该方法同样能对微震台网进行定位能力评价,并用于确定最佳台网布局.但是由于重复实验次数较多,而鲍威尔法用于定位计算又需要耗费一定的时间,将导致优化台网布设的计算成本过高.

6 实验和现场应用 6.1 实验台网数值分析仿真与理论结果对比

Steinberg和Rabinowitz[12]推荐了一种有效监测核试验的地震台网监测方案.为验证数值分析仿真的有效性,对选取的实验台网方案(见表 2图 3a)进行震源和震中定位误差分析,并与理论分析的结果进行对比.

图 3 实验台网定位误差理论计算结果(a)台网布设方案A;(b)震中标准差等值线图(公式12);(c)震中标准差等值线图(公式13);(d)震源标准差等值线图(公式14) Fig. 3 The expected location error of experimental network in theory (a) The station arrangement planA; (b) Contour plot of standard epicenter location error calculated by equation 12;(c) Contour plot of standard epicenter location error calculated by equation 13;(d) Contour plot of standard hypocenter location error calculated by equation 14
表 2 用于监测核试验的有效地震台网布设方案 Table 2 The effective seismic network used for the inspection of the Nuclear Test

仿真试验和理论分析计算过程中,取P波期望速度为4000m/s,波速方差为50 m/s,P波首次到时读入时间方差为0.005s.仿真试验中重复试验次数Nm为2000,且假设震动能被所有台站记录.分析区域内X方向区间为[-2000 2000],间距50m,Y方向区间为[-2000 2000],间距50m,平面标高-1000m.

利用D值优化理论,结合公式(12)、(13)和(14)分别计算得图 3b图 3c图 3d.理论上,震源定位误差应大于震中定位误差.但是对比图 3b图 3d,在一定范围内,从图中中心点向外,两者差距逐渐缩小.超出此范围后,在边角区域,震中误差大于震源误差,而改为(13)式计算震中标准差后,误差明显降低,消除由子矩阵带来的影响.

根据数值试验方案,计算的震中误差与震源误差等值线图如图 4.在自由度为9998时,图 3c图 4a图 3d图 4b线性相关系数分别为0.689和0.937,具有较好的对应关系.消除了由于置信区间较小和特征值差异较大时,震源和震中定位误差被低估的现象.

图 4 实验台网定位误差数值仿真结果(a)震中定位误差等值线图(公式16);(b)震源定位误差等值线图(公式16). Fig. 4 The location error of experimental network by numerical emulation (a) Contour plot of epicenter location error calculated by equation 16 (b) Contour plot of hypocenter location error calculated by equation 16
6.2 现场应用

某煤矿为深部开采矿井,曾在多个区域发生较大矿震,并导致冲击矿压发生.为实现对冲击矿压的连续监测和预警,计划安装SOS微震监测系统,台网由16个探头组成.初期根据高微震活动区域和台站候选点选择原则,在图 5中选择了3个冲击危险区域.由综合指数法计算的区域Ⅰ冲击危险指数为0.88,区域Ⅱ冲击危险指数为0.57,区域Ⅲ冲击危险指数为0.34,对照表 1,三个区域内发生矿震的概率分别为0.35、0.65和0.85.随后在三个区域的周围确定了312个候选点,共形成2.61×1026种布设组合方案,理论上该组合问题可解,但由于问题规模太大,时间耗费巨大,不适于采用分支定界或穷举法找到最优解,故本论文选取遗传算法[15, 16]进行求解.由现场实测确定P波期望速度为4100m/s,q为1.85,常值λ为2.8×10-3.波速方差为50 m/s,P波首次到时读入时间方差为0.005s.参照波兰冲击矿压防治规范汇编,为保证能够记录和定位能量至少为1000J的震源,确定用于求最优布设方案的能量为2000J,则可探测距离为1460m.由以上参数最终确定的最优台网布设方案见图 5.

图 5 某矿冲击危险区域划分及求解的台网最优布设方案红色虚线包围的区域为划分的冲击危险区域 Fig. 5 The divided rockburst hazard zones and solved optimum distribution of seismic network in a certain coal mine The regions surrounded by red dash lines are dangerous zones

对区域Ⅲ,采用理论和数值仿真试验得到的震源和震中误差等值线见图 6图 7.在自由度为9998时,图 6a图 7a图 6b图 7b线性相关系数分别为0.879和0.82,两种方法具有良好的耦合关系.应用D值优化理论能够快速确定微震台网的最优布设方案,而数值仿真能够更准确的体现台网的实际定位能力.

图 6 区域Ⅲ内-1000水平定位误差理论计算结果(a)震中标准差等值线图;(b)震源标准差等值线图. Fig. 6 The theoretic location error of-1000 level in region Ⅲ (a) Contour plot of standard epicenter location error; (b) Contour plot of standard hypocenter location error.
图 7 区域Ⅲ内-1000水平定位误差数值仿真结果(a)震中误差等值线图;(b)震源误差等值线图. Fig. 7 The location error of-1000 level in region Ⅲ by numerical emulation (a) Contour plot of epicenter location error; (b) Contour plot of hypocenter location error.

SOS微震监测台网布设完成后,在区域Ⅲ巷道掘进过程中,系统记录了大量清晰的爆破震动波形数据,例如在图 8的5个地点中记录的图 9中的震动波形数据.

图 8 区域Ⅲ中5次爆破所在的地点和定位结果 Fig. 8 The five blasting events and corresponding location results in region Ⅲ
图 9 区域Ⅲ中SOS微震系统记录的5次爆破震动波形数据 Fig. 9 The five explosive seismograms recorded by seismological observation system in region Ⅲ

利用sos系统对5次爆破震动数据进行标记和定位计算,定位结果与实际位置对比,得到表 3中的震中误差和震源误差.可以看出,利用台网优化系统确定的布设方案满足现场震中误差与震源误差分别小于50 m和70 m的要求,能够帮助分析人员准确确定放炮的地点,并且利用数值仿真计算的误差分布与实际结果更接近.

表 3 5次爆破定位震中和震源误差 Table 3 The epicenter and hypocenter locationerror for five blasting events
7 结论

(1)D值最优设计理论为微震台网设计提供了可行的方法,在存在大量候选点时,结合D值优化准则,采用遗传算法可以快速、低成本的确定各种因素条件下的台网最优布设方案.

(2)推广震中标准差至震源标准差形式,但由于标准误差椭圆或椭球置信区间只能包含近40%的点,并不能恰当的反应台网的定位能力,特别当特征值差异比较大时,计算的震源标准差偏低.

(3)采用数值仿真实验,在大样本下求震中误差和震源误差的期望值,与D值最优设计理论有很好的耦合关系,相比D值优化准则能够更准确的评价台网的定位能力.

(4)采用综合指数法可确定煤矿矿井中高微震活动区域及区域内发生矿震的概率,保证了台网在重点区域内的监测精度.

(5)根据台站候选点确定原则可快速形成求解最优台网布设的可行点集合,确保最优解的有效性.

(6)实验和现场应用结果表明煤矿微震台网布设优化及评价系统是有效可行的.

参考文献
[1] 窦林名, 何学秋. 冲击矿压防治理论与技术. 徐州: 中国矿业大学出版社, 2001 . Dou L M, He X Q. Theory and Technology of Rock Burst Prevention (in Chinese). Xuzhou: China University of Mining and Technology Press, 2001 .
[2] 齐庆新, 窦林名. 冲击地压理论与技术. 徐州: 中国矿业大学出版社, 2008 . Qi Q X, Dou L M. Theory and Technology of Rock Burst (in Chinese). Xuzhou: China University of Mining and Technology Press, 2008 .
[3] Mendecki a J. Seismic Monitoring in mines. London: Chapman and Hall Press, 1997 .
[4] Gibowicz S J, Kijko A.矿山地震学引论.修济刚译.北京:地震出版社, 1996 Gibowicz S J, Kijko A.An Introduction to Mining Seismology.Translated by Xiu Jigang.Beijing:Scismological Press, 1996
[5] Douglas A. Joint epicentre determination. Nature , 1976(215): 47-48.
[6] Lilwall R C, Douglas A. Estimation of P-wave travel times using the joint epicenter method. Geophys.J.Roy.Astr.Soc , 1970, 19: 165-181. DOI:10.1111/gji.1970.19.issue-2
[7] Sato Y, Skoko D. Optimum distribution of seismic observation points. .IL Bull.of Earthquake Res.Inst , 1965, 43: 451-457.
[8] Kijko A. An algorithm for optimum distribution of a regional seismic network-I. Pageoph , 1977(115): 999-1009.
[9] Kijko A. An algorithm for the optimum distribution of a regional seismic network-Ⅱ.An analysis of the accuracy of location of local earthquakes depending on the number of seismic stations. Pageoph , 1977(115): 1011-1021.
[10] Mendecki A J, Van Aswegen G.A method for the optimal design of mine seismic networks in respect to location errors and its application.Research report, AAC Research and Development Services, Rock Mechanics Department, 56
[11] 唐礼忠, 杨承祥, 潘长良. 大规模深井开采微震监测系统站网布置优化. 岩石力学与工程学报 , 2006, 25(10): 2036–2042. Tang L Z, Yang B X, Pan C L. Optimization of microseiismic monitoring network for large-scale deep well mining. Chinese Journal of Rock Mechanics and Engineering (in Chinese) , 2006, 25(10): 2036-2042.
[12] Steinberg D M, Rabinowitz N. Optimal seismic monitoring for event location with application to On Site Inspection of the Comprehensive Nuclear Test Ban Treaty. Metrika , 2003(58): 31-57.
[13] A Ronald Gallant. Nonlinear Statistical Models. New York: Wiley, 1987 .
[14] 朱永生. 实验物理中的概率和统计. 北京: 科学出版社, 2006 . Zhu Y S. Probability and Statistical in Physics Experiments (in Chinese). Beijing: Science Press, 2006 .
[15] Angus Errington. Sensor Placement for Microseismic Event Location. Saskatoon: University of Saskatchewan, 2006 .
[16] 巩思园, 窦林名, 李志华, 等. 混合遗传算法在矿震定位中的应用.2008全国冲击地压研讨会论文集.. 徐州: 中国矿业大学出版社, 2008 : 137 -141. Gong S Y, Dou L M, Li Z H, et al. Mixed genetic algorithm application on the location of mine tremor.Proc.of the coal mine rockburst (in Chinese). Xuzhou: China University of Mining and Technology Press, 2008 : 137 -141.