左文喆,王斌海,程紫華

(華北理工大學(xué)礦業(yè)工程學(xué)院，河北 唐山　063009)

?

復(fù)雜水文地質(zhì)條件下含水層識別在非穩(wěn)定流抽水試驗(yàn)中的應(yīng)用：以司家營鐵礦南區(qū)為例

左文喆,王斌海,程紫華

(華北理工大學(xué)礦業(yè)工程學(xué)院，河北 唐山063009)

司家營鐵礦南區(qū)水文地質(zhì)條件復(fù)雜，控水構(gòu)造發(fā)育，含水層類型具有明顯的分區(qū)特性，表現(xiàn)出了含水層類型的多樣性。為獲得司家營南區(qū)合適的水文地質(zhì)參數(shù)，提出一種新思路，同時(shí)考慮含水層類型和求參方法，利用含水層識別獲得含水層類型、補(bǔ)給條件、隔水邊界、井筒儲集等含水層性質(zhì)，再結(jié)合相應(yīng)的數(shù)學(xué)模型求參。研究表明，含水層識別在指導(dǎo)水文地質(zhì)求參方面表現(xiàn)出了它的優(yōu)越性，提高了求參精度，是一種指導(dǎo)正確求取水文地質(zhì)參數(shù)的手段。

司家營鐵礦南區(qū);水文地質(zhì)參數(shù);含水層識別；非穩(wěn)定流抽水實(shí)驗(yàn)

0　引　言

司家營鐵礦南區(qū)越流廣泛存在，構(gòu)造發(fā)育對礦體及其圍巖破壞較大，水文地質(zhì)條件復(fù)雜，求參難度較大。前期的水文地質(zhì)工程地質(zhì)工作主要圍繞第四系開展，基巖研究程度低，且在水文地質(zhì)求參上采用了單一的Theis配線法，由于礦區(qū)控水構(gòu)造發(fā)育，含水層類型分區(qū)明顯，求參采用單一的Theis配線法顯然是不合理的，利用該水文地質(zhì)參數(shù)會對地下水資源計(jì)算與評價(jià)、地下水?dāng)?shù)值模擬及預(yù)報(bào)造成誤差。本文提出一種求參新思路，解決含水層性質(zhì)不明條件下的水文地質(zhì)求參問題。

水文地質(zhì)參數(shù)的求解，經(jīng)歷了從穩(wěn)定流階段到非穩(wěn)定流階段，有Dupuit[1]和Thiem[2]提出的地下水流向水井的穩(wěn)定流解析解和Theis[3]提出的非穩(wěn)定流解析解，Theis公式標(biāo)志著多孔介質(zhì)滲流理論發(fā)展到非穩(wěn)定流的新階段[3-6]。除此外，還有利用Theis公式的Jacob直線圖解法和Theis水位恢復(fù)法為水文地質(zhì)求參的線性求解提供了方便[7]；Hantush公式能夠解決越流系統(tǒng)存在的狀況[8]。以上提出的井流方程旨在解決承壓含水層和存在越流條件下承壓含水層的求參問題，而潛水含水層與承壓含水層對井孔擾動(dòng)的響應(yīng)具有完全不同的特征，相應(yīng)的在解決此類求參上有考慮滯后的Boulton模型和考慮垂向流速的Neuman模型[9-10]。因此，根據(jù)井流模型可將含水層大體上劃分為Theis無越流補(bǔ)給承壓含水層、Hantush越流補(bǔ)給型承壓含水層、Boulton或Neuman潛水含水層[7]。

司家營鐵礦南區(qū)位于京山鐵路以南的山前傾斜平原區(qū)，地處灤河沖積扇一級階地中上部，區(qū)內(nèi)地形開闊平坦，地勢北高南低，坡降0.7‰～0.8‰，位于河北省灤縣和灤南縣交界。南區(qū)蓋層和結(jié)晶基底在沉積和構(gòu)造作用下共同構(gòu)成了賦水的孔隙和裂隙介質(zhì)，根據(jù)其孔隙的發(fā)育程度和裂隙的成因劃分出4個(gè)含水層，分別為第四系上部強(qiáng)含水層、第四系下部含水層、風(fēng)化帶裂隙含水層和構(gòu)造裂隙含水層。南區(qū)礦床上覆巨厚的第四系強(qiáng)含水體，使南區(qū)鐵礦成為大水呆滯礦；垂向上第四系下部含水層和黏土隔水層相互疊置、犬牙交錯(cuò)，造成黏土隔水層厚薄不一，局部形成獨(dú)立的含水層儲水系統(tǒng)，多區(qū)域黏土層缺失，越流廣泛存在；礦區(qū)北部發(fā)育F9、F10、F11、F12、F13多條斷層，東部新河斷裂帶形成一條近南北走向、平面北寬南窄、垂向上寬下窄似“花朵”狀的構(gòu)造裂隙含水帶，構(gòu)造裂隙的發(fā)育使不同含水層間在垂向和側(cè)向上產(chǎn)生水力聯(lián)系。天然條件下，該構(gòu)造裂隙帶主要發(fā)揮貯水功能，開采條件下轉(zhuǎn)化為礦床充水通道，第四系水通過底部黏性土層和強(qiáng)風(fēng)化帶垂向越流補(bǔ)給基巖弱風(fēng)化帶，然后沿構(gòu)造裂隙帶進(jìn)入礦坑。第四系水為礦床充水的主要水源，基巖構(gòu)造裂隙帶為礦體充水的主要通道，由此構(gòu)成了南區(qū)復(fù)雜的水文地質(zhì)條件[11]。

1　多孔抽水試驗(yàn)

司家營鐵礦南區(qū)完成了CK01、CK02、CK03、CK04共4個(gè)大孔的多孔抽水試驗(yàn)，各抽水大孔均布置在斷層及其影響帶上。多孔抽水試驗(yàn)為一個(gè)基巖抽水大孔附帶7～14個(gè)觀測孔的抽水試驗(yàn)(表1， 圖1)，抽水試驗(yàn)歷時(shí)72 h，抽水試驗(yàn)結(jié)束后開始停泵進(jìn)行水位恢復(fù)觀測，水位恢復(fù)歷時(shí)60 h。

表1抽水主孔與部分觀測孔徑向距離

Table 1The radial distance of the main holes and the partial observation holes

主孔觀測孔徑向距離/mCK01NK06255.12NK07428.06CGK01703.57CK02GK02 2.10CGK02223.54NK03415.89CK03GK03 4.25CGK03513.27NK12870.65CK04NK18345.15GK04 3.65CGK04260.84

目前，Theis公式求解已成為解決非穩(wěn)定流完整井流的主要求解方法，利用Theis公式計(jì)算含水層參數(shù)的傳統(tǒng)方法主要有配線法、Jacob直線圖解法和水位恢復(fù)法，求參方法上各有優(yōu)缺點(diǎn)，配線法可利用抽水資料的全部資料，避免個(gè)別數(shù)據(jù)的偶然誤差，其產(chǎn)生的誤差主要來源于3方面：抽水時(shí)間過短，早期擬合不理想；抽水后期曲線平緩易受人為操作影響；在配線上，無論早、中、后期，配線法隨意性較大，受人為干擾較大；Jacob直線解法避免了配線法的隨意性，但在徑向距離(r)和抽水時(shí)間(t)上約束較強(qiáng)，要求r較小，而t足夠大[12]；Theis水位恢復(fù)法是利用水位恢復(fù)期數(shù)據(jù)的線性解法，同樣也避免了配線法的隨意性和抽水的波動(dòng)，但恢復(fù)期數(shù)據(jù)處理較麻煩。

圖1　多孔抽水試驗(yàn)主孔及觀測孔布置情況Fig.1　The main and observation hole arrangement in porous pumping experiment

含水層的參數(shù)主要受含水層性質(zhì)和求參方法的影響，然而在推求水文地質(zhì)參數(shù)方面，國內(nèi)的研究重點(diǎn)放在改變求參方法和改進(jìn)數(shù)學(xué)模型上[6]，在含水層類型的認(rèn)知上主要依靠勘探資料，但在控水構(gòu)造發(fā)育的地區(qū)，地下水的賦存和運(yùn)移較復(fù)雜，單純依靠勘探資料是不全面的。因此，為了更精確地獲得水文地質(zhì)參數(shù)，本文提出一種新思路，在根據(jù)抽水試驗(yàn)數(shù)據(jù)求參時(shí)，先進(jìn)行含水層識別，再確定求參方法，然后利用含水層識別來獲得含水層類型、補(bǔ)給條件、隔水邊界、井筒貯集效應(yīng)等含水層性質(zhì)，最后結(jié)合相應(yīng)的數(shù)學(xué)模型求參。

2　含水層類型的識別

含水層的識別[13]可通過s-lgt(降深-時(shí)間)、ds/dt-lgt(降深速率-時(shí)間)曲線與標(biāo)準(zhǔn)曲線對比分析確定。各含水層類型的識別特征曲線描述如下。

(1)承壓含水層。在抽水早期，s-lgt曲線在單對數(shù)坐標(biāo)下表現(xiàn)為非線性段，但在抽水后期，表現(xiàn)為線性變化[12]。ds/dt-lgt曲線在早期表現(xiàn)為斜率漸增的曲線；在后期，水頭趨于等速下降，降深速率曲線平行于橫軸，這一階段為抽水后期的“等幅下降”階段。符合這一曲線特征的含水層類型定為Theis無越流補(bǔ)給承壓含水層，如NK06孔含水層識別曲線(圖2a)。

圖2　不同含水層類型識別曲線Fig.2　The identification curves of different aquifer types

(2)弱透水含水層。在抽水早期，s-lgt曲線與Theis曲線曲率變化一致，但在抽水后期，由于弱透水含水層的越流補(bǔ)給，s-lgt曲線漸漸偏離Theis曲線，并且平行于橫軸，達(dá)到抽水的穩(wěn)定階段。ds/dt-lgt曲線在早期也與承壓含水層具有相同變化特點(diǎn)，在補(bǔ)給段，變化率隨時(shí)間逐漸降低，當(dāng)抽水量達(dá)到補(bǔ)給量，ds/dt-lgt曲線與橫軸相交，變化率變?yōu)?。如GK02孔含水層識別曲線具有類似特點(diǎn)(圖2b)。

(3)隔水邊界的識別曲線特征。在抽水早期，隔水邊界與承壓含水層識別曲線特征仍一致，但在抽水后期由于隔水邊界的存在，水量補(bǔ)給不足，降深幅度增加，漸漸偏離Theis曲線，表現(xiàn)為斜率漸增的s-lgt曲線，而ds/dt-lgt曲線表現(xiàn)出較明顯斜率漸增的特征。如NK12孔含水層識別曲線具有類似特點(diǎn)(圖2c)。

(4)井筒儲集效應(yīng)是由于井筒在抽水早期，井筒中儲集的水量對抽水產(chǎn)生影響，使地下水降落產(chǎn)生滯后效應(yīng)[14]。井筒儲集效應(yīng)影響在抽水早期較大，且對ds/dt-lgt曲線的影響最為明顯，在單對數(shù)坐標(biāo)中主要表現(xiàn)為“單駝峰”式的特點(diǎn)(圖2d)，這種特點(diǎn)反映了井筒儲集從影響到消失的過程。

圖3　觀測孔GK02的Theis及Hantush模型擬合曲線Fig.3　Fitting curves of Theis and Hantush models of observation hole GK02

多孔抽水試驗(yàn)中，主孔的識別曲線均表現(xiàn)出井筒儲集效應(yīng)的現(xiàn)象(圖2d)。通過含水層識別，各水文孔含水層類型及邊界條件見表2。多孔抽水試驗(yàn)各觀測孔含水層識別均顯示，司家營南區(qū)含水層類型多樣，含水層類型有Theis無越流補(bǔ)給承壓含水層、Hantush越流補(bǔ)給型承壓含水層、雙重空隙介質(zhì)[15-17]的承壓含水層、具隔水邊界的承壓含水層和存在井筒儲集效應(yīng)的含水層(圖3)。

表2　含水層類型識別結(jié)果

3　實(shí)例分析

為說明含水層識別與求參的聯(lián)系，本文對部分水文孔的含水層識別進(jìn)行了詳細(xì)分析說明。

司家營鐵礦南區(qū)觀測孔GK02識別曲線(圖3)反映，在抽水后期s-lgt曲線表現(xiàn)出承壓含水層線性變化的特征，未表現(xiàn)補(bǔ)給邊界或越流承壓含水層“折平拐彎”的特征；而ds/dt-lgt識別曲線在抽水早期就表現(xiàn)出弱補(bǔ)給的特征，抽水期整體表現(xiàn)出具越流條件的“拱形式”特點(diǎn)。GK02孔附近雖受F9-F13復(fù)合斷裂的影響，但GK02孔附近黏土層較發(fā)育，含水層表現(xiàn)出越流滯后的現(xiàn)象，同時(shí)也反映出抽水時(shí)間不夠長的問題[18]。

GK02孔抽水前中期s-lgt識別曲線并不能明顯區(qū)別出補(bǔ)給的特性，s-lgt曲線識別特征的不明顯造成含水層類型的定性上偏向于Theis無越流的承壓含水層，而通過ds/dt-lgt曲線特征識別分析，含水層類型定為Hantush越流補(bǔ)給的承壓含水層更符合實(shí)際。

觀測孔GK04識別曲線(圖4a)反映，GK04孔附近含水層為具備孔隙-裂隙兩重性的雙重含水系統(tǒng)。這種含水系統(tǒng)中，孔隙主要發(fā)揮貯存地下水作用，且具有高貯存率、低傳導(dǎo)性的特點(diǎn)；裂隙以導(dǎo)水為主，具有低貯水率、高傳導(dǎo)性的特點(diǎn)，該含水系統(tǒng)在抽水過程中存在由于延遲彈性釋水的滯后現(xiàn)象。雙重介質(zhì)含水系統(tǒng)的s-lgt識別曲線在抽水早后期存在兩組互相平行的直線段，ds/dt-lgt表現(xiàn)為“雙拱式”的識別曲線。在GK04孔區(qū)域存在“天窗”，第四系下部含水層與基巖裂隙含水層存在直接的水力聯(lián)系，這與識別分析的含水層性質(zhì)一致。

根據(jù)GK02和GK04孔含水層識別結(jié)果，在參數(shù)擬合上分別采用了Theis(識別前)和Hantush、Double porosity(識別后)模型，參數(shù)識別結(jié)果見表3。參數(shù)表明含水層識別前后，水文地質(zhì)參數(shù)存在明顯的差異，而這種差異歸結(jié)于對含水層性質(zhì)認(rèn)識的不明確、選用不合理的求參方法導(dǎo)致的。

圖4　觀測孔GK04 Theis及Double porosity模型擬合曲線Fig.4　Fitting curves of Theis and Double porosity models of observation hole GK04

以上分析也暴露出傳統(tǒng)水文地質(zhì)求參的劣勢，傳統(tǒng)水文地質(zhì)求參往往會以任意模型的擬合精度去評價(jià)水文地質(zhì)參數(shù)的準(zhǔn)確性，含水層性質(zhì)認(rèn)識不明，造成模型的假定條件偏離實(shí)際的水文地質(zhì)條件，而結(jié)合含水層識別的水文地質(zhì)求參明確了含水層和數(shù)學(xué)模型的類型，表現(xiàn)出求參過程中的優(yōu)勢。

表3GK02、GK04孔不同求參方法的水文地質(zhì)參數(shù)

Table 3The hydrogeological parameters by different methods in hole GK02 and hole GK04

計(jì)算方法數(shù)學(xué)模型孔號T/(m2/d)K/(m/d)S識別前TheisGK0280.64.58×10-14.24×10-2TheisGK04605.03.672.09×10-4識別后HantushGK0286.44.90×10-18.72×10-2DoubleporosityGK04360.03.712.82×10-4

含水層識別過程中，觀測孔NK12孔在300～400 min表現(xiàn)出補(bǔ)給特性，但補(bǔ)給微弱，300 min之后降深驟增，s-lgt曲線偏離Theis曲線，ds/dt-lgt曲線斜率增大(圖2c)，這些特性反映了隔水邊界的存在，NK12孔附近含水層定為具隔水邊界的承壓含水層，在擬合上采用簡單的Theis擬合顯然與實(shí)際的水文地質(zhì)條件是不相符的。

根據(jù)含水層識別分析，NK12孔參數(shù)識別中應(yīng)運(yùn)用斯托爾曼法(Stallman)，為求參方便，采用了能夠分析特殊邊界條件的Aquifer test軟件[13]。針對NK12孔的含水層識別在300～400 min時(shí)段出現(xiàn)的補(bǔ)給段，將其分為兩個(gè)時(shí)段擬合，分別為0～300 min和400～4 000 min。擬合證明：在抽水段0～300 min和400～4 000 min，用具有隔水邊界的Theis模型擬合效果較好(圖5)，為驗(yàn)證擬合參數(shù)的正確性，將其與Theis水位恢復(fù)法擬合參數(shù)進(jìn)行對比(表4)。

擬合結(jié)果對比分析認(rèn)為，NK12孔0～300 min

圖5　NK12抽水期內(nèi)Theis擬合曲線Fig.5　Fitting curves of Theis of hole NK12 in pumping period

Table 4Hydrogeological parameters of hole using Theis formula

計(jì)算方法時(shí)段/minT/(m2/d)K/(m/d)S配線法 0～3002.38×1041.2×1022.70×10-4400～40001.72×1038.683.19×10-3水位恢復(fù)法4680～86407.85×1023.97-

水文地質(zhì)參數(shù)擬合結(jié)果偏大、不合理，其主要原因?yàn)椋涸?～300 min，抽水降深數(shù)據(jù)過少，雖然擬合效果較好，但主要擬合Theis曲線的后半段，造成井函數(shù)W(u)值偏大，降深偏小，導(dǎo)水系數(shù)和滲透系數(shù)偏大；選用在時(shí)段400～4 000 min的數(shù)據(jù)，用Theis配線法擬合結(jié)果更符合實(shí)際，引入邊界的Theis配線法擬合，根據(jù)含水層識別對數(shù)據(jù)進(jìn)行了篩分，從而提高了擬合的精度，突出了識別分析的重要性。

在水文地質(zhì)參數(shù)計(jì)算中，計(jì)算人員多是根據(jù)經(jīng)驗(yàn)進(jìn)行強(qiáng)制性劃分，選擇相應(yīng)公式計(jì)算，在含水層性質(zhì)不明的條件下用擬合精度來評價(jià)參數(shù)的準(zhǔn)確與否，顯然這種求參方法和結(jié)果是不合理的。因此，如何正確認(rèn)識含水層性質(zhì)，并選用合適的求參方法，對于求參顯得至關(guān)重要。利用抽水試驗(yàn)的s-lgt曲線、ds/dt-lgt曲線進(jìn)行含水層或抽水含水層影響因素識別，正是一種獲得含水層信息和指導(dǎo)使用正確求參方法的手段，并且含水層識別能夠降低求參方法造成的誤差，提高求參精度。

4　結(jié)　論

(1)司家營鐵礦南區(qū)控水構(gòu)造發(fā)育，含水層類型存在區(qū)域的分區(qū)特性，含水層類型有Theis無越流補(bǔ)給承壓含水層、Hantush越流補(bǔ)給型承壓含水層、雙重空隙介質(zhì)的承壓含水層、具隔水邊界的承壓含水層和存在井筒儲集效應(yīng)的含水層。

(2)含水層類型可以通過含水層識別確定，各類型含水層在含水層識別曲線上表現(xiàn)出各自明顯的特征，這些特征可識別含水層的類型、邊界條件、儲集效應(yīng)等含水層特征。

(3)含水層識別在指導(dǎo)水文地質(zhì)求參上表現(xiàn)出了它的優(yōu)越性。含水層識別降低了選用求參方法造成的誤差，提高了求參精度。

[1]DUPUIT J. Theoretical and Practical Studies of Water Flow in Open Channels and through Permeable Terrains[M].Second Edition. Paris: Dunod,1863:304.

[2]THIEM G.Hydrologische Methoden(Hydrological Methods)[M]. Leipzig: Gebhardt,1906:56.

[3]THEIS C V.The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground water storage[J]. Transactions of the American Geophysical Union, 1935,16(2):519-524.

[4]肖長來, 梁秀娟, 崔建銘, 等. 確定含水層參數(shù)的全程曲線擬合法[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2005(6):73-77.

[5]李佩成. 地下水非穩(wěn)定流解析法[M]. 北京：科學(xué)出版社，1990:20-30.

[6]李培月. 非穩(wěn)定流抽水試驗(yàn)確定越流承壓含水層水文地質(zhì)參數(shù)方法對比研究[D]. 西安：長安大學(xué), 2011:2-24.

[7]JACOB C E,LOHMAN O W.Nonsteady flow to a well of constant drawdown in an extensive aquifers[J].Transactions of the American Geophysical Union,1952,33(4):559-569.

[8]HANTUSH M S. Analysis of from pumping tests in leaky aquifers[J]. Transactions of the American Geophysical Union, 1956, 37(6):702-714.

[9]BOULTON N S. The drawdown of the water table under nonsteady conditions near a pumped well in an unconfined aquifer[J].Journal of Institution Civil Engineering,1954, 3(4):564-579.

[10]NEUMAN S P,WITHERSPOON P A.Theory of flowing a confined two aquifer system[J]. Water Resources, 1969,5(4):803-816.

[11]劉大金. 河北省灤南縣司家營鐵礦南區(qū)水文地質(zhì)工程地質(zhì)補(bǔ)充勘探報(bào)告[R]. 石家莊: 華北有色工程勘察院有限公司，2013:20-90.

[12]薛禹群. 地下水動(dòng)力學(xué)[M]. 北京：地質(zhì)出版社，1997:100-164.

[13]Hydrogeologic Incorporated.Aquifer Test v 4.2 User’s Manual[M]. Waterloo: Hydrogeologic Incorporated, 2007:121-175.

[14]劉日, 劉慈群. 考慮井筒存儲和表皮效應(yīng)的有限道路垂直裂縫井的試井分析方法[J]. 油氣井測試, 1993(2):1-10.

[15]羅二輝, 胡永樂, 王磊,等. 雙重介質(zhì)水平井產(chǎn)量遞減分析方法[J]. 現(xiàn)代地質(zhì),2013,27(6):1440-1444.

[16]王洪濤, 謝秀輝, 曹以臨. 北方巖溶含水層雙重介質(zhì)水流模型及其應(yīng)用探討[J]. 現(xiàn)代地質(zhì),1989，3(1):124-135.

[17]李愛芬, 劉照偉, 楊勇. 雙重介質(zhì)中有限導(dǎo)流垂直裂縫井試井模型求解的新方法[J]. 水動(dòng)力學(xué)研究與進(jìn)展, 2006,21(2):217-222.

[18]雷國輝, 葛國昌, 雷國剛, 等.承壓完整井非穩(wěn)定流抽水試驗(yàn)的觀測時(shí)間要求[J]. 水文地質(zhì)工程地質(zhì), 2009(1):40-43.

Application of Aquifer Identification in Complicated Geological Condition by Transient Flow Pumping Test：An Example of Southern Area of Sijiaying Iron Mine

ZUO Wenzhe,WANG Binhai,CHENG Zihua

(CollegeofMineEngineering,NorthChinaUniversityofScienceandTechnology，Tangshan,Hebei063009，China)

The southern area of Sijiaying iron mine is of complex hydrogeological conditions, and the water-controlling structures well develop. The aquifer type has an obvious zonal feature, showing the diversity of aquifer type, which undoubtedly increases the difficulty to infer the hydrogeological parameters of aquifer. Hydrogeolo-gical parameters of aquifer have become an important basis for groundwater numerical simulation to forecast and evaluate groundwater resources. Therefore, to verify the obtaining parameters under complex hydrogeological conditions in southern area of Sijiaying iron mine, the authors put forward a new idea, considering the aquifer type and method of obtaining parameters, using aquifer identification of the aquifer type to obtain the aquifer type, recharge condition, impervious boundary,and well storage.Combined with the corresponding mathema-tical model, the parameters are determined. This research shows that the aquifer identification has the superiority in the aspect of guiding to the calculation of hydrogeological parameters. This method increases the calculation accuracy and it is an effective way to obtain the correct parameters.

southern area of Sijiaying iron mine; hydrogeological parameter; aquifer identification; transient flow pumping test

2015-05-25；改回日期：2016-02-03；責(zé)任編輯：戚開靜。

河北省自然科學(xué)基金項(xiàng)目(D2011209071，E2013209328)；唐山市科技計(jì)劃項(xiàng)目(14130245B)。

左文喆，女， 博士，副教授， 1969年出生，水文地質(zhì)學(xué)專業(yè)，主要從事水文地質(zhì)及環(huán)境地質(zhì)研究。

Email: zuowenzhej@sina.com。

P641.4

A

1000-8527(2016)02-0478-06