吴超，刘建华，张东清，陈小峰，赵伟杰

（1.中国石油化工股份有限公司石油工程技术研究院；2.中国石油大学（北京）石油天然气工程学院）

基于地震波阻抗的预探井随钻井壁稳定预测

吴超1，刘建华1，张东清1，陈小峰1，赵伟杰2

（1.中国石油化工股份有限公司石油工程技术研究院；2.中国石油大学（北京）石油天然气工程学院）

常规井壁稳定预测方法应用于预探井时，由于资料缺乏及使用条件限制、计算参数误差、运算过程复杂等因素会严重影响其预测效果，针对这些问题提出了一种适用于预探井的随钻井壁稳定预测方法。通过考察波阻抗与孔隙压力、构造地应力、岩石强度之间的定量关系，建立包含波阻抗与井壁稳定力学参数的非线性模型，在此基础上利用神经网络算法识别波阻抗与地层三压力（孔隙压力、坍塌压力与破裂压力）之间的映射关系。通过神经网络分层建模及录井资料实时分析，实际钻进中利用地震反演波阻抗数据随钻预测钻头前方的井壁稳定性，及时优化钻井液密度以控制井壁稳定性。现场预探井的应用结果表明，新方法较之常规技术有效提高了适应性与运算速度，操作流程更为便捷，预测精度能满足工程需求。图5表1参21

预探井；井壁稳定性；随钻预测；地震波阻抗；神经网络算法

0 引言

钻井工程中由于井壁失稳引起的复杂故障会造成严重的钻井质量与安全事故，准确描述井壁围岩力学环境、深入探索井壁失稳机理、建立合理的井壁稳定计算模型，是实现安全、优质、高效钻井的关键。目前井壁稳定钻后定量检测方法已基本成熟[1-5]，它通过钻井与测井资料的综合分析，确定井壁失稳的主导机理，计算并建立地层孔隙压力、坍塌压力及破裂压力剖面，指导后续钻井。但从实际工程角度而言，要求对井壁失稳问题尽可能做到提前预知、及早发现、尽快处理，所以井壁稳定性预测研究意义更为重大。当前井壁稳定预测方法总体包括两类，首先提出的是钻前预测方法，它利用地震特征参数分析、地震反演等方法预测得到待钻井全井段的声波速度等参数，利用其预测井壁稳定性并用于钻井设计[6-9]。对探井而言，由于钻前对其地质情况了解有限，钻前预测精度常不够理想，需要在实际钻进中根据实钻情况对钻前预测模型进行实时调整与更新，井壁稳定随钻预测方法因此发展起来[10-12]。这类方法基于地震属性与测井数据之间的相关性，综合利用地震、测井、录井及钻井资料，对钻遇的每一层系的井壁稳定性进行精细分析评价，实钻中及时修正钻前预测数据。常规随钻预测方法虽然有效提高了预测精度及效率，但也存在一些问题：对于资料相对缺乏的预探井，其井壁稳定计算模型参数存在较大误差，多步运算过程中的误差累计会严重降低最终预测结果的准确度；另一方面对于待预测的预探井尤其是区域探井，其与参照井可能位于不同的物探测网，这也会影响常规预测方法的应用效果。

针对上述情况，本文提出了一种新的井壁稳定随钻预测方法，该方法根据波阻抗与岩石力学参数、孔隙压力、地应力等井壁稳定参数之间存在的映射关系，综合利用地震反演波阻抗、智能计算工具及实钻信息进行井壁稳定随钻预测。

1 波阻抗与井壁稳定参数的关系

1.1 概述

地震反演是以地震勘探资料为基础，以已知地质特征以及钻井、测井资料为约束，对地下岩层空间的岩石物理参数进行成像的过程。地震资料中包含丰富的岩性、构造、物性信息，经过地震反演可以把地震信息转换为岩石物理信息，建立起适用于勘探开发的地质、力学和储集层模型。

按所使用的地震资料来划分，地震反演主要包括叠前反演和叠后反演两大类。叠前反演主要包括基于旅行时的层析成像技术和基于振幅的AVO分析技术；叠后反演主要指基于振幅的波阻抗及其他岩石物理参数反演技术。其中波阻抗反演技术发展得最为成熟，具有严密的物理和数学理论基础[13-14]，较好地将反演精度与运算效率结合起来，反演结果集地震剖面的横向连续性、测井资料垂向高分辨率以及地质剖面的直观性于一体，综合来看最适合钻井工程师使用。

波阻抗是重要的地层特征参数，为地震波在地层中传播速度和岩石密度的乘积。以目前最常用的叠后纵波阻抗反演为例，纵波阻抗公式为：

由公式（1）可见，地层波阻抗与岩石的性质，包括岩性、岩石物理、岩石力学、埋藏深度、地质年代、构造特征、含油气性等均关系密切。因此，特定的波阻抗分布规律包含丰富的地层信息，能不同程度地反映地层的岩石力学特征，这为利用波阻抗预测井壁稳定性提供了依据。

岩石力学理论表明，导致井壁垮塌和破裂的根本原因是井筒围岩发生剪切和拉伸破坏，岩石强度参数、地层孔隙压力、构造地应力是决定井壁稳定性的关键地质参数。当以波阻抗数据为基础预测井壁稳定性时，必须充分考察井壁稳定参数与波阻抗之间的关系，以反映孔隙压力、构造应力及岩石强度（包括粘聚力、内摩擦角、抗拉强度等）等参数对井壁稳定性的影响，在此基础上优选井壁稳定预测模型。

1.2 波阻抗与地层孔隙压力之间关系

岩石物理实验表明，纵波速度与有效应力之间存在如下关系[15]：

根据岩石物理实验，泥质含量和纵、横波速度之间存在以下经验关系[16]：

综合考察以上各式可见，地层孔隙压力与地震波速度、密度之间存在一定函数关系，波阻抗是地震波速度与密度的乘积，且地震波速度与密度之间也存在一定相关性。考虑到地震反演波阻抗相比分别反演速度和密度可信度更高，因此建立孔隙压力与波阻抗之间的映射关系模型更适于预测孔隙压力及井壁稳定性，即

1.3 波阻抗与构造应力之间关系

勘探开发实践表明，地震波在岩层中的传播速度整体上随地质过程中构造作用力的增大而增大，即波阻抗与构造运动关系密切。目前国内外已提出多种构造地应力计算模型，其适应性取决于目标探区的地质情况。以其中使用最为广泛的黄氏模型为例[17]：

由弹性动力学理论可知，岩石弹性参数与地震波速度及密度关系密切[17]：

综合考察以上各式可见，构造应力与地震波速度、岩石密度之间存在一定的函数关系，可建立构造地应力与波阻抗之间的映射关系，即：

1.4 波阻抗与岩石强度之间关系

通过岩石力学实验数据分析，可以建立适用于钻井工程的多种岩石强度经验计算模型，例如岩石抗压和抗拉强度可以使用以下经验公式计算[17]：

而基于实验资料统计得到的岩石粘聚力与内摩擦角的常用计算公式[18]为：

综合考察以上各式，认为岩石强度与波阻抗之间存在一定的映射关系，再考虑到岩石强度受围压影响较大，结合上文中构造应力计算模型，可以建立包含岩石强度、波阻抗等因素的关系模型：

2 井壁稳定随钻预测模型

2.1 预测原理

按照通常的研究思路，在通过地震反演得到待钻井的波阻抗数据之后，分别利用以上各种理论或经验公式求取孔隙压力、构造应力、粘聚力、内摩擦角、抗拉强度等数据，再利用剪切及拉伸破坏准则分别求取岩石坍塌压力与破裂压力，通过地层三压力（孔隙压力、坍塌压力与破裂压力）综合分析，对井壁稳定性进行预测。但对于预探井而言，由于探区前期完钻井数量少，用于统计分析的测井、测试、地质等资料较为匮乏，导致常规预测模型中的大量系数和经验参数存在较大误差；同时，常规井壁稳定预测方法计算步骤繁多，致使误差逐步累积，最终计算结果的预测精度也将会受到很大影响。其次，由于预探井地质条件复杂，常规方法中复杂的运算流程易导致操作出错。此外，由于地震属性分析的需要，常规预测方法要求参照井和待预测井必须位于同一地震测网内，对于预探井尤其是区域探井而言，其与已完钻的参照井常具有不同的观测系统，这也会严重影响常规方法的预测效果，所以需研发新方法用于预探井的井壁稳定预测。

从另一个角度来看，波阻抗和上覆压力分别与孔隙压力、构造应力、岩石强度之间构成隐函数关系，进一步根据井壁稳定力学理论，波阻抗和坍塌压力、破裂压力之间也存在一定的非线性映射关系，即

由于钻井地质环境的复杂性和测量手段的限制，这种函数关系具有非线性、多解性、模糊性、随机性等特征，因此直接确定或应用这种非线性关系存在一定困难。可使用智能计算工具处理这类问题，其特点是无需直接确定非线性映射模型的系数，而通过智能算法自动识别各类地质参数之间的关联性。目前在理论和应用上都比较成熟的智能算法有神经网络、遗传算法、蚁群算法等，针对预探井井壁稳定预测问题的实际要求，本文采用神经网络算法实现这一过程，即首先利用完钻井资料进行神经网络学习建模，再利用其预测目标井的井壁稳定性。

2.2 操作流程

2.2.1 钻前阶段工作

①收集、整理与分析已完钻参照井的测井、钻井、地质、测试等资料，对于本区已有完钻井的预探井，最好选取同区邻井作为参照井；对于同区无完钻井的区域探井，则选取邻近区块具有代表性的完钻探井作为参照井。②以测井数据为基础，结合其他资料，求取参照井的波阻抗、上覆压力、孔隙压力、坍塌压力及破裂压力数据。需要说明的是，在对区域探井进行预测时，由于资料有限，利用参照井密度测井资料建立的上覆压力梯度模型也用于待预测井。③利用神经网络算法，分层系建立参照井各层段波阻抗、上覆压力、孔隙压力、坍塌压力、破裂压力的网络模型，以供实钻计算使用。④收集待预测井所在工区的地震资料，利用广义线性算法反演得到待预测井的波阻抗数据，以供实钻计算使用。

2.2.2 实钻过程中的随钻预测工作

①实时检测分析待预测井的录井资料，确定当前钻遇地层的类型。②从已建立的分层神经网络模型中选择对应当前钻遇地层类型的模型。③将待预测井钻头前方一定深度范围内的波阻抗、上覆压力数据输入神经网络，经过计算可以从网络输出待钻地层的孔隙压力、坍塌压力及破裂压力，进而预测待钻地层的井壁稳定性。在以上步骤中，最关键的是建立参照井的神经网络模型并利用其计算待预测井的地层三压力。下面重点论述这一过程。

神经网络算法种类很多，根据井壁稳定预测过程具有先学习再应用的特点，同时随钻预测操作要求学习速度快、学习精度高，这种情况下较适宜采用有监督学习的网络。有监督学习的神经网络主要包括BP网络、径向基函数网络、小波网络等。笔者结合实际问题对上述几种常用算法进行了反复试算比较，结果表明在利用神经网络建立波阻抗与地层压力之间的非线性映射关系模型时，运用径向基函数学习算法能在保证较高预测精度的同时，还具有较快的计算速度，因此在进行井壁稳定随钻预测时使用径向基函数神经网络效率最高。除了神经网络，笔者在研究中也尝试应用了其他一些主流的智能算法，如遗传算法、蚁群算法、粒子群算法，应用结果表明对于本文提出的随钻预测模型，以上算法的综合性能均不如径向基函数神经网络，主要表现为各类算法在计算精度方面差别不大，但遗传算法、蚁群算法等算法的操作简易度、运算快捷性均逊于径向基函数算法。

笔者使用的径向基函数神经网络是以一种函数逼近理论为基础的3层前向网络，网络的隐层特性函数采用非线性径向基函数，通过非线性变换将输入向量映射到新的高维空间，输出层对隐层输入进行线性加权叠加。本文使用的径向基函数神经网络的输入层有2个节点，分别为波阻抗与上覆地层压力，输出层为3个节点，分别为待预测的坍塌压力、孔隙压力与破裂压力。

通常选取径向基函数为高斯函数，其网络输出为：

径向基函数网络的学习过程包括隐层单元学习和输出层单元学习两个阶段，第1阶段通常采用无监督的自组织聚类算法进行学习，关键是确定基函数的中心。在本阶段学习中，首先给出径向基函数的初始中心向量ck(0)，以ck(l)为中心，利用最邻近划分准则将xi划归m个聚类区。

令ck(l+1)为内各xi的中心，则有：

计算编码误差参数d(l)为：

判断d(l)是否小于给定的门限误差，若小于误差则此时的ck(l)即为基函数的中心向量。得到聚类中心之后，通常根据各中心的距离分布按经验方法确定宽度。

确定径向基函数的中心与宽度以后就进入第2阶段学习，其目的是要确定隐含层到输出层的各权值ωkj。与BP神经网络算法基本一致，本阶段是利用样本训练集按常规的梯度下降法进行有监督的学习而得到ωkj，具体算法可参阅文献[19-21]。在钻前阶段，需要完成网络学习各阶段的工作，实际钻进过程中将钻头前方待钻地层的地震反演波阻抗及上覆压力数据输入已经学习好的神经网络，网络即在输出端产生对应地层所需预测的坍塌压力、孔隙压力和破裂压力。

3 现场应用

SA区块位于西非地区，探井深度一般在3 000 m以内，井壁失稳问题是该地区钻井施工中的突出问题，尤其是Anguille、Azile等层段井壁垮塌掉块现象比较普遍，同时地层压力较高，井控压力较大。WZ井位于SA5探区，是本工区的第1口探井，因此选取邻近SA3工区中的重点井IE井作为参照井。在WZ井开钻前，首先利用SA5探区的地震资料进行波阻抗反演，图1为该探区的波阻抗反演结果（以过WZ井的二维剖面为例），图2为从波阻抗反演数据体中提取的WZ井波阻抗曲线。

图1 SA5探区地震波阻抗反演成果

钻前阶段的工作还包括对参照井资料的计算分析。图3为利用测井资料得到的IE井波阻抗曲线，同时根据IE井的多种资料计算出该井的坍塌压力、孔隙压力和破裂压力（见图4）。此外根据IE井的密度测井资料得到本区块上覆地层压力经验计算公式：

以本区块RPG、UPC、LPC、Anguille、Azile等地层为基础，按照工程地质特征进一步将其划分为更多的次级层段，作为随钻预测的基本单元。通过径向基函数神经网络学习算法，分别建立IE井波阻抗与上覆压力和坍塌压力、孔隙压力及破裂压力之间的分层神经网络模型。按照实钻录井信息的指示，在WZ井实际钻进到各层段顶部时，分别将该井各层段的波阻抗反演数据输入相应层段的神经网络映射模型，分别预测出WZ井各层段的坍塌压力、孔隙压力和破裂压力（见图5，为节约篇幅，将分别预测得到的各段地层压力曲线组合在一起进行显示）。

图2 从地震波阻抗反演数据中提取的WZ井波阻抗曲线

图3 从测井数据中提取的IE井波阻抗曲线

图4 基于测井数据计算的IE井孔隙压力、坍塌压力及破裂压力曲线

图5 利用波阻抗反演数据预测的WZ井孔隙压力、坍塌压力及破裂压力曲线

在WZ井实际钻进过程中，根据随钻预测的地层三压力曲线进行井壁稳定随钻评价，对钻前设计的钻井液密度进行实时调整优化以满足控制井壁失稳的要求。WZ井完钻后井径测井资料显示，Anguille、Azile等重点层段的井径扩大率仅为7.35%，与邻近探区的完钻井数据对比，井径扩大率减小59.17%，说明基于压力预测成果的钻井液密度优化是有效的。为进一步验证预测精度，在WZ井完钻之后收集了该井现场试验及测试数据，并将地层孔隙压力与破裂压力预测值和实测值进行了对比（见表1），其中破裂压力实测值通过地层破裂压力实验得到，孔隙压力实测值通过地层压力测试得到。由表1可见各测试点压力预测精度均在90%以上，可以满足工程施工的要求，再次证明本文提出的预测方法可行。

表1 WZ井孔隙压力和破裂压力预测值与实测值的对比

4 结论

波阻抗与孔隙压力、地应力、岩石强度乃至坍塌压力及破裂压力之间存在较为复杂的非线性映射关系，可利用神经网络等智能计算工具识别这种关系，在此基础上利用波阻抗反演数据预测井壁稳定性。实际应用中，在钻前进行地震波阻抗反演和神经网络分层学习建模，实钻中利用建立的神经网络分层模型分层预测钻头下方待钻地层的井壁稳定性，可以及时优化钻井液密度以控制井壁稳定性。该方法运算过程简捷、经验参数少，可以减少误差累积，较快的运算速度有利于实现随钻预测，在现场应用中能够达到较高的预测精度，还克服了传统预测方法对地震资料使用条件的限制，适用于资料相对缺乏的预探井尤其是区域探井。

符号注释：

Ip——纵波阻抗，106kg/（m2·s）；vp——纵波速度，km/s；ρ——岩石密度，g/cm3；Vsh——泥质含量，f；σv——上覆地层压力，MPa；ppor——地层孔隙压力，MPa；a1，a2，a3，A1，A2，B0，B1，B2，B3，B4，C1，C2，C3，V1，V2，V3，u1，u2，ξ1，ξ2——模型系数，无量纲；vs——横波速度，km/s；σH，σh——水平方向的最大与最小地应力，MPa；E——弹性模量，MPa；μ——泊松比，无量纲；α——有效应力系数，无量纲；σc——抗压强度，MPa；St——抗拉强度，MPa；C——粘聚力，(MPa)2；φ——内摩擦角，（°）；M——构造应力，MPa；S——岩石强度参数，MPa；pc——坍塌压力，MPa；pf——破裂压力，MPa；yji——对应深度待预测的地层3压力值，MPa；ωkj——隐层到输出层的连接权值，无量纲；xi——参照井的波阻抗与上覆压力数据组成的输入向量；ck，σk——径向基函数的中心和宽度，均为2维向量；Ωk(l) ——聚类区；d(l) ——误差，无量纲；D(l) ——计算编码，无量纲；M(l) ——Ωk(l)中包含的xi数目，无量纲；i——样本编号；n——学习样本对的数目；j——输出层节点编号；k——连接权，无量纲；m——基函数的个数；l——递推步骤数；Go——上覆地层压力梯度当量钻井液密度，g/cm3；h——井深，m；f1，f2，f3，F——表示波阻抗与井壁稳定性参数之间关系的函数，无量纲。

[1] Karstad E,Aadnoy B S.Optimization of borehole stability using 3D stress optimizations[R].SPE 97149,2005.

[2] Hemphill T.Integrated management of the safe operating window:Wellbore stability is more than just fluid density[R].SPE 94732,2005.

[3] Greenwood J A,Brehm A,Van Oort E,et al.Application of real-time wellbore stability monitoring on a deepwater ERD well[R].SPE 92588,2005.

[4] Nes O M,Fjar E,Tronvoll J,et al.Drilling time reduction through integrated rock mechanics analysis[R].SPE 92531,2005.

[5] Bradford I D R,Aldred W A,Cook J M,et al.When rock mechanics met drilling:Effective implementation of real-time wellbore stability control[R].SPE 59121,2000.

[6] 金衍,陈勉,张旭东.钻前井壁稳定预测方法的研究[J].石油学报,2001,22(3):96-99.Jin Yan,Chen Mian,Zhang Xudong.Study on prediction for borehole stability before drilling[J].Acta Petrolei Sinica,2001,22(3):96-99.

[7] 金衍,陈勉,杨小奇.利用层速度钻前预测安全泥浆密度窗口[J].岩石力学与工程学报,2004,23(14):2430-2433.Jin Yan,Chen Mian,Yang Xiaoqi.Prediction of safe mud density window by interval velocity before drilling[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(14):2430-2433.

[8] 金衍,陈勉,吴超.探井二开以下地层井壁稳定性钻前预测方法[J].石油勘探与开发,2008,35(6):742-745.Jin Yan,Chen Mian,Wu Chao.Borehole stability prediction before drilling for formations under the second section in an exploration well[J].Petroleum Exploration and Development,2008,35(6):742-745.

[9] 金衍,陈勉.利用地震记录钻前预测井壁稳定性研究[J].石油学报,2004,25(1):89-92.Jin Yan,Chen Mian.Prediction of borehole stability by seismic records[J].Acta Petrolei Sinica,2004,25(1):89-92.

[10] 吴超,陈勉,金衍.井壁稳定实时预测新方法[J].石油勘探与开发,2008,35(1):80-84.Wu Chao,Chen Mian,Jin Yan.Real-time prediction method of borehole stability[J].Petroleum Exploration and Development,2008,35(1):80-84.

[11] 吴超,陈勉,金衍.基于地震属性技术的井壁稳定随钻预测新方法[J].中国石油大学学报:自然科学版,2007,31(6):141-146.Wu Chao,Chen Mian,Jin Yan.The prediction model of borehole stability based on seismic attribute technology[J].Journal of China University of Petroleum:Edition of Natural Science,2007,31(6):141-146.

[12] Wu Chao,Chen Mian,Jin Yan.A prediction method of bore hole stability based on seismic attribute technology[J].Journal of Petroleum Science and Engineering,2009,65(3/4):208-216.

[13] 李庆忠.论地震约束反演的策略[J].石油地球物理勘探,1998,33(4):423-438.Li Qingzhong.On the strategy of seismic restricted inversion[J].Oil Geophysical Prospecting,1998,33(4):423-438.

[14] Treitel S,Lines L.Past,present and future of geophysical inversion:A new millennium analysis[J].Geophysics,2001,66(1):21-24.

[15] 樊洪海.地层孔隙压力预测检测新方法研究[D].北京:石油大学,2001.Fan Honghai.Study on calculation and prediction methods of formation pore pressure[D].Beijing:University of Petroleum,2001.

[16] 金衍.井壁稳定预测理论及应用研究[D].北京:石油大学,2001.Jin Yan.Study on predicting method of borehole stability and its application[D].Beijing:University of Petroleum,2001.

[17] 陈勉,金衍,张广清.石油工程岩石力学[M].北京:科学出版社,2008.Chen Mian,Jin Yan,Zhang Guangqing.Rock mechanics in petroleum engineering[M].Beijing:Science Press,2008.

[18] 楼一珊,金业权.岩石力学与石油工程[M].北京:石油工业出版社,2006.Lou Yishan,Jin Yequan.Rock mechanics and petroleum engineering[M].Beijing:Petroleum Industry Press,2006.

[19] Elanayar S,Shin Y C.Radial basis function neural network for approximation and estimate on of nonlinear stochastic dynamic system[J].IEEE Transactions on Neural Network,1994,4(5):594-604.

[20] Bors A G,Pitas I.Median radial basis function neural network[J].IEEE Transactions on Neural Network,1996,7(6):1351-1364.

[21] Robert J,Schilling J,Caroll J.Approximation of nonlinear systems with radial basis function neural network[J].IEEE Transactions on Neural Network,2001,12(1):21-28.

（编辑 黄昌武）

A prediction of borehole stability while drilling preliminary prospecting wells based on seismic impedance

Wu Chao1,Liu Jianhua1,Zhang Dongqing1,Chen Xiaofeng1,Zhao Weijie2
(1.SINOPEC Research Institute of Petroleum Engineering,Beijing 100101,China;2.Faculty of Petroleum Engineening,China University of Petroleum,Beijing 102249,China)

When the conventional method of predicting borehole stability is applied to the preliminary prospecting wells,the precision of prediction is often affected severely by the lack of data,limit of applying conditions,error of calculated parameters and complexity of operating courses.This paper presents a method of predicting borehole stability while drilling preliminary prospecting wells and establishes a nonlinear model including wave impedance and borehole stability mechanical parameters by investigating the quantitative relationships between wave impedance and pore pressure,in-situ stress as well as rock strength.Based on the nonlinear model,neural network algorithm is used to identify the relationship between the wave impedance and the three formation pressures (pore pressure,fracture pressure and collapse pressure).Through the establishment of the model by layered neural network and the timely analysis of logging data,the wellbore stability before the bit can be predicted while drilling by using the data of seismic wave impedance.Field application in preliminary prospecting wells shows that the new method has higher adaptability,faster calculation speed,and simpler operational procedure than conventional methods,and its prediction accuracy can meet the requirement of engineering.

preliminary prospecting well;borehole stability;prediction while drilling;wave impedance;neural network

国家科技重大专项“海相油气井井筒环境监测技术”（2011ZX05005-006-001）

TE22

A

1000-0747(2015)03-0390-06

10.11698/PED.2015.03.18

吴超（1976-），男，安徽明光人，博士，中国石化石油工程技术研究院高级工程师，主要从事油气井岩石力学分析、钻井工程优化设计等方面的研究工作。地址：北京市朝阳区北辰东路8号北辰时代大厦916室，邮政编码：100101。E-mail：wuchao9138@163.com

2014-10-22

2015-04-08