范翠松，李桐林，嚴加永

1 吉林大學地球探測科學與技術學院，長春 130026

2 中國地質科學院礦產(chǎn)資源研究所，北京 100037

1 引 言

為了有效地消除電磁耦合和評價激電異常，Pelton等［1］提出了頻譜激電法（SIP）.它利用常規(guī)電阻率法的觀測方式，在若干個頻率（頻率范圍n×0.01～n×102Hz）下測量視復電阻率譜和視相位譜，推斷目標體的真復頻譜參數(shù)特性，從而提高對異?？陀^評價的準確性.Pelton等［1］認為，對均勻巖、礦石，由激電效應引起的復電阻率隨頻率的變化可由Cole－Cole模型描述，并通過對大量巖、礦石標本的廣域頻譜復電阻率測量，驗證了該模型的正確性，并指出可以根據(jù)巖、礦石的Cole－Cole模型參數(shù)，按結構評價極化異常體.

由于SIP法可觀測的參數(shù)多，多參數(shù)組合的解釋能為評價激電異常源的性質提供更豐富的信息，且勘探成本低，周期短，抗干擾能力強，近年已在礦產(chǎn)勘查、油氣藏探測及環(huán)境調查等方面得到了廣泛應用［2-4］.尤其在隱伏多金屬礦勘探方面，復電阻率法能較好地解決礦與非礦異常的區(qū)分問題，提高了尋找隱伏礦的能力.

隨著SIP法的理論研究、儀器設備及施工技術日益成熟，與之相應的反演算法也得到了快速發(fā)展.我國學者羅延鐘提出了由視譜參數(shù)解析計算極化體真譜參數(shù)的方法［5］.張桂青等［6］提出由視譜直接反演極化體真譜的思想，并建立了稀釋系數(shù)為常數(shù)及隨頻率變化下的反演方法.劉崧等［7］提出了聯(lián)合譜激電反演真Cole－Cole參數(shù)的思想，并實現(xiàn)了橢球極化體的真譜參數(shù)反演.蔡軍濤等［8］在直流電阻率法數(shù)值模擬的基礎上，提出了利用最優(yōu)化反演和遞推反演相結合求取真頻參數(shù)的反演方法［9］.Loke等［10］同樣在不考慮電磁感應的情況下，以近似反演得到的非均勻半空間作為初始模型，實現(xiàn)了理論數(shù)據(jù)的SIP二維反演.就以上SIP反演方法而言，在面向實際應用時存在著一些難以避免的問題.利用視譜直接擬合真譜的方法，在建立視譜和真譜的近似關系式時，難以對普遍的情況使用統(tǒng)一的復電阻率模型進行描述，在解釋復雜構造異常體及真譜參數(shù)的空間分布上存在著一定的局限性，并且將電磁效應僅用一個基于Cole－Cole模型的簡單數(shù)學表達式描述時，也很可能造成解釋上的偏差.而以直流電阻率法正演為基礎的SIP反演算法，一般均建立在可忽略電磁耦合這一假設之上，而實測的SIP數(shù)據(jù)中卻一定包含了此種效應，特別是在觀測頻率較高時，電磁效應更為顯著，這意味著電磁耦合“噪聲”在此類反演中將起著更大的干擾作用.

近年來，由頻率域Maxwell方程出發(fā)建立的SIP反演方法得到了廣泛研究.張輝等［11］基于體積分方程的SIP正演模擬，首先提出并完成了基于阻尼最小二乘法的SIP三維反演，李建平［12］在此基礎上進行了改進，實現(xiàn)了起伏地形下的SIP三維反演.Ahmad Ghorbani［13］編寫了包含電磁效應的SIP一維反演代碼，并驗證了其正確性.徐凱軍［14］將解耦的波數(shù)域電磁場偏微分方程用有限元法實現(xiàn)了正演模擬，并利用阻尼最小二乘法通過對電場的擬合，研究了復雜地形的2.5維SIP局部反演.梁盛軍［15］將共軛梯度法和阻尼最小二乘法按先后順序應用于三維SIP反演中，并對理論模型進行了局部反演計算.可以看出，雖然以上反演研究既考慮了電磁效應和激電效應并存的情況，又能直接反演出目標體的多種復電阻率參數(shù)、位置和幾何特征.但是，由于反演單元存在多個Cole－Cole模型參數(shù)，在二維和三維全區(qū)反演時，模型參數(shù)量遠多于觀測數(shù)據(jù)量，導致反演問題的欠定性嚴重.因此，該方面的研究目前也僅停滯在理論研究及小規(guī)模的模型試算階段，尚未達到實際應用的程度.

針對頻譜激電反演中存在的問題，尤其是目前最符合實際應用的2.5維反演問題，本文提出并實現(xiàn)了利用多個排列的視電阻率和視相位數(shù)據(jù)的聯(lián)合反演方法，在考慮電磁效應的情況下，能夠同時反演出二維地質斷面上所有單元的復電阻率參數(shù)，解決了SIP反演欠定性嚴重的問題.算法中利用最小二乘原理構建了反演目標函數(shù)，并加入了Occam反演［16-17］的光滑模型約束，以增加反演的穩(wěn)定性.對于反演中的雅克比矩陣的計算問題，本文借助了電場的偏導數(shù)形式，推導出了靈敏度元素的解析表達式，并應用互換定理［18-20］對其進行直接求解.最后，對安徽某地區(qū)的SIP實測數(shù)據(jù)進行了2.5維反演成像，通過與已知鉆孔資料以及和CSAMT反演結果的對比分析，驗證了該反演算法的應用效果.

2 2.5維復電阻率反演原理

SIP反演問題可歸結為求解以下泛函極值的問題：

其中，x為模型的復電阻率參數(shù)向量；ε為數(shù)據(jù)擬合誤差；f為實測數(shù)據(jù)向量；F（x）為正演函數(shù).考慮到電磁效應對觀測數(shù)據(jù)的影響，正演模擬是以波數(shù)-頻率域的電磁耦合方程［21］為理論基礎，將實電阻率變量用（2）式的Cole－Cole模型替換后，再用等參單元的有限元法［22］完成了離散化處理，最終可得到由式（3）表示的正演線性方程組［23］：

式中，ρ（iω）為復電阻率；ω 是角頻率；ρ0、m、τ、c為Cole－Cole模型參數(shù)，分別表示零頻電阻率、極化率、時間常數(shù)、頻率相關系數(shù).

式中，K是對稱、帶狀的系數(shù)矩陣；F為剖分網(wǎng)格節(jié)點處的波數(shù)域待求場向量；B為發(fā)射源項.當求解出空間域電場值后，利用以下近似公式就能計算出測點處視電阻率及視相位值.

式中，ρa和φa分別為觀測點處的視電阻率值與視相位值；K是裝置系數(shù)；E為觀測偶極中心處的電場值；ΔL為觀測偶極長度；I為發(fā)射電流.

2.1 反演方程組的建立

為了充分利用視電阻率和視相位數(shù)據(jù)包含的地質信息，本文采用了聯(lián)合反演的策略.從公式（4）中不難發(fā)現(xiàn)，復電阻率法的反演擬合，本質上是基于電場的振幅與相位的擬合，因此，文中給出的反演算法均針對電場的振幅及相位進行推導.考慮到四種SIP參數(shù)的量級不同，為了提高反演的穩(wěn)定性，對參數(shù)進行了對數(shù)歸一化處理，即令x ＝ （ln（ρ0），ln（m），ln（τ），ln（c））.遵循以上規(guī)則，以 Occam 原理構建了單一排列下的反演目標函數(shù)：

式中，a和φ分別為觀測電場的振幅向量與相位向量；A（x）和φ（x）分別為電場振幅與相位的正演函數(shù)；Wa和Wφ為歸一化對角陣；λ是阻尼因子；μ為縮放系數(shù)；R為模型的二階粗糙度矩陣.

如圖4，由∠BPA′=30°，∠ABO=60°，得∠BEP=90°，∠A′EO=90°.所以所以由得于是點P坐標為

將（5）式在k次迭代模型xk鄰域線性展開，并極小化，可得到以下實系數(shù)反演方程組：

其中，ΔdA和Δdφ分別為電場振幅和相位的絕對擬合差向量；JA和Jφ分別為電場振幅與相位的靈敏度矩陣；Δxk＝xk＋1－xk為模型修正向量.通過求解以上方程組，就能得到本次迭代更新的模型參數(shù)向量xk＋1，再將其作為初始模型進行下一次迭代計算，直到擬合誤差收斂至設定的閾值.

當利用M個排列的數(shù)據(jù)反演時，首先求出各排列對應于（6）式中左端的靈敏度相關項P1，P2…，PM，及右端項S1，S2…，SM，再令其線性相加生成新的靈敏度相關項Q＝P1＋P2＋…＋PM及右端項T＝S1＋S2＋…＋SM，此時方程組（6）變?yōu)椋?/p>

通過上式便可實現(xiàn)利用多個排列數(shù)據(jù)的SIP反演計算.

2.2 靈敏度矩陣推導

靈敏度矩陣的建立是反演的重要環(huán)節(jié)，它直接決定了反演成像的時間和精度.在（6）式中可以看出，反演方程中包含了JA和Jφ兩種類型的靈敏度矩陣，在反演過程中需要分別求取，其解析表達可借助場的靈敏度形式獲得.在頻率域中，將空間分布的電場表示成以下復數(shù)形式：

再根據(jù)復合函數(shù)求導法則，便可推導出以下各復電阻率參數(shù)的靈敏度解析表達式：

在（10）式中，exj ＝exp（xj）＝ （ρ0，m，τ，c）；ω 為角頻率；ρ為復電阻率變量，與σ互為倒數(shù)為j單元的復電阻率對其Cole－Cole模型參數(shù)的偏導數(shù)，可解析求?。?4］.這樣，靈敏度元素JAj和Jφj的計算，最終便歸結到電場對第j塊單元復電阻率的偏導數(shù)計算問題，即求解，為節(jié)約計算成本，應用互換定理，并通過偽正演對該項進行計算.鑒于互換定理在2.5維反演中應用的復雜性，限于篇幅，擬另文發(fā)表.

3 安徽某斑巖銅礦上的應用試驗

3.1 礦床地質

安徽某斑巖銅礦床位于廬樅火山巖盆地西北緣，礦床的形成與中國東部燕山期巖漿侵入與噴發(fā)活動有關，是東部地區(qū)一個典型的巖漿熱液型礦床.礦區(qū)地層簡單，主要有第四系（Q）、白堊系楊灣組（K1y）、志留系高家邊組（S1g）和墳頭組（S2f），早侏羅統(tǒng)磨山組（J1m）和中侏羅統(tǒng)羅嶺組（J2l）地層以及早白堊紀龍門院組和浮山組火山巖系.

礦區(qū)構造主要表現(xiàn)為褶皺和斷裂兩種類型.高家邊組和墳頭組地層組成NNE向背斜，背斜上斷裂交匯處為銅礦的有利富集部位.區(qū)內巖漿活動強烈，形成了一套以石英閃長斑巖和黑云母石英閃長斑巖為主的鈣堿性系列的中酸性巖體，總體呈北東向分布.侵入于志留系和侏羅系地層中的石英閃長斑巖、黑云母石英閃長斑巖等為主要容礦巖體，總體呈北東—北北東向沿背斜核部分布［25］.

礦石礦物成分簡單，主要金屬礦物有黃鐵礦、黃銅礦、斑銅礦、輝鉬礦、磁鐵礦、輝銅礦；非金屬礦物除原巖蝕變礦物及交代殘余礦物外，還有石英和長石等.礦石構造以浸染狀、細脈狀和細脈浸染狀為主，其中含銅斑巖型礦石普遍為浸染狀礦化疊加疏密不等的細脈狀礦化.

3.2 野外數(shù)據(jù)采集

數(shù)據(jù)采集時選擇了一條有鉆井控制的剖面，使用的測量儀器是加拿大Phoenix地球物理公司研制生產(chǎn)的V8多功能電法測量系統(tǒng)，觀測方式為偶極-偶極.共進行了16組排列的數(shù)據(jù)采集，每排列12個觀測道，發(fā)射偶極長度為200m，接收偶極長度為50m，共進行25個頻率（0.0313～128Hz）的觀測.反演前，對觀測數(shù)據(jù)進行了篩選，剔除了部分干擾較強的排列、頻點及觀測道.

3.3 反演結果與分析

利用本文提出并實現(xiàn)的2.5維SIP反演算法，對該地區(qū)實測數(shù)據(jù)進行了反演成像.反演擬合情況見圖1～4.其中，圖1和圖2分別為1Hz的視電阻率和視相位擬合斷面圖；圖3和圖4分別為16Hz的視電阻率和視相位擬合斷面圖.經(jīng)對比可以發(fā)現(xiàn)，兩個頻率數(shù)據(jù)的反演擬合情況大致相當，總體擬合程度較好，僅局部存在微小的差異.但視相位較視電阻率的擬合程度稍差，這是因為視相位對SIP參數(shù)的靈敏度較視電阻率高，同時也說明了視相位數(shù)據(jù)提供了更多的SIP異常信息，在反演中起到了更為關鍵的作用.

四個SIP參數(shù)的反演結果見圖5，圖中白線圍成的封閉區(qū)域為鉆井控制的礦體范圍.圖6是該剖面的CSAMT法反演結果.由反演結果不難看出，SIP反演的零頻電阻率與CSAMT法反演的電阻率剖面吻合較好，其電阻率結構可以反映出礦區(qū)的地層、巖體及背斜構造.主體高阻區(qū)是含礦巖體的反映，而分支的高阻體應該是巖支的表現(xiàn)，而低阻區(qū)域在主體上應該是志留系地層的表現(xiàn).對于礦體而言，主體礦體具有較高的電阻率，這與礦床的類型和含礦體的電阻率測量結果是一致的（表1）.斑巖型礦體一般具有高的電阻率，特別是石英斑巖具有較其它斑巖更大的電阻率值.SIP反演結果顯示，礦體具有較小的時間常數(shù)，這與侵染狀或細脈狀礦石的時間常數(shù)較小的特征較為一致.在礦體范圍，頻率相關系數(shù)具有由淺至深逐漸變大的趨勢，這也同該參數(shù)的物理意義一致，因為該參數(shù)主要反映礦物顆粒的均一性特征.礦體的淺部由于埋藏淺，不均勻性嚴重，因此頻率相關系數(shù)較小，而深部壓力大，形成的顆粒更均質，頻率相關系數(shù)會相對較大.

圖2 1Hz視相位擬合斷面圖（a）觀測數(shù)據(jù)；（b）正演數(shù)據(jù).Fig.2 Fitting section of apparent phase with 1Hz（a）Measured data；（b）Forward data.

圖3 16Hz視電阻率擬合斷面圖（a）觀測數(shù)據(jù)；（b）正演數(shù)據(jù).Fig.3 Fitting section of apparent resistivity with 16Hz（a）Measured data；（b）Forward data.

圖4 16Hz視相位擬合斷面圖（a）觀測數(shù)據(jù)；（b）正演數(shù)據(jù).Fig.4 Fitting section of apparent phase with 16Hz（a）Measured data；（b）Forward data.

圖5 SIP反演斷面圖（a）零頻電阻率；（b）極化率；（c）頻率相關系數(shù)；（d）時間常數(shù).圖中由白色實線圍成的封閉區(qū)域為鉆孔控制的礦體范圍，下同.Fig.5 Cross－section of SIP inversion（a）Result of zero frequency resistivity；（b）Result of polarization coefficient；（c）Result of frequency correlation coefficient；（d）Result of time constant.The closed area surrounded by the solid white line is the range of orebody，which is controlled by drilling.It＇s the same in the following figure.

圖6 CSAMT反演斷面圖Fig.6 Cross－section of CSAMT inversion

表1 礦區(qū)巖石物性參數(shù)表Table 1 Rock physical property parameters in study area

反演的極化率在礦體上約為20.0%，與標本測量結果基本一致，但并沒有出現(xiàn)人們期待的高極化異常.相反，志留系地層則表現(xiàn)出較高的極化率特征，這同粘土質粉砂巖標本相對較低的極化率（7.0%）相矛盾.我們認為，可能是粘土質粉砂巖標本的代表性出現(xiàn)了問題，志留系砂泥巖可能存在較多的碳質成分，在后期巖體侵入過程中接觸變質，因此應該表現(xiàn)為較成礦巖體具有更高的極化率特征.

就該區(qū)的斑巖型銅礦而言，不難看出：該礦體具有高電阻率，相對低極化，中等頻率相關系數(shù)和較小的時間常數(shù)的特點.基于如上認識，推測反演圖右側（水平坐標2600，深度約200m）的異常可能是礦致異常，同樣，左側（水平坐標1200，深度約200m）也可能是相同類型的礦致異常.

對于剖面上其它參數(shù)組合的異常，我們也進行了分析和推斷.在水平坐標1500m，深度約200m處存在的異常，表現(xiàn)為電阻率低，極化率較低，頻率相關系數(shù)高和時間常數(shù)高的特點，推斷為隱伏角礫巖.淺地表的相對高阻（相對于圍巖地層），高極化率，低頻率相關系數(shù)和低時間常數(shù)的異常，推斷為斑巖型銅礦的外圍似千枚巖化蝕變帶（黃鐵礦化，石英和絹云母化）.

4 結 論

針對復電阻率法反演在實際應用中存在的問題，特別是生產(chǎn)中急需的2.5維反演問題，本文提出并實現(xiàn)了2.5維復電阻率法反演，并進行了應用試驗，得到了如下結論：

（1）本文提出并實現(xiàn)的復電阻率反演算法，不對電磁效應做任何假定，是建立在頻率域Maxwell方程基礎上的完全意義上的擬合反演.不但能同時反演出全區(qū)的零頻電阻率、極化率、頻率相關系數(shù)和時間常數(shù)，而且能反演出異常體的幾何參數(shù).

（2）算法中同時利用了視電阻率和視相位數(shù)據(jù)進行聯(lián)合反演擬合，在增加數(shù)據(jù)信息的基礎上，有效的提高了反演的分辨率.通過擬合斷面圖的對比發(fā)現(xiàn)，視相位對模型參數(shù)的靈敏性高于視電阻率，同時也說明了視相位數(shù)據(jù)提供了更多的SIP異常信息，在反演中起到了更為關鍵的作用.

（3）利用反演程序對已知礦區(qū)的實測數(shù)據(jù)進行了SIP反演實驗，通過與CSAMT的反演結果、鉆探資料和物性資料的對比驗證，表明了該反演方法能成功地圈定礦區(qū)的地質構造及礦體，使按結構區(qū)分礦與非礦成為可能，具有良好的應用效果.

（4）從反演實例可以看出，即使SIP的多參數(shù)反演能夠提供更為豐富的異常信息，但在實際生產(chǎn)中，地質情況往往更加復雜，僅依靠反演結果難以實現(xiàn)精確的地質解釋.因此，需緊密結合實際地質資料進行綜合分析，才能更準確地解決對礦床的定位與預測問題.

（References）

［1］Pelton W H，Ward S H，Hallof P G，et al.Mineral discrimination and removal of inductive coupling with multifrequency IP.Geophysics，1978，43（3）：588－609.

［2］Vanhala H，Soininen H，Kukkonen I.Detecting organic chemical contaminants by spectral－induced polarization method in glacial till environment.Geophysics，1992，57（8）：1014－1017.

［3］Weller A，Borner F D.Measurements of spectral induced polarization for environmental purposes.Environmental Geology，1996，27（4）：329－334.

［4］Vanhala H.Mapping oil－contaminated sand and till with the spectral induced polarization （SIP）method.Geophysical Prospecting，1997，45（2）：303－326.

［5］羅延鐘，方勝.視復電阻率頻譜的一種近似反演方法.地球科學-武漢地質學院學報，1986，11（1）：93－102.Luo T Z，F(xiàn)ang S.An approximate inversion of the apparent complex resistivity spectrum. Earth Science－Journal of Wuhan College of Geology （in Chinese），1986，11（1）：93－102.

［6］張桂青，崔先文，羅延鐘.一種反演頻譜激電法視頻譜求取真參數(shù)的方法.地質與勘探，1987，23（4）：48－54.Zhang G Q，Cui X W，Lou T Z.The determination of intrinsic parameters by inversing IP apparent spectrum.Geology and Prospecting （in Chinese），1987，23（4）：48－54.

［7］劉崧，官善友，高鵬飛.求極化橢球體真Cole－Cole參數(shù)的聯(lián)合譜激電反演.地球物理學報，1994，37（增刊2）：542－551.Liu S，Guan S Y，Gao P F.Joint SIP inversion for estimation of intrinsic Cole－Cole parameters of a polarizable ellipsoid.Chinese J.Geophys.（Acta Geophysica Sinica）（in Chinese），1994，37（SuppⅡ）：542－551.

［8］蔡軍濤，阮百堯，趙國澤等.復電阻率法二維有限元數(shù)值模擬.地球物理學報，2007，50（6）：1869－1876.Cai J T，Ruan B R，Zhao G Z，et al.Two－dimensional modeling of complex resistivity using finite element method.Chinese J.Geophys.（in Chinese），2007，50（6）：1869－1876.

［9］蔡軍濤，阮百堯，羅潤林.一種快速準確求取復電阻率真頻參數(shù)的反演方法.工程地球物理學報，2005，2（5）：338－342.Cai J T，Ruan B R，Luo R L.A rapid and precise inversion method for the determination of intrinsic complex resistivity parameters.Chinese J.Engineering Geophysics （in Chinese），2005，2（5）：338－342.

［10］Loke M H，Chambers J E，Ogilvy R D.Inversion of 2D spectral induced polarization imaging data.Geophysical Prospecting，2006，54（3）：287－301.

［11］張輝.復電阻率三維電磁場正反演研究［博士論文］.長春：吉林大學地球探測科學與技術學院，2006.Zhang H.Research of complex resistivity 3Delectromagnetic forward and inversion ［Ph.D.thesis］（in Chinese）.Changchun：College of Geo－Exploration Science and Technology，Jilin University，2006.

［12］李建平.帶地形的三維復電阻率電磁場正反演研究［博士論文］.長春：吉林大學地球探測科學與技術學院，2008.Li J P.Research on electromagnetic modeling of 3Dcomplex resistivity with topography［Ph.D.thesis］（in Chinese）.Changchun：College of Geo－exploration Science and Technology，Jilin University.2008.

［13］Ghorbani A，Camerlynck C，F(xiàn)lorsch N.CR1Dinv：A Matlab program to invert 1Dspectral induced polarization data for the Cole－Cole model including electromagnetic effects.Computers＆Geosciences，2009，35（2）：255－266.

［14］徐凱軍.2.5維復電阻率電磁場正反演研究［博士論文］.長春：吉林大學地球探測科學與技術學院，2007.Xu K J.Study on 2.5DComplex Resistivity Electromagnetic Forward and Inversion［Ph.D.thesis］（in Chinese）.Changchun：College of Geo－exploration Science and Technology，Jilin University，2007.

［15］梁盛軍.復電阻率法三維正反演問題研究［博士論文］.北京：中國地質大學，2011.Liang S J.Research on complex resistivity 3Dmodeling and inversing problem［Ph.D.thesis］（in Chinese）.Beijing：China University of Geosciences，2011.

［16］Constable S C，Parker R L，Constable C G.Occam′inversion：A practical algorithm for generating a smooth models form electromagnetic sounding data.Geophysics，1987，52（3）：289－300.

［17］De Groot－Hedlin C，Constable S.Occam′s inversion to generate smooth，two－dimensional models from magnetotelluric data.Geophysics，1990，55（12）：1613－1624.

［18］De－Lugao P P，Wannamaker P E.Calculating the twodimensional magnetotelluric jacobian in finite elements using reciprocity.Geophysical Journal International，1996，127（3）：806－810.

［19］Farquharson C G，Oldenburg D W.Approximate sensitivities for the electromagnetic inverse problem.Geophysical Journal International，1996，126（1）：235－252.

［20］Chen J P，Oldenburga D W，Haber E.Reciprocity in electromagnetics：application to modelling marine magnetometric resistivity data.Physics of the Earth and Planetary Interiors，2005，105：45－61.

［21］Mitsuhata Y.2－D electromagnetic modeling by finite－element method with a dipole source and topography.Geophysics，2000，65（2）：465－475.

［22］徐世浙.地球物理中的有限單元法.北京：科學出版社，1994：29－47.Xu S Z.FEM in Geophysics（in Chinese）.Beijing：Science Press，1994：29－47.

［23］Fan C S，Li T L.The simulation of 2.5Dcomplex resistivity model with finite element method.GEM Beijing 2011.Global Mtg.Expanded Abstracts，2011，15：47－50.

［24］羅延鐘，張桂青.頻率域激電法原理.北京：地質出版社，1988.Luo Y Z，Zhang G Q.Theory and Application of SIP （in Chinese）.Beijing：Geological Publishing House，1988.

［25］許文藝，徐兆文，顧連興等.安徽沙溪斑巖銅（金）礦床成巖成礦熱歷史探討.地質評論，1999，45（4）：361－367.Xu W Y，Xu Z W，F(xiàn)u L X，et al.Heat evolution from intrusion to mineralization in Shaxi porphyry copper（gold）deposits，Anhui Province.Geological Review （in Chinese），1999，45（4）：361－367.