考虑应力敏感影响的改进Blasingame产量递减分析方法

陈民锋，王兆琪，孙贺东，王言，陈璐，张琪琛

1 中国石油大学(北京)石油工程学院，北京 102249

2 中国石油勘探开发研究院廊坊分院，廊坊 065000

3 新疆油田勘探开发研究院，克拉玛依 834000

考虑应力敏感影响的改进Blasingame产量递减分析方法

陈民锋1*，王兆琪1，孙贺东2，王言3，陈璐1，张琪琛1

1 中国石油大学(北京)石油工程学院，北京 102249

2 中国石油勘探开发研究院廊坊分院，廊坊 065000

3 新疆油田勘探开发研究院，克拉玛依 834000

Blasingame产量递减分析方法，可用于分析变井底压力和变产量条件下的产量递减问题，但Blasingame方法并没有考虑储层应力敏感的影响。本文基于低渗透应力敏感储层渗流特征，建立考虑储层应力敏感影响的定产生产模型，并采用不等距差分的数值方法进行求解，并在计算过程中将定产生产条件转化为内边界条件、显示处理应力敏感对储层渗透率的影响；基于Blasingame方法的无因次变化量，制作应力敏感储层的产量递减曲线图版，并分析了不同条件下的产量变化规律。研究表明，考虑应力敏感的改进Blasingame产量递减曲线，在边界控制流阶段曲线的形态变化互相平行，而随着无因次距离的增加，曲线向下移动；随着应力敏感系数的增加，改进Blasingame产量递减曲线向下移动。实际应用表明，应用改进的Blasingame分析方法，解释结果更符合低渗透、应力敏感储层的开发实际状况。

应力敏感；不等距差分；改进Blasingame方法；产量递减曲线；解释图版

0 引言

对于应力敏感储层，在衰竭式开发中，压力变化对渗透率的影响是不可忽略的。目前针对应力敏感的实验研究进行了很多[1-3]，得到了明确的渗透率随压力变化的关系，然而针对应力敏感储层进行产量递减分析的解释方法目前还相对较少[4-6]。

目前大多数学者用解析方法处理包含应力敏感的数学模型时，主要采用Pedrosa变换线性化处理，再结合摄动技术在拉氏空间解析求解。2001年段永刚[7]和2009年孙来喜[8-9]，分别利用该方法研究了应力敏感油藏在无限大地层中、定压生产的产量递减规律，得到对于定压生产的油气井在边界控制流阶段产量呈指数递减。而在实际生产过程中，该阶段油气井产量递减变慢，并不符合指数递减规律，而是服从双曲递减规律。因此定压生产条件下的产量递减规律，并不适用于油气井在整个开发阶段的产量递减分析。在油气藏衰竭开发阶段，实际油气井常根据自身特点，在不同开发生产阶段往往采取不同的、相对定产的生产制度，因此研究定产条件下、考虑应力敏感影响的产量递减规律，是油气藏精细化开发的重要基础。

在实际生产过程中，油气井产量是不断递减变化的，Blasingame在1993年提出了物质平衡时间来处理实际时间与实际产量的变化关系，可将变产量的问题转化为定产量生产的问题进行分析。因此可以根据定产数学模型的理论解，制作不同条件下生产曲线的理论图版，进行产量递减规律分析。

在现代产量递减分析中[10-12]，以Blasingame方法为代表的双对数拟合分析法应用较广[13-14]。Blasingame方法的复合图版包含3条曲线，即规整化产量曲线、规整化产量积分曲线、规整化产量积分导数曲线，这种方法不仅可以解释变井底流压下生产情况，也可以处理衰竭开发过程中流体性质的变化。而β导数分析方法[15-17]，也是用来预测在不稳定渗流和边界控制流阶段的油气井产量简化模型，不仅可以应用在边界控制流阶段的直井，也可以应用在不稳定渗流阶段的多级压裂水平井中。

Blasingame分析方法没有考虑储层应力敏感的影响，在分析低渗透应力敏感油藏的产量变化和参数解释时，其应用受到一定的限制。本文在前人研究的基础上，建立应力敏感储层定产条件下的数学模型并进行数值求解，得到了改进的Blasingame方法，可以深入研究低渗透应力敏感油藏的产量递减规律。

1 定产模型的建立及求解

在供给半径为re、外边界封闭的圆形地层中间，一口井以恒定产量qw进行生产，见图1。

其他假设条件作如下：油气藏均匀、各向同性，上、下为不渗透边界；储层孔隙中充满了单相微可压缩流体，流动符合达西定律；忽略重力的影响；不考虑表皮效应和井筒储存效应的影响。

1.1 模型参数变化描述

1.1.1 储层渗透率变化数学描述

当储层存在应力敏感时，以生产井为中心，储层中渗透率和压力径向分布对应关系见图2。

储层渗透率随压力的变化而变化，其关系式可表示为：

图1 圆形地层中心一口直井Fig. 1 A vertical well in the center of a circular reservoir

图2 在应力敏感储层中渗透率和压力径向分布示意图Fig. 2 Radial distribution of permeability and pressure in stress sensitive reservoir

式中：α为介质应力敏感系数，1/MPa；k、k0为油气藏渗透率和初始渗透率，10-3μm2；p、pe为油气藏压力和原始压力，MPa。

1.1.2 时间的无因次处理

在应力敏感储层中，对于不同时刻、第i个空间单元的无因次时间定义为：

式中：tDi为不同时刻第i个空间单元无因次时间；ki为不同时刻第i个空间单元渗透率，10-3μm2；φ为孔隙度，%；μ为流体黏度，mPa·s；Ct为综合压缩系数，1/MPa。

在无因次时间中渗透率和时间都是变量，为消除渗透率对无因次时间的影响，可进一步定义为：

式中：tD为不同时刻的无因次时间。因此，对于不同时刻、第i个空间单元的无因次时间tDi可表示为：

式中：ωi为第i个空间单元初始渗透率与某一时刻渗透率比值，

1.1.3 压力的无因次处理

在应力敏感储层中，对于不同时刻第i个空间单元无因次压力定义为：

式中：pDi为不同时刻第i个空间单元的无因次地层压力；pe为原始地层压力，MPa；qw为油气井产量，m3/d；B为体积系数。

为消除渗透率对无因次产量的影响，可进一步定义为：

式中：pD为不同时刻的无因次压力。

因此，对于不同时刻、第i个空间单元的无因次压力pDi可表示为：

1.1.4 应力敏感系数的无因次处理

定义无因次应力敏感系数为：

1.2 应力敏感储层定产非稳态数学模型

具有应力敏感的储层在定产生产的过程中，对于不同时刻任意空间单元的无因次渗流方程如下：

特别地，对于近井筒的第1个空间单元的内边界条件为：

对于靠近封闭边界第N个空间单元的外边界条件为：

其中：i∈[1,2,…N]；rDi为第i个空间单元的无因次距离，；rD0为井壁处的无因次距离，rD0=1。

1.3 数学模型的求解

1.3.1 模型数值求解处理过程

将封闭边界油藏的整个开采过程，处理成多个在很小时间单元内和很小空间单元的渗流问题来求解。在某个时间步长内，每个空间单元体内渗透率是不变的，其渗透率大小等于每个空间单元的一侧的渗透率值。在不稳定生产的早期阶段，由于近井地带的压力降较大，为了计算结果更加精确，减小网格效应，采用不等距网格划分，在近井地带网格划分相对密集，而油气藏边界网格划分相对稀疏；整体采用点中心网格系统。网格处理见图3。

图3 不等距点中心网格Fig. 3 Non isometric point center grid

考虑压敏效应的无因次渗流方程，经过有限差分离散后得到相应的差分方程，则任意时刻第i个空间单元的差分方程为：

在任意时刻，对于近井筒的第1个空间单元的差分方程为：

在任意时刻，对于靠近封闭边界第N个空间单元的差分方程为：

由于应力敏感中的描述参数ωi为Δpi的函数，差分方程(15)、(16)、(17)构成非线性方程组。本文将参数ωi显示处理，Δpi取N时刻的无因次压力值，转化为线性方程组。

1.3.2 差分方程组求解的稳定性、收敛性分析

将差分方程(15)、(16)、(17)显示处理后得到三对角方程组。

在本文中：

2 产量递减分析方法

结合上述理论模型得到的压力解，建立理论变化图版，进行产量递减规律分析；可利用实际生产数据与无因次变量之间的正比关系，对照理论图版进行曲线拟合来解释地层参数。

2.1 产量递减分析的数据预处理

为便于对实际油气井产量的变化进行递减分析，引入参数及处理过程如下：

① 物质平衡时间tc，为目前累积产量与当前产量的比值：

提出物质平衡时间的意义，在于建立了变产量生产和定产量生产之间的关系，经过物质平衡时间处理后，变产量/压力数据就转换为定产量/压力数据。

② 定义规整化产量为：

③ 根据公式(7)，得到无因次压力倒数：

2.2 Blasingame产量递减分析方法的参数定义

Blasingame (1993)引入了规整化产量和物质平衡时间后，建立了典型递减曲线图版。无因次压力和无因次时间采用Fetkovich的递减变量定义，并在此基础上，又提出无因次规整化产量积分和无因次规整化产量积分导数[10]。处理过程如下：

①由公式(3)，Fetkovich定义的无因次时间tDd为：

②由Fetkovich定义的无因次产量qDd为：

③无因次累积产量NpDd定义为：

④无因次规整化产量积分qDdi定义为：

⑤无因次规整化产量积分导数qDdid定义为：

2.3 β导数的定义

为进一步分析无因次时间与无因次压力的在不同阶段的变化特点，定义β导数[15-17]，有：

根据以上参数定义，结合理论模型得到的压力解(考虑应力敏感影响下)，计算并编制不同条件下的产量递减复合图版，为深入分析不同条件下产量递减变化规律及参数解释奠定基础。

3 应力敏感储层产量递减规律分析

3.1 Blasingame方法图版理论曲线变化规律

当介质应力敏感系数α=0时(不考虑储层应力敏感的影响)，即为Blasingame基础分析方法。Blasingame产量递减复合图版由四簇曲线构成，见图4。

从图4中可以看出：β导数曲线，在不稳定渗流(对应无因次时间段0.001～0.1)与边界控制流阶段(对应无因次时间段0.1～100)曲线斜率明显不同，在边界控制流阶段β值趋近于1；无因次规整化产量曲线qDd，在边界控制流阶段趋近于斜率为-1的直线；在边界控制流阶段，无因次规整化产量曲线qDd与无因次规整化产量积分曲线qDdi的变化趋势出现明显差别，开始分叉并且越来越明显；无因次规整化产量积分导数曲线qDdid，整体形态像一把倒扣的“勺头”。

在不稳定渗流阶段，随着无因次半径的增加，图版中的无因次规整化产量曲线、无因次规整化产量积分曲线、无因次规整化产量积分导数曲线和β导数曲线整体向下移动；在边界控制流阶段，不同半径的每一簇曲线逐渐重合，逐渐迭合成一条曲线。

3.2 应力敏感影响下的产量递减理论变化规律

当考虑介质应力敏感的影响时，即为改进的Blasingame分析方法，其产量递减复合图版见图5，计算中采用的介质应力敏感系数α＝0.05 MPa-1。

对比图4和图5可以看出：当考虑介质应力敏感的影响时，产量递减图版曲线形态与原Blasingame符合图版大体相近，但整体变化趋势、变化幅度和变化形态存在一定的差异。其不同点在于：无因次规整化产量曲线簇、无因次规整化产量积分曲线簇、无因次规整化产量积分导数曲线簇，这三簇曲线在边界控制流阶段已经不再重合，而是互相平行；并且在每一簇曲线中，随着无因次距离的增加，对应的曲线向下方移动；而β导数曲线簇，则整体稍向上移动。其物理意义在于，当考虑介质应力敏感的影响时，导致地层压力下降速度和幅度更快，产量开始递减的时间更早。

图4 Blasingame产量递减复合图版Fig. 4 Blasingame production decline composite plate

图5 考虑应力敏感影响的产量递减复合图版Fig. 5 Production decline composite plate considering stress sensitivity

由于β导数曲线明显与其它三簇曲线在变化形态、变化趋势等方面存在明显不同，可以作为分析产量递减规律及参数解释的重要补充。

3.3 应力敏感对参数变化的影响

在相同无因次距离reD=1000条件下，不同无因次应力敏感系数对产量递减的影响，见图6～图9。

从图6～图9中可以看出：随着介质应力敏感系数的增加，无因次规整化产量曲线、无因次规整化产量积分曲线和无因次规整化产量积分导数曲线向下移动，而且这三簇曲线的弯曲程度逐渐减小；在边界控制流稳定后，不同应力敏感系数下这三条曲线基本平行。随着压敏系数的增加，β导数曲线上移，但曲线的弯曲程度减小；当边界控制流稳定后，不同压敏系数的β导数值都趋近于1。

通过以上的对比分析，在不稳定渗流(对应时间段0.001～0.1)与边界控制流阶段(对应时间段0.1～100)的过渡阶段，即在无因次时间0.1～1.0段，不同条件下各簇曲线的变化差异较大。在实际开发过程中，需要对该阶段的曲线变化进行深入分析。

图6 无因次规整化产量曲线Fig. 6 Dimensionless normalized production curve

图7 无因次规整化产量积分曲线Fig. 7 Dimensionless normalized production integral curve

图8 无因次规整化产量积分导数曲线Fig. 8 Dimensionless normalized production integral derivative curve

图9 β导数曲线Fig. 9 β derivative curve

4 实例计算对比分析

4.1 不同方法生成曲线结果对比

在封闭圆形气藏中有一口气井A，以变产量、变井底流压方式进行生产。原始地层压力为30.0 MPa，地层温度为80 ℃，地层有效厚度为10 m，孔隙度为10%，天然气相对密度为0.6，拟临界温度为195.697 K，拟临界压力为4.668 75 MPa，流体黏度为0.0228 mPa·s，岩石压缩系数为4.351 13×10-4MPa-1，原始条件下气体压缩系数为0.024 6 MPa-1，原始条件下气体体积系数为0.003 949 7。A井的生产数据(采用书《油气井现代产量递减分析方法及应用》[10]中数据)，见图10。

根据A井前期阶段(不稳定渗流阶段)的生产数据，利用Blasingame方法、本文方法计算无因次规整化产量曲线和无因次规整化产量积分曲线，分别预测后期阶段(边界控制流阶段)的曲线变化趋势，并和A井实际数据进行对比，结果见图11。

在图11中，无因次应力敏感系数为0.01，Blasingame方法处理曲线为虚线，本文方法处理曲线为实线；为了避免一条曲线在拟合过程中出现多解性，因此采用了规整化产量和规整化产量积分两条曲线。

与A井实际数据相比，在生产前期阶段(不稳定渗流阶段)，Blasingame方法、本文方法处理曲线基本重合；在生产后期阶段(边界控制流阶段)，Blasingame方法、本文方法处理曲线出现差别，考虑应力敏感的无因次曲线下降速度稍快，而本文方法处理的曲线变化趋势与A井实际更为符合。

4.2 解释结果分析

根据A井的生产数据，依据本文方法分别计算其规整化产量、规整化产量积分、规整化产量积分导数、β导数，并与无因次应力敏感系数0.01的理论图版进行拟合，如图12所示。

在实际拟合过程中，β导数曲线在拟合过程只需要左右平移，不需要上下平移。将β导数曲线与其他三条曲线画在一起，β导数曲线的纵坐标0.48对应图版的纵坐标1的位置，因此β导数曲线需要乘以系数0.48进行拟合。

根据拟合结果记录，无因次井控半径reD=50，选择任何一个拟合点，为方便起见，记录拟合点(tcaDd,qDd)M为(1, 0.1)M，相应的实际拟合点为(108, 4.8)M。

根据产量拟合点确定渗透率的大小，若用法定单位，有：

图10 A井生产曲线Fig. 10 Production curves of well A

图11 Blasingame方法和本文方法拟合曲线对比Fig. 11 Comparison of Blasingame method and the method of this paper

图12 A井拟合曲线Fig. 12 Fitting curves of well A

取应力敏感系数取二者的平均值：

将本文方法与原Blasingame方法的解释结果进行对比，其中原Blasingame图版解释过程在《油气井现代产量递减分析方法及应用》书中已详细说明，解释结果对比见表1。

表1 A井产量递减解释结果对比Table 1 Comparison of analytical results of production decline of A well

可以看出，在相同计算条件下，相对Blasingame解释方法，本文方法可以多解释出应力敏感系数这一参数，渗透率、井控半径解释结果相对较大，有效井径解释结果相对较小，更符合低渗透、应力敏感储层的生产特点。综上可看出改进Blasingame解释方法适用范围更广，结果更为准确。

5 主要认识

通过以上研究，可以得出以下认识：

(1) 在建立考虑储层应力敏感影响的定产生产模型的基础上，采用不等距差分方法和显示处理渗透率的数值解法求解模型，参考Blasingame方法，建立了改进的产量递减分析方法。

(2) 基于改进Blasingame方法，考虑储层应力敏感影响时，无因次规整化产量曲线簇，无因次规整化产量积分曲线簇，无因次规整化产量积分导数曲线簇这三簇曲线，在曲线变化趋势、曲线形态上存在较大的差别；增加的β导数曲线图版，丰富了原有分析方法。

(3) 利用应力敏感产量递减图版进行产量递减分析，对于低渗透、存在应力敏感影响的油气藏，解释结果更符合实际情况，本文方法可为此类油气藏的精细评价提供技术支持。

参考文献

[1] 李传亮. 渗透率的应力敏感效应分析方法[J]. 新疆石油地质, 2006, 27(3): 348-350.[LI C L. A study on the stress sensibility of permeability of reservoir rocks[J]. Xinjiang Petroleum Geology, 2006, 27(3): 348-350.]

[2] 郭平. 采用两种实验方法进行气藏岩心应力敏感研究[J].西南石油大学学报(自然科学版), 2009, 29(2): 7-9.[GUO P. Study on core stress sensitivity for gas reservoir with two experiment method[J]. Journal of Southwest Petroleum Institute, 2009, 29(2):7-9.]

[3] 郭小哲, 周长沙. 考虑应力敏感的页岩气藏压裂水平井渗流模型建立与分析[J]. 长江大学学报(自然科学版), 2015, 1(1): 44-50. [GUO X Z, ZHOU C S. Establishment of a seepage model considering the stress sensitivity in fractured horizon shale gas reservoirs and its analysis[J]. Journal of Yangtze University (Natural Science Edition), 2015, 1(1): 44-50.]

[4] 徐梦雅, 廖新维, 刘姣姣. 储层应力敏感性对致密气藏压裂水平井试井分析的影响[J]. 陕西科技大学学报, 2012, 30(5): 57-60. [XU M Y, LIAO X W, LIU J J. In fl uence of stress sensitivity on well test analysis for fractured horizontal well in tight gas reservoir[J]. Journal of Shanxi University of Science and technology, 2012, 30(5): 57-60.]

[5] 牛丽娟. 应力敏感性对低渗透油藏弹性产能影响[J]. 科学技术与工程, 2014, 3(14): 137-140.[NIU L J. Effect of pressure sensitivity on the elastic productivity of low permeability reservoir[J]. Science Technologies and Engineering, 2014, 3 (14): 137-140.]

[6] 田冷, 肖聪, 顾岱鸿. 考虑应力敏感与非达西效应的页岩气产能模型[J]. 天然气工业, 2014, 34(12): 70-75.[TIAN L, XIAO C, GU D H. A shale gas reservoir productivity model considering stress sensitivity and non-Darcy fl ow[J]. Natual Gas Industry, 2014, 34(12): 70-75.]

[7] 段永刚, 黄诚, 陈伟, 唐丽萍, 甯飞. 应力敏感裂缝性油藏不稳态压力动态分析[J]. 西南石油学院学报, 2001, 23(5): 19-22. [DUAN Y G, HUANG C, CHEN W, TANG L P, NING F. Pressure transient analysis for naturally fractured reservoirs with stress sensitivity[J]. Journal of southwest petroleum institute, 2001, 23(5): 19-22.]

[8] 孙来喜, 李成勇, 李成, 朱绍鹏. 低渗透气藏应力敏感与气井产量分析[J]. 天然气工业, 2009, 29(4): 74-76.[SUN L X, LI C Y, LI C, ZHU S P. Stress sensitivity effect of low-permeability gas reservoirs and production analysis of gas wells[J]. Natual Gas Industry, 2009, 29(4): 74-76.]

[9] OZKEN E, RAGHAVAN R. New solutions for well-test-analysis problems: Part 2 Computational considerations and applications[J]. SPE Formation Evaluation, 1991, 6(3): 369-378.

[10] 孙贺东. 油气井现代产量递减分析方法及应用[M]．北京: 石油工业出版社, 2013.[SUN H D. Advanced production decline analysis and application[M]. Petroleum University Press, 2013.]

[11] 程林松. 高等渗流力学[M]. 石油大学出版社, 2011.[CHENG L S. The higher seepage mechanics[M]. Petroleum University Press, 2011.]

[12] 张建国. 油气层渗流力学[M]. 石油大学出版社, 2002.[ZHANG J G. Reservoir seepage mechanics[M]. Petroleum University Press, 2002.

[13] BLASINGAME T A, LEE W J. Properties of homogeneous reservoirs, naturally fractured reservoirs, and hydraulically fractured reservoirs from decline curve analysis[C]// Society of Petroleum Engineers. SPE 15028 presented at the Permian Basin Oil and Gas Recovery Conference, 13-15 March 1986, Midland, Texas.

[14] BLASINGAME T A, LEE W J. Variable rate reservoir limits testing[C]// Society of Petroleum Engineers. SPE 15018 presented at the Permian Basin Oil and Gas Recovery Conference, 13-15 March 1986, Midland, Texas.

[15] ILK D, HOSSEINPOUR-ZONOOZI N, AMINI S, BLASINGAME T A. Application of the β-integral derivative function to production analysis[C]// Society of Petroleum Engineers. SPE 107967 presented at the 2007 SPE Rocky Mountain Oil & Gas Technology Symposium, 16-18 April 2007, Denver, Colorado.

[16] SHAHAMAT M S, MATTAR L, AGUILERA R. Analysis of decline curves based on beta derivative[C]// Society of Petroleum Engineers. SPE 169570 presented at the SPE Western North American and Rocky Mountain Joint Meeting, 17-18 April 2014, Denver, Colorado.

[17] IDORENYIN E, OKOUMA V, MATTAR L. Analysis of production data using the beta-derivative[C]// Society of Petroleum Engineers. SPE 149361 presented at the Canadian Unconventional Resources Conference, 15-17 November 2011, Calgary, Alberta, Canada.

Improved Blasingame production-decline analysis method considering stress sensitivity

CHEN Minfeng1, WANG Zhaoqi1, SUN Hedong2, WANG Yan3, CHEN Lu1, ZHANG Qichen1
1 College of Petroleum Engineering, China University of Petroleum-Beijing, Beijing 102249, China
2 Research Institute of Petroleum Exploration and Development (RIPED)-Langfang Branch, Langfang 065000, China
3 Research Institute of Exploration and Development, Xinjiang Oil fi eld Company, Karamay 834000, China

The typical Blasingame decline curve has enhanced production data analysis, which can not only account for variations in bottom-hole pressure, but also handle changing fl uid properties upon reservoir depletion more con fi dently. However, reservoir stress sensitivity is not considered by the Blasingame method. This paper introduces the seepage characteristics of stress sensitive reservoirs and builds up the mathematical model which is then solved by numerical methods at a constant-rate condition. In the solution, the non-isometric difference is used and the constant-rate condition is transformed into the inner boundary condition of the fi rst space unit, the effect of stress sensitivity on permeability is also handled through explicit methods. On this basis, production decline considering stress sensitivity was obtained according to Blasingame non-dimensional de fi nitions. The calculation results show production decline curves of stress sensitive reservoirs are parallel with each other in each group during the period of boundary control fl ow. Besides, the corresponding curve moves down with the increase of dimensionless distance. As the stress sensitive coefficient increases, improved Blasingame production decline curves all move down. Case analysis shows that the interpretation results of the improved Blasingame method are more consistent with the low permeability and stress sensitive reservoir.

stress sensitivity; non-isometric difference; improved Blasingame method; production decline curve; interpretation plate

10.3969/j.issn.2096-1693.2017.01.006

(编辑 马桂霞)

*通信作者, 416382884@qq.com

2016-08-30

国家重点基础研究发展计划(973计划)课题《致密油高效开发油藏工程理论与方法研究》资助

陈民锋, 王兆琪, 孙贺东, 王言, 陈璐, 张琪琛.考虑应力敏感影响的改进Blasingame产量递减分析方法. 石油科学通报, 2017, 01: 53-63

CHEN Minfeng, WANG Zhaoqi, SUN Hedong, WANG Yan, CHEN Lu, ZHANG Qichen. Improved Blasingame production-decline analysis method considering stress sensitivity. Petroleum Science Bulletin, 2016, 03: 53-63. doi: 10.3969/j.issn.2096-1693.2017.01.006