隧道施工中，常因为施工条件、施工工艺以及地质环境因素(岩溶空腔)的影响，在硐室壁上形成空洞缺陷，导致隧道开挖轮廓面呈现不规则、不光滑的特点[1-2]。与理想设计硐室相比，存在缺陷的不光滑硐室会导致围岩应力发生变异，从而导致岩爆、崩塌等地质灾害或支护与衬砌受力条件的恶化，并对隧道稳定性造成负面影响[3-5]。因此，针对此类硐室围岩应力的求解及分析具有非常重要的意义。
目前硐室围岩应力的计算方法可分为数值法与解析法两类。与数值法相比，解析法具有运算速度快、在应力高度集中的拐角区域精度较高和易进行参数分析的特点[6-7]。目前基于不同工程背景的围岩应力的解析方法不断涌现。针对一般圆形硐室的围岩应力问题，可利用弹性力学理论及叠加原理解决[8]。焦春茂等[9]采用积分算子解法避开了Laplace逆变换的困难，计算分析了含有弹性支护的圆形硐室围岩应力分布情况。Wang等[10]利用Laplace变化和Schwarz交替法得到了两个圆形硐室相互影响下的围岩应力解析解。以上研究均为不同条件下圆形硐室的应力解析解，而针对非圆形硐室，复变函数理论是非常有效的计算工具。Muskhelishvili首次将复变函数理论引入到求解硐室围岩应力的应用中[11]。随后，很多学者基于复变函数理论对围岩应力求解开展了大量研究。Zhao等[12]利用复变函数理论求解了矩形硐室的围岩应力解析解。饶军应等[13]利用保角映射将椭圆外域映射为单位圆外域，求得了含有内压的椭圆硐室围岩应力解析解。Evangelia Pelli等[14]通过贝塞尔型无穷级数展开式估算波势，并结合复变函数理论，得到SV波荷载下硐室应力分布情况。以上研究在分析围岩应力时，均以理想的设计硐形为计算模型，并未考虑在硐室开挖过程中开挖轮廓存在缺陷的情况。而缺陷的存在势必对围岩应力的分布情况产生较大的影响[3]。
本文基于弹性力学和复变函数理论，利用逐次渐进法建立含缺陷硐室外域与单位圆外域之间的映射关系，结合柯西积分公式，最终得到含空洞缺陷圆形硐室围岩应力的解析解。以此分析硐室围岩的应力分布情况，并进行参数分析，讨论缺陷位置、凸出程度以及场地应力对硐室周围围岩应力分布的影响，为工程开挖支护提供一定的理论基础。
1 问题描述与弹性力学基本解
1.1 基本假设
1) 埋深远大于硐室直径时，可以忽略重力梯度的影响，将重力作用简化为无穷远处最大主应力和最小主应力。
2) 围岩介质为均质、各向同性的弹性体。
3) 将空洞缺陷边界简化为弧线。
1.2 计算模型简介
取硐室某一横截面为研究对象，主应力σ1与水平方向的夹角为α。真实平面(Z平面)上的孔洞边界线为L，内域为S+,外域为S-，L的走向以物体位于边界线左侧为正，即顺时针方向为正。含空洞缺陷的圆形硐室围岩应力计算模型如图1所示。
图1 含空洞缺陷的圆形硐室围岩应力计算模型
Fig.1 Calculation model for surrounding rock of a circular chamber with a void defect
1.3 映射函数
利用映射函数[15-18]，将Z平面上的孔洞外域与复平面单位圆外域建立映射关系，函数形式下：
式中，ζ代表单位圆外域上的点，在复平面上采用极坐标表示时，ζ=ρeiθ 。
采取逐次渐进法[19]计算映射函数，步骤如下：
1) 分别令真实孔洞和单位圆与水平轴正向的交点为z0、ζ0，并按逆时针在孔边取m个等分点，且单位圆上的点ζk(1,θk )与孔洞上的点zk(xk,yk)对应，如图2所示。
2) 由单位圆和孔洞上一一对应的两组点坐标，可以得到一个映射函数，具体公式如下：
式中，
为Z平面上将坐标轴逆时针转动θn后，点zk在新坐标系中的X轴坐标值。
3) 将单位圆上的m个点代入由第2)步得到的映射函数，可得到Z平面上的一组点
并对其进行幅角修正，若点
不在孔洞边界上，则用边界上与该点幅角相同的点
替代，由此可得到单位圆和孔边上一一对应的两组点。
重复步骤2)、步骤3)直至精度达到要求。
图2 真实平面与复平面上各点对应关系
Fig.2 Correspondence between real plane and complex plane
1.4 硐室围岩应力的复式表达式
本文重点研究含缺陷硐室孔边以及围岩内部的应力情况，在极坐标系下硐室应力的复式表达式如下式[20]：
式中：σθ为切向正应力；σρ为径向正应力；ϕ(ζ)和ψ(ζ)为复函数。
2 复变函数解
2.1 映射函数求解
当洞形关于x轴对称时，其映射函数的各项系数均为实数[21]，所以使真实洞形关于x轴对称。在单位圆和孔洞边界上分别取9000个点，映射函数项数取42项，迭代20次时，由映射函数得到的近似孔洞与真实孔洞形状重合，如图3所示。
图3 近似孔洞与真实孔洞的对比情况
Fig.3 Comparison between approximate hole and real hole
2.2 应力边界条件
在边界线L上，边界应力值fx和fy为边界线正向的右部对左部的表面作用力。设某边界点A为参考点，B为边界上任意点，则力边界条件在Z平面的复式表达式为[22]:
设ζ平面上的边界线为Γ，σ为Γ上的任意一点，利用映射函数式(1)可将Z平面上的边界条件转到ζ平面上，可得到：
2.3 柯西积分公式
若f(z)在内域上解析，在内域和边界上连续，则柯西积分公式为[23]：
若f(z)在外域上解析，在外域和边界上连续，则柯西积分公式为：
2.4 复函数求解
柯西积分解法是求解复函数的常见方法，其一般性求解过程如下[24]。
复函数的基本形式为：
其中：
因系数A1、B1由孔边外荷载确定，而在本文中孔边不存在外荷载，所以均为0。
将式(8)代入应力边界条件，并对每项都乘以
并沿Γ正向做回路积分，利用柯西积分公式整理可得：
将式(8)代入应力边界条件后，对每项均取共轭，并乘以
沿Γ正向做回路积分，整理可得：
其中：
为计算方便，可将映射函数展开为：
所以，将式(12)代入式(9)右端项，得：
为求解式(13)作如下函数：
可以看出，f0(σ)是f1(ζ1)的边界值。其中，
在内域中只有
一个n级极点，所以可展开为：
式(15)中各项系数可求出，所以结合柯西积分公式，可将式(9)的右端项化简为：
由式(8)可得到：
所以，式(9)左端的积分项可展开为：
为求解式(18)作以下函数：
f2(ζ1)在内域中只有ζ1=0一个(n-2)级极点，所以其洛朗级数展开式为：
所以结合柯西积分公式，式(18)可化简为：
因为ϕ0(ζ)未知，不能直接求出式(21)中的bk，将式(16)和式(21)代入式(9)，根据同幂次项系数相等可得到如下方程组：
由式(15)可知，
将式(23)、式(17)代入式(20)中，通过ζ的幂次项对比，可得到一组αk和bk的关系，与式(22)联立，得到αk和bk。将其带入式(8)即可求出ϕ(ζ)。
同理，式(10)的右端项利用柯西积分公式可化简为：
式(10)左端的积分项，利用柯西积分公式可化简为：
将式(24)和式(25)代入式(10)可求出ψ0(ζ)。
将式(26)代入式(8)可得到ψ(ζ)。将ψ(ζ)和ϕ(ζ)代入式(2)即可求出复平面上ζ点的应力，将ζ带入映射函数中，即可求出该点在真实平面上的位置。
3 解析解的验证
通过特定算例对解析解进行验证。隧道断面如图4所示，硐室半径5 m，其圆心位于原点，缺陷的曲率半径为2 m，且其圆心位于(4,0)点，场地应力情况为σ1=σ2=1 MPa，主应力σ2方向与水平方向的夹角为0°。
图4 解析解验证模型
Fig.4 Verification model of analytical solution
通过与数值解进行对比，对本文所求解析解进行验证。采用有限元分析软件ANSYS建立二维平面应变模型，模型计算范围：水平方向左、右两边各取10倍硐径宽，上、下也各取10倍硐宽。模型上边界和左边界分别施加1 MPa，固定下边界的竖直位移和右边界的水平位移，计算模型如图5所示。
图5 有限元模型
Fig.5 Finite element model
孔边环向应力解析解与数值解的对比情况如图6所示，二者变化趋势一致，解析解与数值解基本吻合。
图6 含缺陷硐室围岩应力解析解与数值解对比
Fig.6 Comparison between analytical solution and numerical solution
4 含空洞缺陷硐室围岩应力分布
仍考虑各向等压(σ1=σ2=1 MPa)条件下硐室围岩中的应力情况，缺陷位置和凸出程度与图4中的模型保持一致。含缺陷硐室孔边的环向应力与不含缺陷硐室孔边环向应力对比情况如图7所示。由图7可以发现，在空洞缺陷部位的环向应力与无缺陷情况相比具有明显差异，而远离缺陷部位应力分布情况与无缺陷情况相同。在圆形硐室与缺陷相交处环向应力较小接近于0，而缺陷中部则发生应力集中现象，环向应力最大值达到了3.3 MPa。
图7 含缺陷与无缺陷硐室在孔边的环向应力对比
Fig.7 Comparison of circumferential stress around holes between chambers with and without a defect
因在硐室表面径向应力为零，图8给出了1.2倍硐径处有无缺陷时径向应力分布情况。从图8中可以看出缺陷中间部位的径向应力具有一定的应力集中现象，其径向应力最大值是无缺陷硐室围岩中径向应力的2.2倍。
图8 含缺陷与无缺陷硐室径向应力对比
Fig.8 Comparison of radial stress between chambers with and without a defect
为分析含缺陷硐室围岩应力由孔边向围岩内部的变化规律，求解围岩中不同倍数硐径上的应力情况。极坐标系下，将孔洞边界上各点的极径放大数倍，即可将孔洞边界线放大相应倍数，如图9所示。
图9 放大后的边界线
Fig.9 Enlarged boundary lines
利用本文得到的应力解析解便可得到孔洞外部不同倍数边界线上的应力分布情况。各曲线上的径向应力以及环向应力分布情况如图10所示。从图10中可以看出向围岩内部延伸时，环向应力与径向应力在缺陷中部的集中现象逐渐减弱。环向应力在1.7倍硐径时基本趋于均匀分布，而径向应力在1.7倍硐径处仍然具有明显的应力集中现象，说明缺陷对硐周径向应力的影响范围明显大于对环向应力的影响。
图10 围岩内部应力分布情况
Fig.10 Stress distribution in surrounding rock
5 参数分析
目前的灾害判据常常根据孔边环向应力对硐室进行评价[25]，本节将分析缺陷凸出程度、缺陷相对场地主应力方向的位置以及场地应力比这三个因素对孔边环向应力的影响。
5.1 缺陷凸出程度对孔边应力的影响
在各向等压(σ1=σ2=1 MPa)情况下讨论缺陷凸出程度对孔边应力的影响。令硐室半径保持5 m不变，缺陷圆弧的曲率半径为2 m，凸出程度由相对半径比值T表征，即由缺陷上的点到硐室圆心的最大距离(R)与硐室半径(R1)的比值。计算模型如图11所示。
图11 讨论缺陷凸出程度对孔边应力影响的计算模型
Fig.11 Calculation model for considering influence of defect protrusion
利用ANSYS软件以及本文得到的应力解析解，分别求解T=1.02、T=1.06、T=1.1以及T=1.2时硐室孔边环向应力。不同凸出程度下的孔边环向应力如图12所示。由图12可知，在缺陷的相对半径比值T仅为1.02时，硐室最大环向应力达到2.48 MPa，远高于远离缺陷部位的2 MPa，说明硐室上的轻微缺陷便能引起较大的应力集中现象。另外，对比不同凸出距离下的环向应力可以发现，随着缺陷凸出程度的增大，缺陷中部的应力集中现象则越来越强烈，而硐室与空洞缺陷相交处的环向应力急剧降低并向零应力的情况逼近。
图12 不同凸出程度下的孔边环向应力分布图
Fig.12 Circumferential stress distribution at different protrusion levels
5.2 缺陷位置对孔边应力的影响
在缺陷凸出程度及曲率不变并且双向受压的情况下，讨论缺陷与主应力的相对位置对孔边应力分布的影响。令σ1=1 MPa、σ2=0.6 MPa，主应力σ2与水平轴夹角为0°，缺陷与主应力的相对位置利用缺陷圆心与原点的连线同X轴正向的夹角来表征，计算模型如图13所示。
图13 讨论缺陷位置对孔边应力影响的计算模型
Fig.13 Calculation model for considering influence of defect position
利用本文得到的应力解析解以及ANSYS软件，分别求解ϕ=0°、ϕ=45°、ϕ=90°以及不存在缺陷时的孔边环向应力。缺陷处于不同位置时的孔边环向应力与无缺陷硐室应力对比情况如图14所示。由图14可以发现，解析解与数值解的计算结果基本一致。当σ1=1 MPa、σ2=0.6 MPa时，无缺陷硐室的孔边环向应力呈“花生形”分布，环向应力最大值出现在硐室两侧，应力值为2.4 MPa，最小值则出现在硐室顶部和底部，应力值为0.8 MPa。而当硐室存在空洞缺陷时，缺陷中部出现应力集中现象，远离缺陷的部位，应力分布与无缺陷硐室情况相同。
通过对比图14(b)～图14(d)可知，ϕ=0°时，缺陷处的环向应力最大值为3.77 MPa，ϕ=45°时缺陷处应力最大值降至2.67 MPa，而当ϕ=90°时，缺陷处应力最大值为1.57 MPa，远低于远离缺陷部位的应力最大值2.4 MPa。说明随着缺陷凸出方向与主应力σ2方向之间夹角(ϕ)的增大，缺陷部位应力集中现象逐渐减弱。当缺陷凸出方向与主应力σ2方向平行时，应力集中现象最强烈，而两者处于垂直状态时，孔边应力分布与无缺陷情况相似，仅在缺陷部位有较弱的应力集中现象。
图14 缺陷处于不同位置时的孔边环向应力与无缺陷硐室情况对比
Fig.14 Comparison of stress distribution between chambers with a defect at different positions and chambers without defects
5.3 场地应力比对孔边应力的影响
在缺陷凸出程度及曲率不变的情况下，讨论缺陷处于不同位置时，场地应力比σ2/σ1对孔边应力的影响。令σ1=1 MPa、σ2=λσ1，λ值分别取1、0.8、0.6、0.4、0.3和0.2，主应力σ2方向与水平轴夹角为0°。计算模型如图15所示。
图15 讨论应力比λ对孔边应力影响的计算模型
Fig.15 Calculation model for considering influence of stress ratio
利用本文得到的应力解析解，分别计算ϕ=0°、ϕ=45°、ϕ=90°时不同应力比下的孔边环向应力。不同应力比下的孔边环向应力如图16所示。由图16可知，随着场地应力比λ的减小，硐室顶部和底部环向应力逐渐由压应力向拉应力转化，而硐室两侧压应力值逐渐增加。
当ϕ=0°时随着λ的减小，缺陷中部的环向应力由3.34 MPa逐渐增大最终提高至4.21 MPa，而当ϕ=90°时随着λ的减小，缺陷中部的环向应力由3.34 MPa逐渐减小，最终降低至0.22 MPa。说明场地应力比对硐室缺陷部位应力的影响与缺陷位置有关：当缺陷位于主应力σ1方向附近时，随着λ的降低缺陷处的压应力集中现象逐步削弱，并逐渐向拉应力转化。而当缺陷位于主应力σ2方向附近时，随着λ的降低缺陷处的应力集中现象进一步增强。
图16 不同应力比下孔边环向应力汇总
Fig.16 Summary of circumferential stress around hole at different stress ratios
6 结论
基于平面弹性复变理论，将含空洞缺陷圆形硐室的围岩应力问题，简化为求解不规则孔洞的双向受压无限平面问题，建立了含缺陷硐室外域和单位圆外域之间的映射关系，并得到含空洞缺陷圆形硐室围岩应力解析解，分析了围岩应力的分布规律以及缺陷位置、凸出程度以及场地应力比对孔边应力分布情况的影响，具体结论如下：
(1) 含空洞缺陷硐室的孔边应力在缺陷中部产生应力集中现象，而在圆形硐室与缺陷相交处应力较小。如σ1=σ2=1 MPa时，缺陷中部的切向应力达到3.3 MPa,而硐室与缺陷相交处的应力趋近于0 MPa。
(2) 较轻微的凸出程度便能引起较大的应力集中现象。在缺陷的相对半径比值T仅为1.02时，硐室最大环向应力达到2.48 MPa，远高于远离缺陷部位的2 MPa。而且随缺陷凸出程度的增加，缺陷中部应力集中现象明显加剧。
(3) 随着缺陷凸出方向与主应力σ2方向夹角的减小，缺陷中部的应力集中现象加剧。当凸出方向与主应力σ2方向垂直时，应力集中程度最弱，且随着应力比的降低，集中现象进一步削弱。而当两者垂直时，应力集中现象最强，随着应力比的降低，应力集中现象进一步加强。
本文仅考虑了单个空洞缺陷圆形硐室的平面弹性应力状态。
[1]王飞, 赵晓勇, 马宏韬.水平岩层隧道超挖原因及合理超挖数量分析[J].铁道工程学报, 2018, 35(4): 75―80.Wang Fei, Zhao Xiaoyong, Ma Hongtao.Analysis of the overbreaking reason and reasonable overbreaking quantity in the horizontal rock stratum tunnel [J].Journal of Railway Engineering Society, 2018, 35(4): 75―80.(in Chinese)
[2]杨柳, 耿占勇, 唐海军.隧道超欠挖原因分析及控制措施[J].价值工程, 2017, 36(32): 170―171.Yangliu, Geng Zhanyong, Tang Haijun.Cause analysis and control measures of tunnel over and under excavation[J].Value Engineering, 2017, 36(32): 170―171.(in Chinese)
[3]管维威.钻爆法隧道施工中的弹塑性力学问题分析[J].科技创业月刊, 2012, 25(7): 176―178.Guan Weiwei.Analysis of elastoplastic mechanics in tunnel construction by drilling and blasting [J].Pioneering with Science and Technology Monthly, 2012,25(7): 176―178.(in Chinese)
[4]孙振宇, 张项立, 房倩, 等.基于超前加固的深埋隧道围岩力学特性研究[J].工程力学, 2018, 35(2): 92―104.Sun Zhenyu, Zhang Xiangli, Fang Qian, et al.Research on the mechanical property of the surrounding rock of deep-buried tunnel based on the advanced reinforcement[J].Engineering Mechanics, 2018, 35(2): 92―104.(in Chinese)
[5]陈仁朋, 王诚杰, 鲁立, 等.开挖对地铁盾构隧道影响及控制措施[J].工程力学, 2017, 34(12): 1―13.Chen Renpeng, Wang Chengjie, Lu Li, et al.Influence of excavation on exist metro shield tunnel and control measures [J].Engineering Mechanics, 2017, 34(12): 1―13.(in Chinese)
[6]Kargar A R, Rahmannejad R, Hajabasi M A.A semi-analytical elastic solution for stress field of lined non-circular tunnels at great depth using complex variable method [J].International Journal of Solids &Structures, 2014, 51(6): 1475―1482.
[7]Wang H N, Zeng G S, Jiang M G, Analytical stress and displacement around non-circular tunnels in semi-infinite ground [J].Applied Mathematical Modelling, 2018,63(11): 303―328.
[8]蔡美峰.岩石力学与工程[M].第2版.北京: 科学出版社, 2013: 283―302.Cai Meifeng.Rock mechanics and engineering [M].Beijing: Science Press, 2013: 283―302.(in Chinese)
[9]焦春茂, 吕爱钟.粘弹性圆形巷道支护结构上的荷载及其围岩应力的解析解[J].岩土力学, 2004, 25(增刊1):103―106.Jiao Chunmao, Lv Aizhong.Analytical solution of loads on supporting structure for circular tunnel and stresses in viscoelastic surrounding rock [J].Rock and Soil Mechanics, 2004, 25 (Suppl 1): 103―106.(in Chinese)
[10]Wang H N, Zeng G S, Utili S, et al.Analytical solutions of stresses and displacements for deeply buried twin tunnels in viscoelastic rock [J].International Journal of Rock Mechanics and Mining Sciences, 2017, 93(3): 13―29.
[11]Muskhelishvili N I.Some basic problems of the mathematical theory of elasticity [M].Leyden:Noordhoff International Publishing, 1953: 28―164.
[12]Zhao G, Yang S.Analytical solutions for rock stress around square tunnels using complex variable theory [J].International Journal of Rock Mechanics & Mining Sciences, 2015, 80(14): 302―307.
[13]饶军应, 傅鹤林, 刘运思, 等.管状填充性椭圆形溶洞的围岩应力弹性解析分析[J].中南大学学报(自然科学版), 2015, 46(7): 2605―2612.Rao Junying, Fu Helin, Liu Yunsi, et al.Stress analysis of rocks surrounding a tubular filled elliptical karst cave with complex function of elastic mechanics [J].Journal of Central South University (Science and Technology),2015, 46(7): 2605―2612.(in Chinese)
[14]Evangelia Pelli, Alexandros I.Sofianos, Analytical calculation of the half space stress field around tunnels under seismic loading of SV waves [J].Tunnelling and Underground Space Technology, 2018, 79(9): 150―174.
[15]皇甫鹏鹏, 伍法权, 郭松峰, 等.基于边界点搜索的洞室外域映射函数求解法[J].岩土力学, 2011, 32(5):1418―1424.Huangfu Pengpeng, Wu Faquan, Guo Songfeng, et al.A new method for calculating mapping function of external area of cavern with arbitrary shape based on searching points on boundary [J].Rock and Soil Mechanics, 2011,32(5): 1418―1424.(in Chinese)
[16]Wang H N, Zeng G S, Jiang M J.Analytical stress and displacement around non-circular tunnels in semi-infinite ground [J].Applied Mathematical Modelling, 2018,63(11): 303―328.
[17]祝江鸿, 杨建辉, 施高萍, 等.单位圆外域到任意开挖断面隧洞外域共形映射的计算方法[J].岩土力学,2014, 35(1): 175―183.Zhu Jianghong, Yang Jianhui, Shi Gaoping, et al.Calculating method for conformal mapping from exterior of unit circle to exterior of cavern with arbitrary excavation cross-section [J].Rock and Soil Mechanics,2014, 35(1): 175―183.(in Chinese)
[18]Lu A, Zhang N, Kuang L.Analytic solutions of stress and displacement for a non-circular tunnel at great depth including support delay [J].International Journal of Rock Mechanics and Mining Sciences, 2014, 70(6): 69―81.
[19]王桂芳.隧道计算:共同变形理论[M].成都: 成都科技大学出版社, 1992.Wang Guifang.Tunnel computing: Theory of joint deformation [M].Chengdu: Chengdu University of Science and Technology press, 1992.(in Chinese)
[20]Xu M F, Wu S C, Gao Y T, et al.Analytical elastic stress solution and plastic zone estimation for a pressure-relief circular tunnel using complex variable methods [J].Tunnelling and Underground Space Technology, 2019,84(2): 381―398.
[21]陈峰宾, 张顶立, 张素磊.采用逐次渐进法求解非圆形地下洞室映射函数[J].中国铁道科学, 2013, 34(2):80―86.Chen Fengbin, Zhang Dingli, Zhang Sulei.Solving the mapping function of noncircular underground cavern with the method of successive approximation [J].China Railway Science, 2013, 34(2): 80―86.(in Chinese)
[22]陆明万, 罗学富.弹性理论基础.下册-第2版[M].北京: 清华大学出版社, 2001: 1―36.Lu Mingwan, Luo Xuefu.Fundamentals of elasticity theory.Volume II - Second Edition [M].Beijing:Tsinghua University Press, 2001: 1―36.(in Chinese)
[23]Мусхелишвили.数学弹性力学的几个基本问题[M].北京: 科学出版社, 1958.Мусхелишвили.A few basic problems in mathematical elasticity [M].Beijing: Science Press, 1958.(in Chinese)
[24]陈子荫.围岩力学分析中的解析方法[M].北京: 煤炭工业出版社, 1994: 43―87.Chen Zi Yin.Analytical method for mechanical analysis of surrounding rock [M].Beijing: Coal Industry Press,1994: 43―87.(in Chinese)
[25]张镜剑, 傅冰骏.岩爆及其判据和防治[J].岩石力学与工程学报, 2008, 27(10): 2034―2042.Zhang Jingjian, Fu Bingjun.Rockburst and its criteria and control [J].Chinese Journal of Rock Mechanics and Engineering, 2008, 27(10): 2034―2042.(in Chinese)