武勝男， 張來斌， 鄧金根， 曹硯鋒， 文 敏

(1.中國石油大學(北京)機械與儲運工程學院，北京 102249；2.油氣資源與探測國家重點實驗室(中國石油大學(北京)，北京 102249；3.中海油研究總院，北京 100028)

?

基于蒙特卡洛模擬的極限注水壓力不確定性分析

武勝男1， 張來斌1， 鄧金根2， 曹硯鋒3， 文 敏3

(1.中國石油大學(北京)機械與儲運工程學院，北京 102249；2.油氣資源與探測國家重點實驗室(中國石油大學(北京)，北京 102249；3.中海油研究總院，北京 100028)

海上油田注水過程中，為防止由于注水壓力不合理而導致蓋層完整性被破壞，需要對極限注水壓力的不確定性及敏感因素對其的影響規(guī)律進行分析。以PL19-3油田某注水井為研究對象，利用FLAC3D軟件建立了儲蓋層模型，模擬計算了蓋層開始發(fā)生破壞時的注水壓力即極限注水壓力，分析了滲透率、內(nèi)聚力和內(nèi)摩擦角等參數(shù)對極限注水壓力的影響，并采用蒙特卡洛模擬方法分析了極限注水壓力的不確定性。結果表明：當注水壓力達到極限注水壓力時，蓋層開始發(fā)生破壞；隨著滲透率增大，極限注水壓力呈現(xiàn)出先增大后減小再緩慢增大的趨勢；隨著內(nèi)聚力和內(nèi)摩擦角的增大，極限注水壓力分別呈現(xiàn)線性增大和二次函數(shù)增大的趨勢。確定性極限注水壓力只是不確定性分析結果的一部分，而通過不確定性分析可以得到蓋層發(fā)生破壞可能性較大的極限注水壓力區(qū)間，能夠對極限注水壓力進行全面評估，可以為現(xiàn)場注水安全管理提供理論指導。

注水壓力；蒙特卡洛模擬；不確定性；滲透率；內(nèi)聚力；內(nèi)摩擦角

在低滲透油田注水開發(fā)過程中，通常通過提高注水壓力在地層堵塞處形成微裂縫，增大注水量，從而達到提高采收率的目的。但注水壓力不可能無限制地提高，如果注水井的井底壓力很高，注入水可能突破蓋層，竄入其他層位造成注水層間干擾，影響注水開發(fā)效果；也可能進入泥頁巖層，引起泥頁巖水化膨脹,從而擠毀套管。例如蓬萊19-3B油田注水過程中由于注入壓力過高形成地層超壓乃至破裂，導致嚴重的溢油事故，造成了嚴重的經(jīng)濟損失，同時給環(huán)境和社會帶來了巨大的危害[1]。因此，確定注水壓力界限時應充分考慮蓋層的完整性和有效的微觀封閉能力[2-3]。

一些學者開展了關于不同地質條件下注水壓力界限的研究和探討[4-5]，但對注水過程中導致蓋層破壞的注水壓力界限和地層不確定性因素對注水壓力影響的研究相對較少。傳統(tǒng)確定注水壓力界限的方法都基于地層破裂壓力，注水井井底壓力必須嚴格控制在破裂壓力以下，或取注水井地層破裂壓力的 80%～90%作為注水井最大允許井底壓力[6]，但該方法確定的注水壓力并不一定能夠確保蓋層不被破壞，并且不能評估地層不確定因素對注水壓力的影響。因此，筆者定義海上油氣田注水過程中致使蓋層開始發(fā)生破壞的注水壓力為極限注水壓力(extreme injection pressure，EIP)，即注水井可以達到最高的井底壓力，利用FLAC3D軟件建模計算極限注水壓力，采用蒙特卡洛模擬(Monte Carlo simulation，MCS)方法分析地層不確定因素對極限注水壓力的影響。

1 極限注入壓力的確定

FLAC3D是一種基于拉格朗日差分法的三維顯式有限差分程序，已經(jīng)廣泛應用于工程地質、巖土力學以及構造地質學和礦山工程領域，適宜處理大尺度、大變形工程和地質問題[7]。該軟件能夠實現(xiàn)多種材料模式的求解和不規(guī)則邊界條件的連續(xù)性問題求解，可以精確模擬材料的屈服和塑性流動狀態(tài)。注水過程中，蓋層初始發(fā)生的破壞為剪切破壞，因此，應用FLAC3D軟件研究蓋層開始發(fā)生剪切破壞時的井底壓力。

1.1 計算模型及參數(shù)的確定

1) 在分析注水井資料的基礎上，利用FLAC3D軟件提供的網(wǎng)格生成器建立儲蓋層應變模型。

2) 由于計算結果的可靠程度與巖石力學參數(shù)和原巖應力等參數(shù)有關，因此根據(jù)模擬計算區(qū)域的巖層結構和巖石力學參數(shù)確定模型的基本參數(shù)。巖體變形參數(shù)采用反映材料基本特性的體積模量和剪切模量[8]，可以將彈性模量(或變形模量)和泊松比轉化為體積模量和剪切模量，公式如下：

(1)

(2)

式中：K為體積模量，MPa；E為彈性模量，MPa；μ為泊松比；G為剪切模量，MPa。

3) 強度準則的確定。蓋層屬于彈塑性體，當荷載達到屈服強度后蓋層將發(fā)生破壞、弱化。筆者選取摩爾-庫侖準則和彈塑性本構關系作為強度準則。摩爾-庫侖準則為：

(3)

式中：fs為剪切破壞閾值，Pa；σ1和σ3分別為最大主應力和最小主應力，Pa；C為巖石的內(nèi)聚力，Pa；φ為巖石的內(nèi)摩擦角，(°)。

當fs大于0時，巖石不發(fā)生剪切破壞；當fs小于0時，巖石將發(fā)生剪切破壞。當巖石單元體達到屈服極限后，在穩(wěn)定應力下產(chǎn)生塑性變形；在拉應力狀態(tài)下，當拉應力超過巖體抗拉強度時，巖體將發(fā)生拉伸破壞。

1.2 模型邊界條件的確定

采用軟件建立儲蓋層模型時，需要確定的邊界條件為：

1) 施加初始地應力和初始孔隙壓力；

2) 儲層部位以恒定壓力注水；

3) 蓋層內(nèi)邊界為非滲透邊界，外邊界為定壓邊界；

4) 固定模型四周的法向位移。

將參數(shù)輸入模型模擬計算不同井底注水壓力時蓋層的破壞狀態(tài)，得到蓋層發(fā)生初始破壞時的井底注水壓力，即極限注水壓力。該方法確定的極限注水壓力必須不大于該區(qū)塊的地層破裂壓力。

2 不確定性分析方法

2.1 蒙特卡洛模擬

MCS方法[9]以概率統(tǒng)計理論為基礎，能準確分析各種不確定性因素的影響。該方法通過確定隨機變量的概率密度函數(shù)，從概率分布中重復抽樣建立輸出變量的分布，對輸出結果進行統(tǒng)計，并擬合輸出結果的概率分布情況。

2.1.1 概率密度函數(shù)的確定

在MCS抽樣中，實際資料統(tǒng)計表明，地質數(shù)據(jù)多數(shù)情況下呈三角形分布或正態(tài)分布[10]。三角形分布通過判斷分析確定變量的最小值a、最可能值c和最大值b來決定整個分布，目前已經(jīng)廣泛應用于概率風險模擬分析之中。服從三角形分布的概率密度函數(shù)為：

(4)

若隨機變量x服從一個數(shù)學期望為μ、方差為σ2的高斯分布，則可以用正態(tài)分布來模擬。正態(tài)分布N(μ,σ2)的概率密度函數(shù)為：

(5)

2.1.2 概率分布模擬

為了解決蓋層破壞時極限注水壓力的概率分布問題，給出一定目標區(qū)的蓋層破壞概率計算公式：

(6)

式中：P為蓋層破壞時極限注水壓力的概率；g(x)為與極限注水壓力有關的函數(shù)；f(x)為概率密度函數(shù)。

由式(6)可知，破壞概率的計算很難利用解析方法解決，因此采用MCS方法解決該類問題。從概率分布f(x)中抽出N個子樣x1，x2，…，xn，則蓋層破壞時極限注水壓力的概率為：

(7)

Crystal Ball 軟件中的MCS模塊可以用來進行各種方案的風險預測與不確定模型的研究[11]。因此，筆者利用Crystal Ball 軟件中的MCS模塊分析地層不確定因素對極限注水壓力的影響。

2.2 參數(shù)不確定性分析步驟

由于部分儲蓋層物性參數(shù)和巖石力學參數(shù)存在不確定性，因此蓋層的破壞以及極限注水壓力也存在不確定性。參數(shù)的不確定性[12]主要來源于測量誤差、采樣誤差、地質環(huán)境改變等，筆者在此主要討論計算極限注水壓力時所涉及到的參數(shù)，并進行不確定性研究。不確定性分析步驟為：

1) 分析儲蓋層的巖石力學參數(shù)的不確定性，確定對計算極限注水壓力影響較大的巖石參數(shù)。

2) 分析不確定性參數(shù)對極限注水壓力的影響規(guī)律，通過多項式擬合建立極限注水壓力與不同參數(shù)之間的函數(shù)關系式。利用決定系數(shù)R2表征擬合方程的擬合程度，R2越接近1，表明極限注水壓力與不同參數(shù)之間的擬合程度越高。

3) 確定不同巖石參數(shù)的概率密度函數(shù)，并開展參數(shù)的不確定性分析。

4) 確定模擬次數(shù)，運用Crystal Ball 軟件進行模擬計算。

5) 統(tǒng)計計算結果，繪制概率分布圖，預測不同參數(shù)下極限注水壓力的概率分布區(qū)間。

3 應用實例

選取蓬萊19-3油田某注水井作為研究對象，其注水層位為L40—L70油組。蓬萊19-3油田注水井井深1 100.00～1 300.00 m，不考慮各層的差異性，儲蓋層參數(shù)取該油田一般值。該油田上覆巖層壓力為24.0 MPa，原始地層孔隙壓力為11.0 MPa，最大水平主應力為22.2 MPa，最小水平主應力為16.5 MPa。

根據(jù)PL19-3油田某井現(xiàn)場鉆井完井報告、油田地層物性及巖石力學分析結果，假設蓋層是完整的(沒有原生裂縫或斷層)，計算突破蓋層時的注水壓力。儲蓋層的參數(shù)見表1。

表1 PL19-3油田某井儲蓋層的參數(shù)

Table 1 Reservoir-caprock parameters of a well in PL19-3 Oilfield

參數(shù)儲層蓋層單軸抗壓強度/MPa3．850抗拉強度/MPa0．320．42彈性模量/GPa46內(nèi)聚力/MPa1．231．62孔隙度0．240．10滲透率/mD170010內(nèi)摩擦角/(°)2424

3.1 極限注水壓力的計算

3.1.1 模型的建立

利用FLAC3D軟件建立儲蓋層模型如圖1所示，模型尺寸(長、高、厚)50.0 m×0.2 m×30.0 m，包括上蓋層、儲層、下蓋層各10個單元。不考慮井眼影響，儲層部位以恒定壓力注水，蓋層內(nèi)邊界為非滲透邊界，外邊界為定壓邊界，同時固定模型四周的法向位移。

3.1.2 模擬結果分析

圖2為通過模擬計算得到的注水時蓋層發(fā)生剪切破壞的情況(圖中藍色區(qū)域代表儲蓋層未發(fā)生剪切破壞，粉色區(qū)域代表儲蓋層已經(jīng)發(fā)生破壞)。

圖1 儲蓋層模型Fig.1 Caprock model

圖2 不同井底注水壓力下蓋層剪切破壞狀態(tài)Fig.2 Shear failure state of cap for different bottom hole injection pressure

從圖2可以看出：當井底注水壓力為21.2 MPa時，儲層區(qū)域發(fā)生破壞，而蓋層均未發(fā)生破壞；當井底注水壓力為21.3 MPa時，儲層和蓋層的接觸面開始發(fā)生破壞。這表明，隨著井底注水壓力增大，蓋層由最開始的不發(fā)生破壞到破壞區(qū)域逐漸增大；若再繼續(xù)提高注水壓力，蓋層區(qū)域將進一步被破壞，可能出現(xiàn)突破蓋層的風險。從圖2還可以看出，極限注水壓力為21.3 MPa，由于該油田的地層破裂壓裂為24.0 MPa，滿足極限注水壓力小于地層破裂壓力的要求。

3.2 不確定性分析

實際注水過程中，極限注水壓力的求取受到地層不確定性參數(shù)的影響，這些不確定參數(shù)包括單軸抗壓強度、抗拉強度、彈性模量、內(nèi)聚力、內(nèi)摩擦角、滲透率和孔隙度等地層參數(shù)。這些參數(shù)隨地質構造在某一范圍內(nèi)隨機變化，很難用具體數(shù)值表征，且部分參數(shù)之間具有相關性，比如內(nèi)聚力受到單軸抗壓強度、抗拉強度、彈性模量等影響，滲透率與孔隙度相關，因此需要開展不確定因素的敏感性分析。筆者在假設其他參數(shù)不變的情況下，分析了內(nèi)聚力、內(nèi)摩擦角和滲透率對極限注水壓力的影響。

3.2.1 不同參數(shù)對極限注水壓力的影響

通過解釋測井資料得知該油田儲層滲透率的變化范圍為500～3 000 mD，內(nèi)聚力變化范圍為1.4～3.0 MPa，內(nèi)摩擦角變化范圍為20°～36°。假設其他參數(shù)不變，在上述范圍內(nèi)分析不同變量的取值對極限注水壓力的影響。

根據(jù)模擬結果，繪制極限注水壓力與滲透率、內(nèi)聚力和內(nèi)摩擦角的關系曲線，如圖3所示。

圖3 極限注水壓力與滲透率、內(nèi)聚力和內(nèi)摩擦角的關系曲線Fig.3 Relationship of extreme injection pressure vs. permeability,cohesive force and internal friction angle respectively

從圖3可以看出：隨著地層滲透率增大，極限注水壓力呈先增大后減小，再緩慢上升最后趨于穩(wěn)定的趨勢；隨著內(nèi)聚力增大，極限注水壓力呈線性增大趨勢；隨著內(nèi)摩擦角增大，極限注水壓力增大。

對圖3中極限注水壓力與滲透率、內(nèi)聚力和內(nèi)摩擦角的關系進行擬合，得到：

10-6K2+3.8×10-3K+19.883

(8)

R2=0.978 9

pEIP=4.358 3C+14.412

(9)

R2=0.999 3

pEIP=-0.014 5φ2+1.100 5φ+3.254 8 (10)

R2=0.991 0

式(8)—式(10)的決定系數(shù)R2都大于0.970 0，表明這3個擬合關系式的擬合程度很高，可以利用這3個關系式預測該油田不同滲透率、內(nèi)聚力和內(nèi)摩擦角下的極限注水壓力，以降低注水過程中發(fā)生突破蓋層的風險。

3.2.2 參數(shù)不確定性分析

極限注水壓力不確定性分析是根據(jù)滲透率、內(nèi)聚力和內(nèi)摩擦角的概率密度函數(shù)隨機取值，帶入式(8)—式(10)中求得極限注水壓力，基于概率統(tǒng)計理論，應用MCS抽樣計算得到具有一定概率分布的計算結果。通過測井資料解釋可知，儲層滲透率的最大值為3 000 mD，最小值為500 mD，最有可能出現(xiàn)的值為1 700 mD，則認為滲透率的概率密度函數(shù)服從三角分布。同理，內(nèi)聚力的概率密度函數(shù)服從區(qū)間(1.4 MPa，1.6 MPa，3.0 MPa)的三角分布。內(nèi)摩擦角隨著地層的抗拉強度的變化而變化，假定其概率密度函數(shù)服從(30°，3°)的正態(tài)分布。確定滲透率、內(nèi)聚力和內(nèi)摩擦角的分布后，把極限注水壓力定義為預測參數(shù)，采用直接抽樣法設定抽樣次數(shù)為10 000次，通過模擬計算統(tǒng)計極限注水壓力的概率分布，結果見圖4。

極限注水壓力計算結果為21.3 MPa，與圖4中滲透率、內(nèi)聚力和內(nèi)摩擦角所對應的概率分別為0.004 93，0.034 91和0.023 77，均在置信水平5%以下，說明蓋層發(fā)生破壞的可能性比較小。由圖4可知，置信水平5%～95%所對應的極限注水壓力區(qū)間分別為(21.51 MPa，22.46 MPa)、(21.47 MPa，26.26 MPa)和(21.73 MPa，24.00 MPa)，即蓋層發(fā)生破壞可能性較大的極限注水壓力分布在上述區(qū)間。鑒于地層參數(shù)的不確定性，極限注水壓力計算結果只是不確定性分析中的一部分，根據(jù)該結果很難準確評估蓋層發(fā)生破壞的可能性，然而不確定性分析卻能夠對極限注水壓力進行全面評估，可以為現(xiàn)場注水安全管理提供指導。

圖4 滲透率、內(nèi)聚力和內(nèi)摩擦角對應極限注水壓力的概率分布Fig.4 Probability distribution of extreme injection pressure corresponding to permeability, cohesive force and internal friction angle

4 結論及建議

1) 海上油田注水過程中蓋層是否發(fā)生破壞與極限注水壓力密切相關，當注水壓力達到極限注水壓力時，蓋層開始發(fā)生破壞，若再繼續(xù)提高注水壓力，蓋層將會被進一步破壞，可能出現(xiàn)突破蓋層的風險。

2) 極限注水壓力計算結果只是不確定性分析中的一部分，根據(jù)其很難準確評估蓋層發(fā)生破壞的可能性，而不確定性分析可以得到蓋層發(fā)生破壞可能性較大的極限注水壓力區(qū)間，能夠對極限注水壓力進行全面評估，可以為現(xiàn)場注水安全管理提供理論指導。

3) 分析極限注水壓力時，由于缺少大量的地質數(shù)據(jù)資料，只考慮了滲透率、內(nèi)聚力和內(nèi)摩擦角3個參數(shù)，進行了簡單的分析。如果能獲得詳細的地質數(shù)據(jù)資料，可以擬合出更加符合實際的參數(shù)分布概率密度函數(shù)，進行更全面、準確的分析。

4) 建議將考慮蓋層破壞的基于蒙特卡洛模擬分析極限注水壓力不確定性的方法，推廣到考慮斷層開啟、套管破壞的注水模型中，為海上油田注水提供更有價值的決策支持。

[1] WEI Lai，HU Zhuowei，DONG Lin，et al.A damage assessment model of oil spill accident combining historical data and satellite remote sensing information：a case study in Penglai 19-3 oil spill accident of China[J].Marine Pollution Bulletin，2015，91(1)：258-271.

[2] JIN Zhijun，YUAN Yusong，SUN Dongsheng.Models for dynamic evaluation of mudstone/shale cap rocks and their applications in the Lower Paleozoic sequences，Sichuan Basin，SW China [J].Marine and Petroleum Geology，2014，49:121-128.

[3] 黃學，付廣，賴勇，等.斷裂對蓋層封閉性破壞程度定量研究[J].大慶石油地質與開發(fā)，2008，27(6)：5-9.

HUANG Xue，F(xiàn)U Guang，LAI Yong，et al.Quantitative research on destroy degree of caprock seal by faults[J].Petroleum Geology & Oilfield Development in Daqing，2008，27(6)：5-9.

[4] 曹廣勝，譚暢，宋福昌，等.低滲透油田強化注水壓力界限研究[J].科學技術與工程，2012，12(13)：3107-3110，3115.

CAO Guangsheng，TAN Chang，SONG Fuchang，et al.Pressure limit of enhanced water injection in low permeability oilfield[J].Science Technology and Engineering，2012，12(13)：3107-3110，3115.

[5] 趙永勝，陸蔚剛，蘭玉波，等.多層砂巖油田水驅開發(fā)的合理注水壓力[J].大慶石油地質與開發(fā)，2000，19(6)：21-24.

ZHAO Yongsheng，LU Weigang，LAN Yubo，et al.Reasonable water injection pressure of water flooding of multi-layer sandstone oil fields[J].Petroleum Geology & Oilfield Development inDaqing，2000，19(6)：21-24.

[6] 劉洪杰，陳鴻，王海更，等.海上復雜斷塊油田合理注水壓力研究[J].斷塊油氣田，2015，22(2)：228-233.

LIU Hongjie，CHEN Hong，WANG Haigeng，et al.Study on reasonable water injection pressure of offshore complex fault-block oilfield[J].Fault-Block Oil and Gas Field，2015，22(2)：228-233.

[7] 張書偉，林杭，任連偉.FLAC3D在巖土工程中的應用[M].北京：中國水利水電出版社，2011：83-107.

ZHANG Shuwei，LIN Hang，REN Lianwei.The application of FLAC3Dto geotechnical engineering[M].Beijing：China Waterpower Press，2011：83-107.

[8] 孫國權，李娟，胡杏保.基于FLAC3D程序的采空區(qū)穩(wěn)定性分析[J].金屬礦山，2007(2)：29-32.

SUN Guoquan，LI Juan，HU Xingbao.FLAC3D-based stability analysis of mined-out area[J].Metal Mine，2007(2)：29-32.

[9] 李春艷.利用蒙特卡洛模擬與風險分析優(yōu)選水平井完井方式[J].石油鉆探技術，2009，37(5)：52-57.

LI Chunyan.Applications of Monte Carlo simulation method and risk analysis in horizontal well completion optimization[J].Petroleum Drilling Techniques，2009，37(5)：52-57.

[10] 董玉亮，顧煜炯，楊昆.基于蒙特卡洛模擬的發(fā)電廠設備重要度分析[J].中國電機工程學報，2003，23(8)：201-205.

DONG Yuliang，GU Yujiong，YANG Kun.Criticality analysis on equipment in power plant based on Monte Carlo simulation[J].Proceedings of the Chinese Society for Electrical Engineering，2003，23(8)：201-205.

[11] 李帥芳，肖果平.Crystal Ball在項目管理風險分析中的應用[J].項目管理技術，2013，11(4)：40-44.

LI Shuaifang，XIAO Guoping.Application of Crystal Ball in project management risk analysis[J].Project Management Technology，2013，11(4)：40-44.

[12] HENLEY H，LUCIA A.Constant pressure Gibbs ensemble Monte Carlo simulations for the prediction of structure I gas hydrate occupancy[J].Journal of Natural Gas Science and Engineering，2015，26：446-452.

[編輯 劉文臣]

Monte Carlo Simulation-Based Uncertainty Analysis on Extreme Water Injection Pressure

WU Shengnan1, ZHANG Laibin1, DENG Jingen2, CAO Yanfeng3, WEN Min3

(1.College of Mechanical and Transportation Engineering, China University of Petroleum(Beijing), Beijing, 102249, China；2.State Key Laboratory of Petroleum Resource & Prospecting (China University of Petroleum(Beijing)), Beijing, 102249, China;3.CNOOC Research Institute, Beijing, 100028, China)

In order to avoid caprock integrity failure caused by unreasonable water injection pressure during water injection in offshore oil fields, it is necessary to conduct an uncertainty analysis on extreme injection pressure and the controlling influences of sensitive factors on extreme injection pressure. In this paper, a water injector in the PL19-3 Oilfield was taken as a research object and a caprock model was established by using FLAC3Dsoftware. Based on this model, the extreme water injection pressure (i.e. the water injection pressure when the caprock is initially broken) was calculated. Then, the effect of sensitive parameters (e.g. permeability, cohesive force and internal friction angle) on extreme injection pressure was analyzed. Furthermore, the uncertainty of extreme water injection pressure was studied by means of Monte Carlo simulation (MCS). It is shown that cap rocks begin to break when the extreme injection pressure is reached. As the permeability increases, extreme injection pressure varies, first in the trend of increasing, decreasing then increasing slowly. The extreme injection pressure increases linearly as the cohesive force increases, but in the form of quadratic function with the increasing of the internal friction angle. A deterministic result of extreme injection pressure is only a part of the uncertainty analysis. Through uncertainty analysis, the range of extreme injection pressure with higher probability of caprock failure can be revealed and extreme injection pressure can be evaluated comprehensively so as to provide the theoretical guidance in the secure management of water injection.

injection pressure; Monte Carlo simulation; uncertainty; permeability; cohesive force; internal friction angle

2015-07-13；改回日期：2016-02-14。

武勝男(1986—)，女，山西大同人，2010年畢業(yè)于青島科技大學安全工程專業(yè)，2013年獲中國石油大學(北京)安全技術及工程專業(yè)碩士學位，在讀博士研究生，研究方向為海上工藝過程風險與海洋裝備的可靠性。E-mail：wshn19@163.com。

中海石油(中國)有限公司科研攻關項目“注水井堵塞機理及增注降耗方法研究-解堵注水安全評價”(編號：CNOOC-KJ125ZDXM06LTD-08-ZY-12)部分研究內(nèi)容。

?油氣開發(fā)?

10.11911/syztjs.201603

TE921

A

1001-0890(2016)03-0109-06