ABSTRACT

The Southwest China Gas Pipeline Project runs through the eastern and western regions of China, with diverse environments along the way and frequent geological disasters. Therefore, the research on the susceptibility and risk of geological disasters along the pipeline has an important guiding role in improving the safe operation of the Southwest China Gas Pipeline. Based on the Frequency Ratio (FR) and weight of evidence model (WOE), after grasping the geological hazard characteristics and distribution rules along the pipeline laying, this paper divides the geological disaster index factors in the study area and superimposes the pipeline laying factor, and then evaluates the susceptibility and risk of geological disasters in the study area. The evaluation results show that the pipeline laying can be divided into geological disaster high-prone area, medium-prone area, and low-prone area, with the area accounting for 30.38%, 49.02%, and 20.6% of the total area, respectively. After the superimposed pipeline laying method, the pipeline laying area is divided into high-risk area, medium-risk area and low-risk area of geological disasters. According to statistical analysis, the area of each dangerous area accounts for 18.2%, 44.13% and 37.67% of the total area respectively.

Keywords:Susceptibility and Risk, FR, WOE

基于FR和WOE评价模型的西南某管道敷设沿线地质灾害危险性分区研究

张博¹，石晓栊¹，高姣姣²，宗乐斌²

¹中石化川气东送天然气管道有限公司，湖北 武汉

²北京中地华安地质勘查有限公司，北京

收稿日期：2019年11月8日；录用日期：2019年11月22日；发布日期：2019年11月29日

摘 要

西南某管道工程贯穿我国东西部地区，沿途环境多样，地质灾害频发。因此，开展管道沿线的地质灾害易发性与危险性研究对于提高管道的安全运营程度具有重要的指导作用。本研究基于频率比模型和证据权模型的定量计算方法，在掌握管道敷设沿线地质灾害特征与分布规律基础上，对研究区地质灾害影响因子划分并叠加管道敷设因子后，进行研究区地质灾害易发性和危险性评价。评价结果表明，管道敷设沿线可分为地质灾害高易发区、中易发区、及低易发区，其面积分别占总面积的30.38%、49.02%、20.6%，在易发性分区基础上，叠加管道敷设方式后，将管道敷设区域分为地质灾害高危险区、中危险区、低危险区，经统计分析，各危险区的面积分别占总面积的18.2%、44.13%、37.67%。

关键词 :易发性与危险性，频率比模型，证据权模型

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

中国石化西南某管道工程是我国继西气东输之后又一条贯穿东西部地区的天然气能源大动脉，承担着将四川普光气田净化天然气向东部地区输送的重任。管道全长2366 km，途径四川、重庆、湖北、安徽、江苏、浙江、上海、江西六省二市(图1)。管道横贯我国东西，穿越地形地貌类型复杂多样，环境与地质灾害多发，威胁管道安全 [1]。特别是管道沿线山高谷深、河流水系密布的地区，地形地质条件易对管道造成不利影响。如果一旦沿线地质灾害影响到管道的运营和安全，就会直接影响到下游多家重点工矿企业的稳定生产和千家万户的正常生活，造成严重的经济损失、政治影响和社会影响。

通过开展管道沿线地质灾害危险性分区研究，对于提高管道的安全运营程度具有重要的指导作用，在实施西部大开发战略、加快四川及沿线地区经济发展、拉动国民经济增长、调整我国能源结构和充分利用天然气资源、促进我国天然气事业的大发展和国家天然气长输管网建设等方面具有重要的政治意义、社会意义和经济意义。

2. 研究区概况

此次研究在掌握管道敷设沿线地质灾害特征与分布规律基础上，开展地质灾害和地表高程、地形坡度、地形坡向、平面曲率、地层岩性、降雨等环境因子的相关关系分析，提炼区域崩塌滑坡泥石流的评价指标，为建立地质灾害预警预报指标体系提供详实的基础因子空间数据。研究基于ArcGIS进行空间管理和数据处理的基础工作，采用易于定量计算的频率比模型(FR)分析地质灾害易发性和危险性，通过证

图1. 研究区及管道地理位置图

据权模型(WOE)量化影响因子权重的优势，弥补频率比模型(FR)中主观因素对影响因子权重的影响，对研究区进行易发性和危险性评价。

3. 数据与原理

3.1. 易发性分区及危险性分区划分方法

基于易发性及危险性分区的油气管道地质灾害风险评价分析，国外的一些公司根据管道实际运行状况，已经进行将近30年的研究与应用，有着丰富的事故数据和失效处理经验，在降低管道风险上取得了一定的成效。现在已经实现由安全评价向风险评价的转变，由定性风险评价转变为半定量和定量风险评价。风险评价逐渐规范化，系统化和完整化。

意大利的SNAM公司于上世纪70年代建立相应的地质灾害监测网络，目前这个网络仍然处在更新中。该公司近年来用于地质灾害的科研经费高达数百万欧元；1997年美国对西北管道实施地质灾害风险评价，较大地提升了管道完整性管理能力；1998年加拿大对贯山管道实施地质和滑坡灾害风险评价；2005年南美管道实施了地质灾害风险评价；2002年NorAndino管道实施地质灾害风险评价，2005年对风险评价系统给予了升级；2002年安第斯(Andino)管道实施地质灾害风险评价；2003年玻利维亚0SSA-1管道实施地质灾害风险评价。实施地质灾害风险评价可以说给管道运营商带来显著的安全效益，比如南美洲的NorAndino管道实施地质灾害风险评价后，其失效概率由0.64次/(1000km.a)降到0.28次/(1000km.a)。

在对某些单点地质灾害的评价方面，国外的一些做法也值得借鉴，如加拿大AEC管道House River黄土滑坡，1998年对滑坡进行了地质评价，采用数值模拟手段评价管道受力，在此基础上提出了5种防灾减灾方案；1999年5月专门研发了基于风险评估的决策树模型(PDM)，分析防灾减灾投资与管道失效风险的耦合成本，确定减灾方案。加拿大BGG公司从20世纪90年代开始研究管道地质灾害的风险评估技术，以Kent的管道风险模型为基础，结合3S技术，开发地质灾害风险管理(GBM)系统，在实践中得到大量应用。2002年加拿大BGG公司开发的GRM系统，采用了半定量风险评价方法，对黄土湿陷、滑坡、水毁等地质灾害进行半定量风险评价，实现风险排序。

近些年来，随着GIS技术及其应用软件功能的不断强大，促进了GIS在地学领域的广泛应用和普及，使得基于GIS平台进行滑坡的空间分布规律与滑坡影响因子之间的分析已经成为了一个重要的研究方向。国内外不少学者对滑坡等地质灾害敏感因素分析和易发性评价进行了大量研究，提出了多种评价方法，主要有频率比法、神经网络法、模糊综合评判法、层次分析法、逻辑回归模型、信息量法、证据权模型等 [2] - [9]。

本课题基于ArcGIS进行空间管理和数据处理的基础工作，综合对比各种分析方法的优缺点，决定首先采用易于定量计算的频率比模型(FR)分析易发性和危险性，但单纯的采用频率比模型(FR)无法准确量化各影响因子的权重，而证据权模型(WOE)可以很好的量化影响因子权重，弥补频率比模型(FR)中主观因素对影响因子权重的影响，所以最终综合采用频率比模型(FR)和证据权模型(WOE)对研究区进行易发性和危险性评价工作。

3.1.1. 频率比模型(FR)

频率比模型(Frequency Ratio，简称FR)属于一种单变量概率分析方法，它基于每个影响因子子类中滑坡分布和每个影响因子之间的关系进行分析，通过计算每个影响因子二级子类下的滑坡面积占总滑坡面积的百分比和该影响因子二级子类下的面积占总研究区的百分比之间的比率来定量描述滑坡发育与影响因子之间的关系 [10]。计算公式为：

$FR=\frac{\frac{Np\left(L{X}_i\right)}{{\displaystyle \underset{i=1}{\overset{m}{\sum }}Np\left(L{X}_i\right)}}}{\frac{Np\left(X_j\right)}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}Np\left(X_j\right)}}}$ (1)

FR是j因子i类的概率比，Np(LX_i)是因子变量X的i类的滑坡评价单元数。Np(X_j)是因子变量X_j的评价单元数，m是因子变量X_i的二级分类因子数量。n代表因子的数量。上述公式中FR值越大，相应的因子与滑坡失稳相关性越强，反之与滑坡失稳相关性越低。

3.1.2. 证据权模型(WOE)

证据权模型(Weight of Evidence，简称WOE)是一种综合各种证据来支持一种假设的定量方法，该方法以贝叶斯概率统计为基础，最初被用于医疗诊断领域，后由加拿大数学地质学家Agterberg对其进行了修整和发展并引入地质领域，应用于成矿预测与评价中 [11]。近年来，其陆续被应用于自流井位置预测、分析地震与断层的空间关系、开采沉陷造成的岩体失稳等方面。该方法也被用于地质灾害易发性评价研究中，通过对易发性因子进行空间叠加分析，计算出各因子的权重，将各因子权重相加，就得到了研究区的滑坡易发性指数，该方法在一定程度上避免了因子权重赋值的主观性，具有较高的评价精度其计算公式如下：

$\{\begin{array}{l}{W}_i^{+}=\text{ln}\left(\frac{P\left\{B/L\right\}}{P\left\{B/\bar{L}\right\}}\right)=\ln \left(\frac{P\left(B\cap L\right)/P\left(L\right)}{P\left(B\cap \bar{L}\right)/P\left(\bar{L}\right)}\right)=\ln \left[\frac{\frac{N_{pix1}}{N_{pix1}+N_{pix2}}}{\frac{N_{pix3}}{N{}_{pix3}+N{}_{pix4}}}\right]\\ W_i^{-}=\text{ln}\left(\frac{P\left\{\bar{B}/L\right\}}{P\left\{\bar{B}/\bar{L}\right\}}\right)=\ln \left(\frac{P\left(\bar{B}\cap L\right)/P\left(L\right)}{P\left(\bar{B}\cap \bar{L}\right)/P\left(\bar{L}\right)}\right)=\ln \left[\frac{\frac{N_{pix2}}{N_{pix1}+N_{pix2}}}{\frac{N_{pix4}}{N{}_{pix3}+N{}_{pix4}}}\right]\\ W_{fi}=W_i^{+}-W_i^{-}\end{array}$ (2)

式中P表示的是某种事件发生的概率， $P\left\{B\cap L\right\}$ 表示条件概率，表示在L事件发生的条件下B事件发生

的概率，因子B表示二级因子中发生滑坡的评价单元，因子表示二级因子中未发生滑坡的评价单元，因子L表示研究区内滑坡的评价单元，因子表示研究区内未发生滑坡的评价单元，ln是自然对数。N_pix1表示该二级因子内发生滑坡的面积栅格数，N_pix2表示该二级因子外发生滑坡的面积栅格数，N_pix3表示该二级因子内未发生滑坡的栅格数，N_pix4表示该二级因子外未发生滑坡的栅格数。

$P\left\{B\cap L\right\}/P\left\{B\cap L\right\}$ 代表滑坡发生的充分率； $P\left\{B\cap L\right\}/P\left\{B\cap L\right\}$ 代表滑坡发生的必要率。 $W_i^{+}$ 表示在该二级因子影响下发生滑坡的概率，为正相关权重， $W_i^{-}$ 表示没有该二级因子影响下滑坡发生的概率，为负相关权重。当 $W_i^{+}>0$ 和 $W_i^{-}<0$ 时，表示该二级因子与滑坡的发生呈正相关，其值越大，表示其与滑坡的相关性越强、滑坡易发性程度越高；当 $W_i^{+}<0$ 和 $W_i^{-}>0$ 时，表示该二级因子与滑坡的发生呈负相关，其值负的越大，表示其与滑坡的相关性越弱、滑坡易发性程度越低； $W_i^{+}=0$ 或 $W_i^{-}=0$ 时表示该二级因子与滑坡的发生无关。W_fi表示二级影响因子的权重，W_fi越大(大于0)，表明在该二级影响因子对滑坡的发生有促进作用，反之则没有，W_fi为零时表示该二级因子与滑坡的发生无关。

3.2. 管道地质灾害的发育总体特征

本次研究的川渝、鄂西两个管理处辖区的管道线路，其中包括达化专线从四川省宣汉县普光首站至达州化肥厂的610 km和川维支线从重庆市梁平县到重庆市长寿区的158 km。输气管道沿线地质环境条件复杂，地形地貌、地层岩性、地质构造多变，加之暴雨以及人类工程活动等影响，输气管道研究区域发生的地质灾害类型多样且分布发育特征各异。

输气管道沿线地质灾害类型主要包括：不稳定斜坡、崩塌、滑坡、水工损毁、地面塌陷、泥石流等。根据2017年调查结果统计，共有灾害点274处，其中不稳定斜坡130处、崩塌57处、滑坡38处、水工损毁31处、地面塌陷17处、泥石流1处(图2(a))，根据2018年调查结果统计，共有灾害点286处，其中不稳定斜坡135处、崩塌44处、滑坡41处、水工损毁34处、地面塌陷21处、泥石流1处(图2(b))。

3.3. 灾害点面积确定及指标因子的选取

无论是用WOE模型或FR模型，都需要确定各地质灾害点的面积，以便于在Arcgis中统计各地质灾害所占栅格个数，为模型计算提供基本条件。在评价指标的统计过程中，将滑坡与不稳定斜坡统归为滑坡；滑坡和崩塌作为斜坡失稳的两种形式，其影响因素既有联系又有所差异，对滑坡起控制作用的因素同样制约着崩塌的产生，只是各影响因素所起的作用或对不同灾害类型的贡献略有差异 [12]。因此，在对管道沿线进行地质灾害易发程度定量评价区划时，将不稳定斜坡、滑坡与崩塌三种种灾害类型合并统计。

图2. 2017年(a)和2018年(b)灾害类型比例

另外，评价区泥石流仅发育1处，泥石流物源区为崩塌易发区，崩塌产生的松散物质是泥石流的主要物源，因此，本次用信息量模型评价时将泥石流灾害等效成相当规模的崩塌来考虑。因为地面塌陷与水工损毁共占地质灾害点17%，占比较少，且影响因子与不稳定斜坡和滑坡相差较大，定量分析时暂且忽略不计。

根据收集的资料、野外地质灾害调查成果分析，结合对典型灾害点的详细勘察研究，依据前人研究和沿线自然地理特征，最终选取对地质灾害影响较大的6个因子作为易发区定量评价的指标因子，主要包括地表高程、地形坡度、地形坡向、平面曲率、地层岩性、降雨 [13]。

3.3.1. 地表高程

滑坡的发育与其分布高程密切相关。一方面，由于不同的高程其地形坡度具有差异性，从而导致不同高程范围的地表集水能力的差异；另一方面，不同的高程范围其人类工程活动强度不同，使得不同高程的临空面条件具有差异性，因而高程是滑坡孕灾环境的重要因子。研究区内高程范围17~1937 m，利用ArcGIS将其分为4个区域，分别为37~500 m，500~1000 m，1000~1500 m，1500~1974 m，分别统计各个区域的栅格数，并分别提取落在各个区域内的地质灾害点，统计地质灾害点所占栅格数，然后根据分别根据FR模型和WOE模型计算出各个二级因子内所对应的FR值及W_fi值(图3(a))。

3.3.2. 地形坡度

随着坡度的增加，它对地表水径流、地下水补给和排泄、物质搬运与堆积、斜坡体应力分布特征等都具有不同程度的影响。因此，坡度是影响滑坡发生、发展及其形态特征的重要因素之一。研究区内地形坡度范围0~83˚，利用ArcGIS将其分为5个区域，分别为0˚~10˚，10˚~20˚，20˚~30˚，30˚~40˚，40˚~83˚，分别统计各个区域的栅格数，并分别提取落在各个区域内的地质灾害点，统计地质灾害点所占栅格数，然后根据分别根据FR模型和WOE模型计算出各个二级因子内所对应的FR值及W_fi值(图3(b))。

3.3.3. 地形坡向

由于不同坡向上光照、气候等物理条件的不同，会导致不同坡向上的坡体在抗风化能力上存在着差

图3. 不同指标因子的地质灾害划分

异，一般向阳坡和常年受季风影响的斜坡更容易发生滑坡等地质灾害。根据45度一个分区，便于软件计算为原则，将研究区分为8个区域，分别是0˚~45˚，45˚~90˚，90˚~135˚，135˚~180˚，180˚~225˚，225˚~270˚，270˚~315˚，315˚~360˚，分别统计各个区域的栅格数，并分别提取落在各个区域内的地质灾害点，统计地质灾害点所占栅格数，然后根据分别根据FR模型和WOE模型计算出各个二级因子内所对应的FR值及W_fi值(图3(c))。

3.3.4. 平面曲率

平面曲率是坡度或坡向在某个方向的变化率。曲率的大小从侧面反映着地面的起伏形态，反映出地面坡度的大小。曲率为正说明该处地形是向上突起的，地表形态主要是山脊地形，曲率为负说明该处地形是向下凹陷的，地表形态主要是山谷地形，曲率为零表示该处地形地势平坦，重力势能基本为零，地表没有高低起伏。

地面的高低起伏形态是由地壳在内外动力作用下长期形成的，地壳隆升和沉降都会破坏地层最开始的沉积形态，形成各种地质构造(如断层、褶皱)，在长期的物理和化学作用下会使岩层破碎，在坡面上堆积大量的碎裂岩块和碎石土，在地球内部作用下(如火山爆发、地震)使得具有一定势能的斜坡失去稳定性从而发生滑坡等地质灾害。凹坡和凸坡分布区域内滑坡发育面积差不多，但滑坡易发性程度明显不同(图3(d))。

3.3.5. 地层岩性

地层对滑坡的形成和发育起着重要的影响作用，它是孕育滑坡发生的重要内在因素，而不同类型的工程地质岩组对滑坡等地质灾害形成的影响程度有明显差异，从某种程度上来说，工程地质岩组决定了滑坡发生的类型规模特征。结合全国1:5万地质图及前期调查报告，将研究区域的岩性分为泥岩，粉砂质泥岩，泥质岩，碳酸盐岩，第四纪堆积土五类，分别统计各个区域的栅格数，并分别提取落在各个区域内的地质灾害点，统计地质灾害点所占栅格数，然后根据分别根据FR模型和WOE模型计算出各个二级因子内所对应的FR值及W_fi值(图3(e))。

3.3.6. 降雨

降雨过程中，雨水会随着岩层裂隙或破裂面在重力作用下向下渗透，会使破裂面或裂隙两边的岩体之间的摩擦阻力减小，一旦岩层裂隙或破裂面上方的岩体沿破裂面方向上向下的合力大于摩擦阻力时，这时破裂面上方的岩体将会失去原有的稳定性，从而向下滑动形成滑坡等地质灾害。结合全国2017年全国平均降雨量分布图，将研究区域的降雨量分为1000~1200 mm，1200~1600 mm，>2000 mm三个区域，分别统计各个区域的栅格数，并分别提取落在各个区域内的地质灾害点，统计地质灾害点所占栅格数，然后根据分别根据FR模型和WOE模型计算出各个二级因子内所对应的FR值及W_fi值(图3(f))。

3.4. 管道敷设

综合以上6个因素，利用FR模型和WOE模型可以叠加分析出研究区易发性分区。在易发性分区的基础上，需结合工程因素，确定对于油气管道而言的危险性分区。对于油气管道而言，能够直观反应出其危险性高低的指标就是受力大小，而不同的管道敷设方式在很大程度上影响着管道应力的分布，所以在危险性研究时引入管道敷设方式因子。

根据管道走向与滑坡主滑方向的相对关系，可以将管道敷设方式分为三类：

1) 横向敷设：管道走向与主滑方向垂直(横穿，图4(a))，管道受滑动方向的排挤，推力主要来自管经方，管道在滑坡边缘受剪切作用。

横向滑坡作用下，管道变形后的计算简图见图4(b)。

图4. 不同管道敷设方式下的应力分析

由计算简图可知，管道中点处弯矩最大，根据苏联勃洛达夫金《埋设管线》中公式求出管道中点处最大弯矩为：

$M_c=\left(M_0+\frac{q}{k^2}\right)\frac{1}{chk\frac{l}{2}}-\frac{q}{k^2}$ (3)

根据管道弯曲的微分方程及边界条件对应的挠度方程，可以利用图解法求出管道最大应力。在不做简化的情况下也可利用有限元法求出管道最大应力。

2) 纵向敷设：管道走向与主滑方向平行，管道外壁受滑坡下滑力作用，产生轴向摩擦，管道两段分别形成拉、压作用；管道在滑坡边缘收到剪切作用(图4(c))。

在计算管道应力的时候需要考虑的因素很多。综合受力分析可知，管道在与滑坡线交点A处出现弯矩，剪力最大值，是管道最先产生屈服的地方。假定花面正好通过坡脚，此时由相关理论及结构力学中的位移法，可得知管道A点处的M、Q分别为：

$M=\frac{3EIx}{l}+\frac{q{l}^2}{8}$ (4)

$Q=\frac{3EI}{l^2}x\sin \theta +\frac{5ql}{8}$ (5)

3) 斜向敷设：管道走向与主滑方向斜穿，主滑方向与管道轴向存在一定夹角，管道承受的力可分解为平行力和垂直力，可看作横敷和纵敷的组合形式。管道承受滑体的拉压和推挤，同时在滑坡边缘也存在剪切(图4(d))。

在此情况下，把平行管道轴向方向的地滑力产生的轴向应力按纵敷时轴向应力计算方法计算，然年后把垂直管道轴向的分量按横敷时的情况考虑，随后对两者进行叠加。随之可得出管道的最大应力为：

${\sigma }_t=P_{\max }/F+\left|M\right|/W$ (6)

基于三种情况下的相关应力计算方法可知管道横向敷设时最危险，纵向敷设时最不危险，斜向敷设危险性介于二者之间，故对各种敷设形式赋值如下：

1——纵向敷设，2——斜向敷设，3——横向敷设。

将管道沿线以每十公里为一个研究区域，分别统计每个区域内的危险性指数，

$Y_{管道}={\displaystyle \sum Y_{管道i}}$.

4. 地质灾害易发性分区及危险性分区

FR模型和WOE模型计算：利用arcgis栅格统计及位置统计分析工具，提取出各个二级因子的相关计算参数，根据计算结果(图5)，利用arcgis对各个二级因子图层栅格叠加计算，可得出易发性分区 [14] [15]。

图5. 二级因子模型计算结果

在此次研究中，将与地质灾害关联性较大得6个因素进行叠加分析，得到地质灾害易发性分区图(图6)。经统计分析，地质灾害高易发区、中易发区、及低易发区的面积分别占总面积的30.38%、49.02%、20.6% (表1)。

在易发性分区的基础上，地质灾害危险性指数W_危，按下式计算:

$W_{危i}=a_1{Y}_{易发i}+a_2{Y}_{管道i}$

式中，a₁为地质灾害易发性权重，取0.6；a₂为管道敷设方式权重，取0.4；

$Y_{易发i}$ 为地质灾害易发程度评价值； $Y_{管道i}$ 为管道敷设方式评价取值。利用Arcgis叠加分析可得出危险性分区结果(图7)。

在易发性分区基础上，叠加管道敷设方式，经统计分析，地质灾害高危险区、中危险区、低危险区的面积分别占总面积的18.2%、44.13%、37.67% (表2)。

图6. 易发区分布图

表1. 地质灾害易发性分区统计分析表

图7. 危险性分布图

表2. 地质灾害危险性分区统计分析表

5. 结论

1) 管道敷设沿线地质灾害影响较大的易发区定量评价的6个指标因子，主要包括地表高程、地形坡度、地形坡向、平面曲率、地层岩性、降雨。

2) 研究根据三种不同敷设方式分别进行应力计算，经计算分析可知管道横向敷设时最危险，纵向敷设时最不危险，斜向敷设危险性介于二者之间。

3) 将与地质灾害关联性较大的6个因素进行叠加分析，管道敷设沿线可分为地质灾害高易发区、中易发区、及低易发区，其面积分别占总面积的30.38%、49.02%、20.6%。

4) 在易发性分区基础上，叠加管道敷设方式后，将管道敷设区域分为地质灾害高危险区、中危险区、低危险区，经统计分析，各危险区的面积分别占总面积的18.2%、44.13%、37.67%。

文章引用

张 博,石晓栊,高姣姣,宗乐斌. 基于FR和WOE评价模型的西南某管道敷设沿线地质灾害危险性分区研究
Research on the Susceptibility and Risk Zoning of Geological Disasters along the Southwest China Gas Pipeline Based on FR and WOE Evaluation Models[J]. 地球科学前沿, 2019, 09(11): 1192-1204. https://doi.org/10.12677/AG.2019.911125

参考文献