尚建閣,李晨暉,蘇麗娟,李 冰,謝 崇,邢尚鑫

(1.河南省有色金屬礦產(chǎn)探測(cè)工程技術(shù)研究中心,河南 鄭州 450016;2.河南省有色金屬地質(zhì)礦產(chǎn)局第二地質(zhì)大隊(duì),河南 鄭州 450016;3.河南省自然資源科技創(chuàng)新中心(礦山生態(tài)環(huán)境保護(hù)修復(fù)研究),河南 鄭州 450016;4.河南有色地質(zhì)礦產(chǎn)集團(tuán)有限公司,河南 鄭州 450016;5.河南省有色金屬地質(zhì)勘查總院,河南 鄭州 450052;6.河南省有色金屬深部找礦勘查技術(shù)研究重點(diǎn)實(shí)驗(yàn)室,河南 鄭州 450052)

0 引言

隨著全球能源礦產(chǎn)的需求量不斷增加,找礦工作的程度也越來(lái)越深入。直流激發(fā)極化法是一種常規(guī)的物探方法,也是目前應(yīng)用最為廣泛的物探方法之一,通過(guò)該方法可以快速的獲取地球內(nèi)部介質(zhì)的極化率、電阻率等物理參數(shù),但地球表面并不是理論上的一個(gè)光滑的平面,而是一個(gè)地形高低不平、連綿起伏的復(fù)雜曲面,特別是在山區(qū),地形變化更加復(fù)雜,正是由于這些復(fù)雜地形的存在,在使用直流激發(fā)極化法開(kāi)展找礦過(guò)程中,發(fā)現(xiàn)的異常會(huì)產(chǎn)生畸變,可能導(dǎo)致淹沒(méi)地質(zhì)體的有用信息,呈現(xiàn)假異常,使研究者不能做出正確的推斷解釋,降低了物探綜合參數(shù)利用水平[1-3]。因此,復(fù)雜地形條件下電阻率地形改正方法的研究,成為物探工作者急需解決的關(guān)鍵問(wèn)題之一。ANSYS軟件是一款大型的融電磁學(xué)、結(jié)構(gòu)力學(xué)、流體力學(xué)、仿聲學(xué)等于一體的大型通用有限元分析軟件[4]。ANSYS具有豐富的單元庫(kù),提供了對(duì)各種物理場(chǎng)量的分析功能[5-6],為地球物理工作者在進(jìn)行數(shù)值模擬分析時(shí)提供了強(qiáng)大的數(shù)學(xué)分析功能。前人運(yùn)用ANSYS開(kāi)展了一系列的地球物理電法數(shù)據(jù)的數(shù)值模擬,對(duì)典型的地電模型全空間介質(zhì)影響進(jìn)行模擬,驗(yàn)證了ANSYS開(kāi)展電法數(shù)據(jù)模擬的正確性,均取得了較為理想的理論效果[7-9],但是對(duì)復(fù)雜地形條件下的模擬應(yīng)用較少。

文章通過(guò)介紹運(yùn)用ANSYS開(kāi)展電法地形改正數(shù)值模擬,總結(jié)經(jīng)驗(yàn)并運(yùn)用到實(shí)際生產(chǎn)找礦工作中,取得了良好應(yīng)用效果,為廣大地球物理工作者在復(fù)雜地形條件下開(kāi)展直流激發(fā)極化法數(shù)據(jù)解釋推斷提供技術(shù)參考。

1 方法原理

地形影響問(wèn)題是從20世紀(jì)50年代伴隨著直流電法在礦產(chǎn)資源勘查等方面的大規(guī)模應(yīng)用開(kāi)始受到重視的[10]。復(fù)雜的地形對(duì)直流激發(fā)極化法中電阻率參數(shù)的影響是致命的,因此,正確消除起伏地形對(duì)電阻率參數(shù)的影響,對(duì)后期地球物理成果解釋推斷的真實(shí)性起到了至關(guān)重要的作用[11]。

研究地形對(duì)視電阻率的影響,目前主要使用的有三種方法:物理模擬法、解析法和數(shù)值模擬法。所謂直流電法的地形改正,就是通過(guò)上述三種方法將在起伏地形上測(cè)量得到的電阻率數(shù)據(jù)轉(zhuǎn)換為平地形影響下數(shù)據(jù)的過(guò)程,從而消除不同地形所引起的誤差。由于物理模擬法和解析法不適用于場(chǎng)源和物性分布復(fù)雜的正演情況,因此,數(shù)值模擬法則成為電阻率法正演最主要的方法[12]。隨著計(jì)算機(jī)技術(shù)的發(fā)展和物探解釋方法的進(jìn)步,數(shù)值模擬技術(shù)得到了飛躍式發(fā)展。常用的電阻率數(shù)值模擬方法有有限差分法、積分方程法、邊界單元法和有限元法,其中有限元法具有較高的求解精度且理論較為成熟,能夠適應(yīng)任意模型形態(tài)。因此,有限元法是目前處理電法正演問(wèn)題的主流方法。

ANSYS是一款能實(shí)現(xiàn)多場(chǎng)及多場(chǎng)耦合分析、實(shí)現(xiàn)前后處理、求解及多場(chǎng)分析統(tǒng)一數(shù)據(jù)庫(kù)的一體化大型有限元分析軟件,一般而言,ANSYS的基本分析過(guò)程可以分為三步,即創(chuàng)建有限元模型、施加載荷進(jìn)行求解和查看分析結(jié)果。對(duì)應(yīng)軟件結(jié)構(gòu)的三個(gè)程序模塊:前處理模塊(PREP7)、分析求解模塊(SOLUTION)和后處理模塊(POST1和POST26)。前處理模塊(PREP7)的主要任務(wù)是建立結(jié)構(gòu)分析的有限單元模型,它是結(jié)構(gòu)分析的開(kāi)始,一般步驟為分析準(zhǔn)備、設(shè)置單元類(lèi)型、設(shè)定實(shí)常數(shù)、定義材料屬性、創(chuàng)建模型、劃分網(wǎng)格。分析求解模塊(SOLUTION)的主要任務(wù)是進(jìn)行結(jié)構(gòu)的計(jì)算,一般步驟:定義自由度、施加荷載、設(shè)定分析類(lèi)型、給定邊界條件及求解等。后處理模塊(POST1和POST26)的主要任務(wù)是進(jìn)行結(jié)構(gòu)的分析,它包括通用后處理器(POST1)和時(shí)程后處理器(PSOT26)兩種,前者主要用于靜態(tài)分析的結(jié)果后處理,而后者主要用于動(dòng)態(tài)分析結(jié)果的后處理。

2 模擬實(shí)驗(yàn)

利用ANSYS軟件進(jìn)行地形改正,其一般思路:① 根據(jù)電法工作的實(shí)際情況,依照測(cè)點(diǎn)高程建立物理模型,并要根據(jù)計(jì)算精度,對(duì)模型進(jìn)行網(wǎng)格剖分和邊界條件的設(shè)定;② 根據(jù)工作精度的要求,合理選擇單元類(lèi)型,利用ANSYS程序進(jìn)行純地形條件下各測(cè)點(diǎn)的視電阻率值計(jì)算;③ 通過(guò)比值法,對(duì)實(shí)測(cè)視電阻率數(shù)據(jù)進(jìn)行地形影響的壓制。計(jì)算公式采用ρs改=(ρs實(shí)/ρs地形)×ρ0,通過(guò)ANSYS有限元分析得到的純地形視電阻率數(shù)據(jù),對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行比值壓制,最終得到改正后的ρs改數(shù)據(jù)?？偟膩?lái)說(shuō),利用ANSYS軟件進(jìn)行地形改正需要正確記錄實(shí)際地形高程數(shù)據(jù),利用軟件所提供的有限元分析軟件計(jì)算出純地形條件下影響的數(shù)值,再利用實(shí)際測(cè)量獲得的電阻率數(shù)值與之進(jìn)行比值法進(jìn)行地形影響消除得到改正后的真實(shí)數(shù)值,基本處理流程見(jiàn)圖1。

圖1 ANSYS有限元數(shù)據(jù)分析流程圖

--------------------------------------------------------------------------------------------------------------

下面利用ANSYS編寫(xiě)命令流加以說(shuō)明。

(1)!進(jìn)入前處理器,建模

/prep7

/PNUM,area,1

et,1,plane230 ! 定義單元類(lèi)型

mp,1,resistivity ! 定義材料電阻率參數(shù)

*ask,KP_N,8,, ! 定義創(chuàng)建地形模型的點(diǎn)數(shù)

*dim,KP_NX,ARRAY,KP_N,,

*dim,KP_NY,ARRAY,KP_N,, ! 定義地形節(jié)點(diǎn)的XY數(shù)值

k,I,KP_NX(I),KP_NY(I) ! 定義關(guān)健節(jié)點(diǎn)

a,1,2,3,4,5,6,7,8…… ! 由點(diǎn)創(chuàng)建面

(2)進(jìn)入求解器!將模型網(wǎng)格劃分,定義邊界條件和施加荷載

asel,s,,,1,2

aatt,1

asel,s,,,4

aatt,2

asel,s,,,1,4

esize,,20 ! 定義網(wǎng)絡(luò)劃分尺度

mshape,1,2d ! 定義剖分類(lèi)型

amesh,all ! 開(kāi)始網(wǎng)絡(luò)劃分

lsscale,all

nsel,s,ext ! 選擇邊界上的節(jié)點(diǎn)

nsel,u,loc,y

*get,Nnod,node,,count

*get,Nmin,node,,num,min

*do,J,1,Nnod,1

Dist1=SQRT((NX(Nmin)-NAX)**2+(NY(Nmin)-NAY(Nmin))**2)

Dist2=SQRT((NX(Nmin)-NBX)**2+(NY(Nmin)-NBY(Nmin))**2)

V0=Currt*Res1*(1/Dist1-1/Dist2)/(2*Pi)

D,Nmin,VOLT,V0 ! 給邊界上的節(jié)點(diǎn)賦值

Nmin=ndnext(Nmin)

*enddo

fk,2,amps,1

fk,11,amps,-1 ! 施加電流荷載

finish

/sol

antype,static

! /status,solu ! 定義分析類(lèi)型和求解

solve

FINISH

(3)! 進(jìn)入后處理器,查看結(jié)果

/post

*SET,DJ,0.2

*SET,MN,0.2

*DIM,V1,ARRAY,20,1,

*DIM,Ra,TABLE,20,2,

*DO,I,1,20,1

V1(I,1)=VOLT(NODE(2+DJ*(I-1),6,0))-VOLT(NODE(2+DJ*I,6,0)) ! 求模擬電位差

AM=2+DJ*(I-1)

AN=2+DJ*I

BM=6-DJ*(I-1)

BN=6-DJ*I

KK=2*3.14156/(1/AM-1/AN-1/BM+1/BN)

Ra(I,1)=2+I*DJ

Ra(I,2)=V1(I,1)*kk/10

*enddo

FINISH

/GCOL,1,dis

/GCOL,2,Resistivity

*VPLOT,Ra(1,1),Ra(1,2),2

FINISH

根據(jù)ANSYS軟件進(jìn)行地形改正的思路和分析流程,分別開(kāi)展了不同坡度角的山谷地形和山脊地形對(duì)直流激發(fā)極化法電阻率影響的數(shù)值模擬(圖2)。圖2a為山脊模型,θ為坡度角,材料電阻率為100 Ω·m,為了研究山脊純地形對(duì)視電阻率的影響程度,分別進(jìn)行了坡度為0°、10°、20°、30°、40°的模型計(jì)算。由圖2a可見(jiàn),在山脊純地形上得到的視電阻率曲線呈現(xiàn)為低阻特征,并且隨著坡度角的增大,其低阻特征更加明顯。圖2b為山谷模型,θ為坡度角,材料電阻率為100 Ω·m,同樣也進(jìn)行了坡度為0°、10°、20°、30°、40°的模型計(jì)算。由圖2b可見(jiàn),在山谷純地形上得到的視電阻率曲線呈現(xiàn)為高阻特征,并且隨著坡度角的增大,其高阻特征更加明顯。由圖2可見(jiàn),純地形異常主要發(fā)生在角域頂點(diǎn)附近,也就是說(shuō)主要發(fā)生在坡度變化的地方。異常的大小與坡度角的大小有關(guān),即與地形起伏的陡緩程度關(guān)系密切。形態(tài)相同的角域,異常幅度隨坡度的增大而增大。

圖2 不同地形條件影響下模擬電阻率曲線圖

為了更好地研究地形對(duì)異常體電阻率曲線的影響,在數(shù)值模擬器中加入異常體進(jìn)行數(shù)值模擬,圖2c為山脊地形下高阻直立薄板電阻率曲線圖,θ為坡度角,背景電阻率為100 Ω·m,直立高阻體電阻率為1000 Ω·m。由圖2c可見(jiàn),在水平地面時(shí),直立薄板視電阻率異常表現(xiàn)為高阻特征,隨著地形坡度的變大,異常曲線表現(xiàn)為低阻特征,尤其當(dāng)坡度角為30°、40°時(shí),基本為純地形引起的異常曲線特征。圖2d為不同山谷地形下低阻水平薄板電阻率曲線,θ為坡度角,背景電阻率為100 Ω·m,水平薄板的電阻率為10 Ω·m。由圖2d可見(jiàn),在水平地面時(shí),水平薄板電阻率曲線表現(xiàn)為低阻特征,隨著地形坡度的變大,異常曲線表現(xiàn)為高阻特征,尤其當(dāng)坡度角為30°、40°時(shí),基本為純地形引起的異常曲線特征。由此可以看出,地形起伏對(duì)異常體引起的電阻率曲線影響很大,特別是當(dāng)?shù)匦吻懈顕?yán)重時(shí),對(duì)電阻率資料的影響更大,使異常體電阻率異常發(fā)現(xiàn)較大畸變,甚至出現(xiàn)相反的結(jié)果。圖2e為高阻直立薄板山脊地形改正前后的電阻率曲線。由圖2e可見(jiàn),藍(lán)色曲線為地改前的電阻率曲線,當(dāng)存在地形影響時(shí),高阻體也表現(xiàn)為低阻特征;紅色曲線為地改后的電阻率曲線,通過(guò)比值法改正后,電阻率曲線表現(xiàn)出高阻特征,與實(shí)際地下介質(zhì)情況相符。圖2f為高阻直立薄板山脊地形改正前后的電阻率曲線。由圖2f可見(jiàn),藍(lán)色曲線為地改前的電阻率曲線,當(dāng)存在地形影響時(shí),由于地形的影響,低阻體在電阻率曲線上也表現(xiàn)為高阻特征;紅色曲線為地改后的電阻率曲線,通過(guò)比值法改正后,電阻率曲線表現(xiàn)出低阻特征,與實(shí)際地?cái)嗝娼橘|(zhì)情況相符。通過(guò)電流源有限元分析計(jì)算,其結(jié)果充分表明,地形對(duì)電阻率的影響很大,特別是地形起伏較大地區(qū),影響程度更大。通過(guò)地形模擬實(shí)驗(yàn)可以總結(jié)得出地形對(duì)電阻率能產(chǎn)生較大影響,具有正地形展現(xiàn)出低阻特性,負(fù)地形展現(xiàn)出高阻特性,這樣的影響隨著地形的復(fù)雜程度不斷提高顯現(xiàn)的越發(fā)明顯,甚至能抵消地下目標(biāo)體所產(chǎn)生的真實(shí)電阻率影響,使得目標(biāo)體的真實(shí)特征不能被探獲。利用ANSYS軟件進(jìn)行快速地形改正有效地解決了地形起伏對(duì)電阻率的影響,使反映的介質(zhì)電阻率回歸為真實(shí)。

在使用ANSYS軟件進(jìn)行數(shù)值分析的過(guò)程中,網(wǎng)格剖分的大小直接影響著計(jì)算結(jié)果的精度,所以網(wǎng)格剖分越細(xì),計(jì)算精度也就有所提高。但是,網(wǎng)格劃分的大小也影響著計(jì)算工作量的大小,直接影響著計(jì)算速度。因此,在應(yīng)用中應(yīng)權(quán)衡這兩種因素,做到綜合考慮。圖3為不同網(wǎng)格剖分密度下的山脊地形電阻率曲線。由圖3可見(jiàn),隨著剖分尺度的大小變化,曲線的光滑程度也發(fā)生變化,且剖分尺度越大,曲線振蕩程度越小,說(shuō)明計(jì)算精度越高。在ANSYS穩(wěn)態(tài)電流有限元分析中,程序沒(méi)有提供遠(yuǎn)邊界單元類(lèi)型,在實(shí)際應(yīng)用時(shí),須構(gòu)建一個(gè)遠(yuǎn)邊界模型,對(duì)研究區(qū)和邊界區(qū)進(jìn)行不同的網(wǎng)格劃分。

圖3 二維純山脊地形電阻率曲線圖(不同剖分尺度)

3 應(yīng)用實(shí)例

河南省西部山區(qū)是我國(guó)重要的有色金屬成礦區(qū)帶之一,該區(qū)域分布有多個(gè)大大小小的有色貴金屬礦床,為河南省礦產(chǎn)資源需求提供了保障[13-14]。研究區(qū)位于河南省西部洛寧縣,熊耳山北坡,其大地構(gòu)造位置處于中朝準(zhǔn)地臺(tái)南緣、華熊臺(tái)緣凹陷、崤山—魯山拱褶斷束中部(圖4)。地層出露主要為太古宇太華群的角閃斜長(zhǎng)片麻巖、熊耳群的安山巖和第四系。區(qū)內(nèi)巖漿活動(dòng)頻繁,巖漿巖種類(lèi)較多,其中以燕山期酸性巖最為發(fā)育,并與本區(qū)內(nèi)生金屬礦產(chǎn)的形成關(guān)系極為密切。斷裂構(gòu)造十分發(fā)育,并以NE—NNE向斷裂為主,且規(guī)模較大。研究區(qū)圍巖蝕變主要為硅化、絹云母化、綠泥石化、黃鐵礦化、泥化等。硅化表現(xiàn)為石英沿破碎帶充填,形成巨厚層的硅化碎裂巖帶出現(xiàn)。構(gòu)造蝕變帶內(nèi)礦化與蝕變主要有黃鐵礦化、方鉛礦化、碳酸鹽化、高嶺土化、褐鐵礦化、硅化、絹云母化等,金屬礦物表現(xiàn)形式為星點(diǎn)狀、浸染狀、細(xì)脈狀、細(xì)脈侵染狀、團(tuán)塊狀。

圖4 區(qū)域地質(zhì)礦產(chǎn)圖(a)與研究區(qū)地質(zhì)簡(jiǎn)圖(b)

表1為研究區(qū)電性參數(shù)統(tǒng)計(jì)結(jié)果,由表1可見(jiàn),研究區(qū)礦石(含礦化體)與圍巖存在著較大的物性差異,目標(biāo)體電性表現(xiàn)特征為低阻高極化特征,為在該研究區(qū)開(kāi)展物探工作提供物性前提。

表1 研究區(qū)電性參數(shù)統(tǒng)計(jì)結(jié)果

本次物探研究工作是在試驗(yàn)區(qū)K9號(hào)礦脈上進(jìn)行,采用時(shí)間域激電偶極-偶極測(cè)深裝置,點(diǎn)距20 m,隔離系數(shù)為5,MN極距40 m。儀器為V8,最大供電電流為14 A,6個(gè)電通道方式。圖5為試驗(yàn)剖面測(cè)量結(jié)果,由圖5a可見(jiàn),極化率斷面圖,極化率高值異常與K9號(hào)礦脈基本對(duì)應(yīng),反映出地形對(duì)極化率測(cè)量結(jié)果影響較小,而圖5b中所反映電阻率異常與礦脈不對(duì)應(yīng),并且展現(xiàn)出了高電阻率的特征,這一高電阻率特征是地下礦脈真實(shí)電阻率與地形影響產(chǎn)生的電阻率疊加影響產(chǎn)生的綜合反映,其低阻異常向左側(cè)地形相對(duì)低緩的一側(cè)進(jìn)行偏離。這是由于地形切割較大,山谷地形純電阻率異常反映為高阻,致使礦體引起的低電阻異常被掩蓋,發(fā)生偏離現(xiàn)象。圖5c為采用ANSYS軟件計(jì)算出該地形所產(chǎn)生的響應(yīng)值與實(shí)際探測(cè)的電阻率值的比值進(jìn)行地形改正計(jì)算,其結(jié)果可見(jiàn)低電阻異常位置向右側(cè)偏移,并與K9號(hào)礦脈基本對(duì)應(yīng),徹底消除了復(fù)雜地形對(duì)探測(cè)結(jié)果的影響效應(yīng)。結(jié)合圖5a和圖5c結(jié)果顯示與巖礦石物性結(jié)果一致,表現(xiàn)為低阻高極率特征,證明采用ANSYS軟件進(jìn)行電阻率地形改正使測(cè)量結(jié)果回歸真實(shí)反映,效果顯著。

4 結(jié)論

(1)激發(fā)極化法是地質(zhì)找礦中應(yīng)用最為廣泛的地球物理方法之一,但是其電阻率數(shù)據(jù)收復(fù)雜地形影響也最為嚴(yán)重,負(fù)地形將形成高電阻率反映,而正地形將產(chǎn)生低電阻率反映,從而實(shí)際物探測(cè)量所獲得電阻率帶有地形加載的影響,使得大地介質(zhì)的真實(shí)電阻率得不到有效體現(xiàn)。

(2)ANSYS軟件是一款能實(shí)現(xiàn)多場(chǎng)及多場(chǎng)耦合分析、實(shí)現(xiàn)前后處理、求解及多場(chǎng)分析統(tǒng)一數(shù)據(jù)庫(kù)的一體化大型有限元分析軟件,能快速準(zhǔn)確進(jìn)行多種場(chǎng)景數(shù)值模擬計(jì)算,能夠?qū)?fù)雜地形條件下不同裝置類(lèi)型的電阻率開(kāi)展二維、三維地形改正,利用該軟件有限元分析功能模擬計(jì)算出純地形下產(chǎn)生的電阻率響應(yīng),通過(guò)比值法將野外所獲取的電阻率值中地形影響消除,獲得大地介質(zhì)真實(shí)電阻率,實(shí)現(xiàn)復(fù)雜地形對(duì)電阻率值的零影響效應(yīng),為物探解釋的可靠性提供保障。

(3)應(yīng)用ANSYS軟件在研究區(qū)復(fù)雜地形條件下開(kāi)展電阻率地形改正實(shí)際應(yīng)用,對(duì)比結(jié)果顯示,運(yùn)用該軟件進(jìn)行地形改正方法可行,效果顯著。