葛新民，范卓穎，范宜仁，李潤澤，徐擁軍

?

基于核主成分?小波能譜分析的復(fù)雜儲層油水界面預(yù)測

葛新民1, 2，范卓穎1，范宜仁1, 2，李潤澤3，徐擁軍1, 2

(1. 中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東青島，266580；2. 中國石油大學(xué) CNPC測井重點實驗室，山東青島，266580；3. 中國石油青海油田天然氣開發(fā)公司，青海格爾木，816000)

針對復(fù)雜儲層油水界面確定困難的問題，研究基于核主成分?小波能譜分析的確定方法，分析油水界面特征并確定最優(yōu)小波基函數(shù)。將深感應(yīng)電阻率、中子、密度、聲波和自然伽馬等常規(guī)曲線作為輸入信號，應(yīng)用核主成分方法提取第一核主成分信號進(jìn)行小波變換得到近似和細(xì)節(jié)信號后分別計算不同頻帶的能譜，基于相鄰兩級能譜差異最大原理確定單井的油水界面。在單井分析的基礎(chǔ)上確定整個油藏的油水界面。研究結(jié)果表明：Coiflet3基函數(shù)在研究工區(qū)中應(yīng)用效果最好；將第一核主成分進(jìn)行6級分解，并根據(jù)第5級和第6級能譜最大差法確定的油水界面與生產(chǎn)測試資料十分吻合；油水界面高度由西向東呈先增大后減小的趨勢，與儲層巖性、物性的平面變化特征一致。

油水界面；核主成分；能譜分析；小波基函數(shù)；Coiflet3基函數(shù)

油水界面對油氣藏流體識別、儲量計算及油氣田井網(wǎng)布置、施工方案設(shè)計等均具有重要意義。復(fù)雜儲層的巖性、物性和孔隙結(jié)構(gòu)變化大，不同區(qū)域發(fā)育特征各異，流體分異差異大，給油水界面識別帶來極大困難。實際油藏并不存在油水截然分開的界面，由于重力分異作用，油藏總體呈現(xiàn)上油下水的趨勢，地質(zhì)工作者和石油公司總是希望鉆遇更多油層且避開水層。油水界面確定可分為地震、測井、壓力和地化法等。地震法方面，劉傳虎等[1]分析了三維地震資料屬性，利用油層、水層在振幅能量上的明顯差異性識別原始油水界面；王學(xué)忠等[2]應(yīng)用油水層地震振幅差異進(jìn)行了春光油田排2油藏的油水界面識別；Berle等[3]應(yīng)用4D地震資料和隨鉆電阻率確定了Troll油田的油水界面。測井法方面，Stashin[4]應(yīng)用完井、生產(chǎn)動態(tài)和測井資料得到了Utikumua油田的原始及現(xiàn)今油水界面；魯國明等[5]提出無需提供油藏條件的油水界面張力和接觸角的測井解釋與毛管壓力相結(jié)合的油水界面確定方法；馬勇等[6]應(yīng)用聲波和三孔隙度差比值等確定了塔河南油田凝析油氣藏的氣水界面；李靜等[7]用自然電位、電阻率和中子伽馬分析了鄂爾多斯盆地中部南梁油田N64井區(qū)的油水界面；趙軍龍等[8]提出基于時間推移測井的油水界面確定方法。壓力法方面，Harris等[9]提出了基于孔、滲、飽和毛管壓力的油水界面擬合公式；肖忠祥等[10]利用三次樣條插值將MDT所得壓力進(jìn)行等間距采樣得到連續(xù)的地層壓力梯度線，進(jìn)而提取油水界面；李興麗等[11]應(yīng)用擬毛管壓力、含油飽和度和相對滲透率進(jìn)行了油水界面的預(yù)測。地化法方面，張春明等[12]將熱解分析與氣相色譜技術(shù)結(jié)合確定了塔里木盆地YH7X-1井寒武系碳酸鹽巖儲層的油水界面；姜振學(xué)等[13]應(yīng)用包裹體顆粒指數(shù)和定量顆粒熒光技術(shù)對塔北典型油氣藏的油水界面進(jìn)行了恢復(fù)。無論是地震、測井還是地化等法，油水界面識別的核心在于敏感信息的提取。隨著勘探開發(fā)的不斷深入，常規(guī)方法在復(fù)雜儲層油水界面識別中存在越來越多局限性[14?15]。本文作者將小波能譜分析引入油水界面預(yù)測中，通過核主成分分析得到流體最佳表征信息，對第一核主成分信號進(jìn)行小波分解及能譜計算，從不同頻帶能譜數(shù)據(jù)的差異中提取出合適的油水界面。

1 小波變換及能譜分析

小波變換是在傅里葉變換基礎(chǔ)上發(fā)展起來的一種信號分析方法，通過可變的時域窗口實現(xiàn)信號多分辨率分析，能提供信號各頻帶信息，也稱為自適應(yīng)時?頻局部化分析方法。它繼承了傅里葉變換的時?頻局部化特性，克服了傅里葉變換在局部時間段頻域特征信息提取困難的問題，已廣泛應(yīng)用于地震和測井信號處理、流體識別及沉積旋回分析等石油勘探領(lǐng)域。連續(xù)小波變換可表示為[16]：

式中：()為輸入連續(xù)信號，為子小波；為尺度參數(shù)；為位移參數(shù)。

測井?dāng)?shù)據(jù)是隨深度變化的離散點，對尺度參數(shù)和位移參數(shù)離散后即可得到離散小波變換。令尺度參數(shù)為，位移參數(shù)為，可得離散尺度函數(shù)和小波函數(shù)[17]：

Mallat算法是離散小波變換的有效方法之一[18]，應(yīng)用高通和低通濾波器對信號進(jìn)行分解得到細(xì)節(jié)和近似信號。圖1所示為對信號進(jìn)行三級離散小波變換的示意圖。若離散信號的長度為(=0, 1, 2, …,?1)，則離散小波變換的近似和細(xì)節(jié)信號可寫為：

對測井信息進(jìn)行離散小波變換得到近似和細(xì)節(jié)信號后，可計算不同頻帶的能量：

式中：E(=1, 2, …,)為不同頻帶的能量；為信號分解的最佳尺度；d為第個頻帶的第細(xì)節(jié)信息；為每個頻帶的細(xì)節(jié)信息個數(shù)。

每個頻帶的能量所占信號比例(即能量譜)可寫為

式中：()為原始離散信號；為離散信號的長度；E為第個信號在第個頻帶上的能量譜。

應(yīng)用小波能量譜識別油水界面的基本原理為：油、水層在小波能量譜中存在突變性，當(dāng)相鄰兩尺度間的能量譜變化最大時，則認(rèn)為該處是油水界面[19]。

圖1 一維離散信號三級小波分解

2 主成分分析及核映射

對于復(fù)雜儲層，通過單一測井信息難以有效表征流體響應(yīng)。主成分分析是一種良好的數(shù)據(jù)融合方法，通過降維技術(shù)將各變量間互相關(guān)聯(lián)的復(fù)雜關(guān)系簡化，進(jìn)而提取若干表征流體信息的綜合信號。將個測井參數(shù)記為1，2，…，X，將其轉(zhuǎn)換成某種線性組合的新指標(biāo)1，2，…，F(≤)，這些新指標(biāo)互相獨立，可表示為[20]

系數(shù)矩陣滿足：

各主成分之間相互獨立且無重疊信息，即

主成分方差依次遞減，即

主成分分析的本質(zhì)是對特征向量進(jìn)行線性映射，而測井參數(shù)由于測量原理等不同，在本質(zhì)上存在明顯的非線性，通過主成分分析難以得到有效、綜合性強的參數(shù)，甚至可能導(dǎo)致有用信息丟失。核主成分法將特征向量進(jìn)行非線性映射，能很好地解決該問題。若是原空間的一個訓(xùn)練樣本集合，是維向量，為原空間到高維Hilbert空間的映射函數(shù)，則有[21]

高維Hilbert空間的協(xié)方差矩陣為

協(xié)方差矩陣符合特征方程：

特征向量由特征空間數(shù)據(jù)集組成：

將式(17)和式(16)代入式(15)可得

式中：為系數(shù)；為Gram矩陣，記為，。

由于Gram矩陣是半正定的，對Gram矩陣進(jìn)行分解得到特征向量后用式(14)計算核主成分。任意樣本的非線性特征可表示為

核主成分分析的核心是核函數(shù)的選擇，較為常用的核函數(shù)有多項式函數(shù)、徑向基函數(shù)(高斯函數(shù))和Sigmoid函數(shù)[21?22]：

實際應(yīng)用中需要將核矩陣中心化再進(jìn)行主成分提取，核矩陣中心化方程為

將深感應(yīng)電阻率、中子、密度、聲波和自然伽馬曲線作為輸入信號，核主成分分析中采用徑向基函數(shù)作為核函數(shù)，通過多次試算確定為10。圖2所示為對我國東部油田X-1井分別進(jìn)行主成分和核主成分分析所得的主成分貢獻(xiàn)率分布。從圖2可知：核主成分分析的效果十分明顯，第一核主成分信號所占貢獻(xiàn)率為0.66，而第一主成分信號所占貢獻(xiàn)率僅為0.37。

(a) 主成分貢獻(xiàn)率；(b) 核主成分貢獻(xiàn)率

3 油水界面確定方法及特征分析

對測井?dāng)?shù)據(jù)進(jìn)行核主成分分析后，提取第一核主成分信號即可進(jìn)行小波能譜分析。X-1井測井曲線及第一核主成分信號圖如圖3所示。圖3中第5道是對X-1井的測井曲線進(jìn)行分析后提取的第一核主成分。

圖3 X-1井測井曲線及第一核主成分信號圖

從圖3可看出：第一核主成分在形態(tài)上既有電阻率曲線的總體平滑信息，也有三孔隙度和自然伽馬曲線的局部鋸齒狀信息。與單一測井信息相比，能較好地表征流體性質(zhì)。

小波基函數(shù)對小波能譜分析的效果影響嚴(yán)重，一般需要考察其連續(xù)性、對稱性、正交性、消失矩、線性相位、時頻窗口的中心和半徑以及時頻窗的面積等性質(zhì)[23]。表1所示為幾種常見的小波基函數(shù)的基本性質(zhì)[24]。從表1可知：Coiflets，Symlets和Daibechies基函數(shù)具有較完整的性質(zhì)。

表1 常見的幾種小波基函數(shù)及其性質(zhì)

注：√表示具有該性質(zhì)；~√表示近似具有該性質(zhì)；—表示不具有該性質(zhì)；為小波函數(shù)的消失矩；r和d分別為重構(gòu)和分解小波的消失矩。

分別對Coiflets，Symlets和Daibechies基進(jìn)行實驗優(yōu)選，根據(jù)小波分解后的重構(gòu)信號與真實信號的誤差確定最優(yōu)小波基函數(shù)[25]。經(jīng)研究，選用Coiflet3作為基函數(shù)并對第一核主成分信號進(jìn)行6級分解即可達(dá)到較好效果。圖4所示為對X-1井進(jìn)行6級小波分解所得的小波能譜(只顯示了4，5，6級能譜)。從圖4可知：當(dāng)界面高度為2 288 m時，5級能量譜與6級能量譜差異最大，將該處定義為油水界面。對該井段進(jìn)行分段試油，射開2 288.5~2 292 m時，油管出水，無油；射開2 284~2 288.5 m時，油管為油水同出，該層日產(chǎn)油140 m3，日產(chǎn)水287 m3，小波分析確定的油水界面與試油結(jié)論基本一致。對研究區(qū)解釋為油水同層的92個層段進(jìn)行小波分析確定了油水界面，對44個具有試油結(jié)論的層位進(jìn)行統(tǒng)計分析，有39個層的小波分析結(jié)果與試油結(jié)論基本一致，符合率達(dá)88.6%?？梢?，應(yīng)用核主成分?小波能譜分析所確定的油水界面精度較高。

圖4 X-1井核主成分信號的三級能量譜(4, 5, 6級能譜)

圖5所示為該油藏西—東方向1條剖面的油水界面分布示意圖，圖6所示為該剖面的平面圖。從圖5可知：自西向東，油水界面高度先增大后減小的趨勢。根據(jù)巖性及物性資料測試結(jié)果，自東向西，儲層泥質(zhì)含量呈先增大后減小，孔隙度、滲透率呈先減小后增大的特征。根據(jù)壓汞毛管壓力可知，當(dāng)巖石的物性差、泥質(zhì)含量高時其排驅(qū)壓力大，油水界面高度增大。根據(jù)小波能譜分析預(yù)測的油水界面與巖性、物性分布具有一致性[26]。

圖5 油水界面分布示意圖(西—東剖面)

圖6 井位平面分布圖

4 結(jié)論

1) 核主成分對儲層流體綜合表征參數(shù)的提取有重要意義。第一核主成分信號與第一主成分信號相比，貢獻(xiàn)率有較大提高(由0.37提高至0.66)。

2) Coiflet3小波基在小波能譜分析中的效果最好。通過核主成分?小波能譜分析后，第5級和第6級能譜差異最大處即為油水界面。油水界面預(yù)測結(jié)果與生產(chǎn)測試資料吻合度高，符合率達(dá)88.6%。

3) 研究工區(qū)的油水界面高度自西向東先增大后減小，與巖性、物性分布規(guī)律一致。自西向東，儲層泥質(zhì)含量呈先增大后減小的趨勢，孔隙度和滲透率呈先減小后增大的趨勢。物性變差，排驅(qū)壓力增大，油水界面高度增大。

致謝：

感謝中國石油大學(xué)勝利學(xué)院朱學(xué)娟老師和實驗室碩士研究生牟洋洋提供的幫助。

[1] 劉傳虎, 王學(xué)忠, 喬明全. 高精度三維地震屬性識別淺層稀油油藏油水界面[J]. 天然氣勘探與開發(fā), 2009, 32(4): 24?27. LIU Chuanhu, WANG Xuezhong, QIAO Mingquan. Apply high-precision 3D seismic attribute to identifying oil-water interface of shallow and thin oil reservoir[J]. Natural Gas Exploration and Development, 2009, 32(4): 24?27.

[2] 王學(xué)忠, 劉傳虎, 王建勇. 應(yīng)用三維地震屬性識別春光油田原始油水界面[J]. 特種油氣藏, 2009, 16(3): 47?48, 52. WANG Xuezhong, HU Chuanhu, WANG Jianyong. Identify original oil-water interface by 3D seismic attributes in Chunguang Oilfield[J]. Special Oil and Gas Reservoirs, 2009, 16(3): 47?48, 52.

[3] Berle A E, Thorsen A K, Gjengedal J A, et al. Determining the oil-water contact while drilling above the transition zone in long horizontal wells-combining 4d seismic and advanced resistivity geosteering for improved wellbore placement[C]// SPWLA 51st Annual Logging Symposium. Perth: SPWLA, 2010: 82870.

[4] Stashin K. An analytical approach to determining oil/water contact rise at Utikuma field[C]//40th Annual Technical Meeting of the Petroleum Society of CIM. Banff: Petroleum Society of Canada Press, 1989: 40?59.

[5] 魯國明, 朱家俊. 利用毛細(xì)管壓力預(yù)測油水界面位置的一種回歸分析方法[J]. 石油勘探與開發(fā), 1995, 22(4): 70?73.LU Guoming, ZHU Jiajun. A regression analytical method to estimate oil-water contact using capillary pressure data[J]. Petroleum Exploration And Development, 1995, 22(4): 70?73.

[6] 馬勇, 趙毅, 章成廣. 塔河南油田氣層識別與氣油界面劃分[J]. 天然氣工業(yè), 2008, 28(6): 59?61. MA Yong, ZHAO Yi, ZHANG Chengguang. Identification of gas zone and demarcation of oil/gas interface in the Southern Tahe Oilfield[J]. Natural Gas Industry, 2008, 28(6): 59?61.

[7] 李靜, 李江江. 油水界面定性識別測井技術(shù)研究[J]. 延安大學(xué)學(xué)報(自然科學(xué)版), 2011, 30(2): 49?51, 53.LI Jing, LI Jiangjiang. Research on qualitative analysis oil-water interface logged identifying[J]. Journal of Yanan University (Natural Science Edition), 2011, 30(2): 49?51, 53.

[8] 趙軍龍, 張寧, 閆博, 等. 基于時間推移的測井信息確定油水界面變化方法研究[J]. 地球物理學(xué)進(jìn)展, 2013, 28(1): 293?300. ZHAO Junlong, ZHANG Ning, YAN Bo, et al. Study on the identifying the oil/water boundary alteration with logging data based on the time-lapse[J]. Progress In Geophysics, 2013, 28(1): 293?300.

[9] Harris T G, Luffel D L, Hawkins J M. Method for determining hydrocarbon/water contact level for oil and gas wells: United States Patent, US005621169A[P]. 1997?04?15.

[10] 肖忠祥, 徐銳. 基于MDT測試資料的儲層油水界面自動提取方法研究[J]. 地球物理學(xué)進(jìn)展, 2008, 23(6): 1872?1877. XIAO Zhongxiang, XU Rui. Study of the method to extract oil-water interface automatically from MDT data[J]. Progress in Geophysics, 2008, 23(6): 1872?1877.

[11] 李興麗, 王培春, 呂洪志, 等. 改進(jìn)擬毛細(xì)管壓力曲線方法及其在預(yù)測油水界面中的應(yīng)用[J]. 石油天然氣學(xué)報, 2012, 34(10): 68?71.LI Xingli, WANG Peichun, Lü Hongzhi, et al. An improved method for obtaining pseudo-capillary curves by using logging- interpreted permeability and its application in predicting oil- water contact[J]. Journal of Oil and Gas Technology, 2012, 34(10): 68?71.

[12] 張春明, 方孝林, 朱俊章. 用熱解和氣相色譜技術(shù)確定碳酸鹽巖儲集層油水界面[J]. 石油勘探與開發(fā), 1998, 25(2): 24?26. ZHANG Chunming, FANG Xiaolin, ZHU Junzhang. Determining oil-water contact in carbonate reservoir using GC and rock eval[J]. Petroleum Exploration and Development, 1998, 25(2): 24?26.

[13] 姜振學(xué), 王顯東, 龐雄奇, 等. 塔北地區(qū)志留系典型油氣藏古油水界面恢復(fù)[J]. 地球科學(xué): 中國地質(zhì)大學(xué)學(xué)報, 2006, 31(2): 201?207. JIANG Zhenxue, WANG Xiandong, PANG Xiongqi, et al. Paleohydrocarbon water contact restoration of typical silurian oil and gas reservoirs in the Northern Tarinl Basin[J]. Earth Science: Journal of China University of Geosciences, 2006, 31(2): 201?207.

[14] 胡曉慶, 趙鵬飛, 武靜. 確定油水界面的新方法: 剩余壓力法[J]. 中國海上油氣, 2011, 23(1): 40?42. HU Xiaoqing, ZHAO Pengfei, WU Jing. A new method to determine oil-water contact: Excess pressure[J]. China Offshore Oil and Gas, 2011, 23(1): 40?42.

[15] Khajooie S, Qassamipour M, Farhani M. Application of a new method to estimate free fluid level and recognition of compositional grading using wireline formation testing data[C]// 6th International Petroleum Technology Conference. Richardson: Society of Petroleum Engineers, 2013: 1057?1066.

[16] Daubechies I. Ten lectures on wavelet[M]. Philadelphia: Capital City Press, 1992: 35?60.

[17] 邸繼征. 小波分析原理[M]. 北京: 科學(xué)出版社, 2010: 12?20. DI Jizheng. Wavelet analysis theory[M]. Beijing: Science Press, 2010: 12?20.

[18] Mallat S G, Hwang W L. Singularity detection and processing with wavelets[J]. IEEE Transactions on Information Theory, 1992, 38(2): 617?642.

[19] Heidary M, Javaherian A. Wavelet analysis in determination of reservoir fluid contacts[J]. Computers & Geosciences, 2013, 52(1): 60?67.

[20] 吳今培. 基于核函數(shù)的主成分分析及應(yīng)用[J]. 系統(tǒng)工程, 2005, 23(2): 117?120. WU Jinpei. Principle component analysis and application based on kernel function[J]. Systems Engineering, 2005, 23(2): 117?120.

[21] 張學(xué)工. 模式識別[M]. 北京: 清華大學(xué)出版社, 2010: 25?40. ZHANG Xuegong. Pattern recognition[M]. Beijing: Tsinghua University Press, 2010: 25?40.

[22] 劉素京, 楊琳, 王從慶. 基于核主成分分析和支持向量機的飛機艙音信號識別[J]. 東南大學(xué)學(xué)報(自然科學(xué)版), 2008, 38(增刊2): 123?127. LIU Sujing, YANG Lin, WANG Congqing. Aircraft cockpit voice signal recognition based on kernel principal component analysis and SVM[J]. Journal of Southeast University (Natural Science Edition), 2008, 38(Suppl 2): 123?127.

[23] 汪新凡. 小波基選擇及其優(yōu)化[J]. 株洲工學(xué)院學(xué)報, 2003, 17(5): 33?35. WANG Xinfan. Selection and optimization of wavelet base[J]. Journal of Zhuzhou Institute of Technology, 2003, 17(5): 33?35.

[24] 曾憲偉, 趙衛(wèi)明, 師海闊, 等. 利用小波包變換對地震信號進(jìn)行時頻分析時小波基函數(shù)的選取[J]. 地震研究, 2010, 33(4): 323?328.ZENG Xianwei, ZHAO Weiming, SHI Haikuo, et al. Selection of wavelet basis function in process of time-frequency analysis of earthquake signals using wavelet packet transform[J]. Journal of Seismological Research, 2010, 33(4): 323?328.

[25] 張華, 陳小宏, 楊海燕. 地震信號去噪的最優(yōu)小波基選取方法[J]. 石油地球物理勘探, 2011, 46(1): 70?75. ZHANG Hua, CHEN Xiaohong, YANG Haiyan. Optimistic wavelet basis selection in seismic signal noise elimination[J]. Oil Geophysical Prospecting, 2011, 46(1): 70?75.

[26] 林景曄, 童英, 王新江. 大慶長垣砂巖儲層構(gòu)造油藏油水界面控制因素研究[J]. 中國石油勘探, 2007, 12(3): 13?16. LIN Jingye, TONG Ying, WANG Xinjiang. Research on controlling factors of oil-water interface of structural oil reservoirs in sandstone reservoirs of Daqing Placanticline Structure[J]. China Petroleum Exploration, 2007, 12(3): 13?16.

Oil/water contact prediction of complex reservoir using kernel principal component analysis and wavelet power spectrum analysis

GE Xinmin1, 2, FAN Zhuoying1, FAN Yiren1, 2, LI Runze3, XU Yongjun1, 2

(1. College of Geosciences, China University of Petroleum, Qingdao 266580, China; 2. CNPC Key Laboratory of Logging, China University of Petroleum, Qingdao 266580, China;3. Qinghai Oilfield Natural Gas Development Company, CNPC, Golmud 816000, China)

To solve the problem of oil-water contact determination for complex reservoir, a new method based on kernel principal component and wavelet power spectrum analysis was carried out. The optimal wavelet base function was analyzed and chosen. The deep induction resistivity, neutron, density, acoustic and natural gamma data were used as the input signals, the kernel principal component (KPCA) method was applied to extract the first kernel principal component signal, and then wavelet transformation was exacted to get the approximate and detail signals. Different frequency band power spectrums were calculated based on the decomposed approximate and detail signals respectively. The oil-water contact of single well was determined based on the principle of the indication of maximum difference between adjacent power spectrum. The oil-water contact distribution of the whole reservoir was determined based on the analysis of single well. The results show that Coiflet3 base function is best applied in the study area; Applying 6 levels of decomposition to the first principal component is the best fit for the determination of the oil-water contact according to the biggest difference method of grade 5 and grade 6 power spectrum value. The result is in good agreement with that of the production and test. The oil-water contact first increases and then decreases from west to east and is consistent with the planar variation characteristics of reservoir lithology and physical properties.

oil water contact; kernel principal component; power spectrum analysis; wavelet base function; Coiflet3 base function

10.11817/j.issn.1672-7207.2015.05.024

P631.8+4

A

1672?7207(2015)05?1747?07

2014?06?30；

2014?09?22

國家自然科學(xué)基金資助項目(41404086)；國家科技重大專項(2011ZX05020-008)；中國博士后基金資助項目(2014M560591)；山東省自然科學(xué)基金資助項目(ZR2014DQ007) (Project(41404086) supported by the National Natural Science Foundation of China; Project(2011ZX05020-008) supported by the National Science and Technology Major Program of China; Project(2014M560591) supported by the China Postdoctoral Science Foundation; Project(ZR2014DQ007) supported by the Natural Science Foundation of Shandong Province, China)

葛新民，博士后，從事測井儲層評價、巖石物理實驗及解釋方面的研究；E-mail: gexinmin2002@163.com

(編輯 趙俊)