程輝， 崔峻卿， 付國紅， 傅崧原

1 湖南科技大學先進礦山裝備教育部工程研究中心， 湘潭 411201 2 桂林理工大學廣西隱伏金屬礦產(chǎn)勘查重點實驗室， 桂林 541006

0 引言

巖礦石電性參數(shù)研究是電法(電磁法)領域的一項重要基礎性工作，其研究程度直接影響到方法理論研究與勘查實踐工作．眾多學者進行了卓有成效的研究工作．Wait(1959a，b)、Madden 和Cantwell(1967)、Ward和Fraser(1976)、Zonge等(1972)、Zonge和Wynn(1975)、Pelton(1977)、Pelton等(1978)、Dias(2000)、Gurin等(2015)、Zhdanov等(2018)、Gurin等(2018)進行了大量露頭、標本以及人造標本的測定工作，先后提出了不同的導電機理與導電模型以及相應的參數(shù)定義，如阻容模型、德拜模型、瓦爾堡模型、Cole-Cole模型、Dias模型、Gemtip模型等；何繼善(2006)、張賽珍等(1984)、王慶乙和徐立忠(2012)、底青云等(2008，2009)、薛國強等(2020)、陳儒軍等(2004)、黃理善等(2014)、柯式鎮(zhèn)等(2006)、李建軍等(2005)、趙云山等(2015)、何蘭芳(2014)、程輝等(2010，2014)、付國紅等(2020)進行了不同的實驗與理論研究工作．

以上研究多是基于變頻法進行的頻率特性測量，然而，使用變頻法觀測有著明顯不足(何繼善，2010)：(1)頻率改變前后的外部條件與干擾因素各不相同，抗干擾能力相對較弱；(2)各個頻點的觀測時間不同，觀測效率較低且增加了觀測數(shù)據(jù)一致性降低的可能性；(3)為了保證必要的穩(wěn)流條件，會增加發(fā)送端的復雜性．

為了達到實驗目的，研究人員使用了各類儀器設備進行觀測，包括阻抗分析儀(Solartron 1250、1260，Agilent4294A，HP4263B等)、LCR測試儀(數(shù)字電橋測試儀)(ZL5，Hioki 3522-50等)、巖心測試設備(LDT-10發(fā)射機與接收機組合形式、GDD公司生產(chǎn)的樣本巖芯IP測量儀SCIP等)、激電儀、電化學工作站等．雖然這些觀測裝置類型不一，但大多是基于變頻法思想逐個頻點觀測．

何繼善(2010)以2n序列偽隨機波形作為信號源應用于電法勘探，提出了雙頻激電法與廣域電磁法．近年來基于相關辨識的檢測技術(shù)應用到了電磁法勘探中，趙璧如等(2006)、萬健如等(2005)、羅維斌等(2012)、武欣等(2015)、羅延鐘等(2015，2020)、王顯祥等(2016)、王若等(2016)將m序列(或逆重復m序列)應用于電法勘探中，并進行相應的理論、實驗與野外對比觀測研究．若將相關辨識技術(shù)應用于巖礦石標本電性參數(shù)實驗，觀測數(shù)據(jù)的抗干擾能力將顯著提高，有利于提高測量精度與測量效率．

在針對巖礦石標本的室內(nèi)電性參數(shù)觀測工作中，巖礦石標本受其巖性、結(jié)構(gòu)與構(gòu)造、導電礦物類型與含量及其導電礦物空間分布與接觸關系、導電礦物與圍巖的接觸關系等多方面因素共同影響，即使同一類巖性的標本，在電性參數(shù)表現(xiàn)上也可能有數(shù)量級的區(qū)別．針對觀測對象的特殊性與復雜性，一般地，選用頻率域信號為激勵源進行實驗時，有以下要求：(1)實驗頻帶要足夠?qū)挘?2)為保證標本的電性參數(shù)為線性響應，流過標本的電流激勵信號幅值最小僅需維持在10-6A量級；(3)信號源要有較強的負載驅(qū)動能力；(4)為了保證觀測數(shù)據(jù)的可靠，要求發(fā)送端的電流輸出模塊必須保證具有較高的輸出穩(wěn)定度．

本文針對以上問題，提出基于逆重復m序列偽隨機電流信號的巖礦石標本電性參數(shù)觀測方案.研發(fā)了μA級編碼雙極性電流源，產(chǎn)生實驗所需的高精度小電流，測量頻帶寬度設置為10-3～104Hz，通過合理改變m序列的編碼參數(shù)，覆蓋整個測量頻帶，能有效抑制周期性干擾，滿足巖礦石標本電性參數(shù)觀測過程中的特殊要求．文中介紹了逆重復m序列的產(chǎn)生方法與頻譜特征，分析了不同階數(shù)序列的測量效率，針對工頻干擾提出編碼參數(shù)選擇要求；圍繞巖礦石標本電性特征，從觀測裝置角度論述了發(fā)送端驅(qū)動能力及接收端阻抗匹配、相關辨識法校正測試系統(tǒng)頻率響應以及高頻微弱電流信號的產(chǎn)生與拾取等關鍵技術(shù)問題，并對所研發(fā)的電流源進行了實驗驗證．

1 逆重復m序列的產(chǎn)生及其頻譜特征

1.1 逆重復m序列偽隨機信號的產(chǎn)生

逆重復m序列是將 2 個周期的m序列隔位取反后得到的一種偽隨機序列，并且與m序列互不相關．因此由m序列可方便得到逆重復m序列．

m序列偽隨機信號在電法勘探中采用二進制、雙極性的信號形式，n階m序列的周期長度為N=2n-1．m序列可由n級線性反饋移位寄存器生成，也可使用線性遞推公式產(chǎn)生：

ai=p1ai-1⊕p2ai-2⊕…⊕pnai-n,i≥0，

(1)

式中，⊕表示異或運算，pi為n階m序列的本原多項式系數(shù)．表1列出部分階數(shù)m序列的遞推公式．

表1 遞推公式表Table 1 Recursive formula

令初始條件a-1,a-2,…,a-n不全為 0 (一般取初始條件全為 1)，則m序列{ai} 可通過式(1)獲得．通過映射關系：

(2)

使用周期為Δt的時鐘控制ui輸出，則可獲得雙極性m序列偽隨機信號，周期為T=NΔt=N/fc(fc=1/Δt，為碼元頻率)．同樣地對逆重復m序列進行處理，可得到周期為2T的逆重復m序列偽隨機信號L(t)，圖1為fc=3.1 Hz的5階L(t)波形圖．

圖1 5階L(t)偽隨機信號波形(fc=3.1 Hz)Fig.1 Waveforms of 5th-order L(t) pseudo-random singal (fc=3.1 Hz)

1.2 逆重復m序列的頻譜特征

L(t)屬于循環(huán)平穩(wěn)隨機過程，其自相關函數(shù)為(鐘延炯，1979)：

RLL(τ)=

(3)

其中：

式(3)中第一項為三角脈沖函數(shù)，第二項為周期三角波函數(shù)，其自相關函數(shù)為三角脈沖函數(shù)與周期三角波之和，歸一化自相關函數(shù)圖如圖2a所示．由式(3)中各項的系數(shù)可知，階數(shù)n越大，周期長度N越大，周期三角波的幅值越小，自相關函數(shù)越接近脈沖函數(shù)，辨識誤差越?。虼诉x取適當大的階數(shù)n，可以提高辨識精度．

根據(jù)維納—辛欽定理，將L(t)的自相關函數(shù)RLL(τ)進行傅里葉變換即可得到L(t)的功率譜：

(4)

將式(3)代入式(4)，求得：

(5)

考察函數(shù) sinc(f/fc)，令sinc(f/fc)=1/2，求得f≈0.44fc，即 -3 dB帶寬截止頻率．

(6)

圖2 5階L(t)自相關函數(shù)圖(a)及理論功率譜圖(b)Fig.2 The 5th-order L(t)-signal autocorrelation function (a) and theoretical power spectrum (b)

kmax=0.44(2n-1)+0.5，

(7)

式中，kmax向左取整.

(8)

L(t)的周期公式為：

(9)

因此兩種方法所用時間之比為：

(10)

根據(jù)式(10)可知，測量效率與階數(shù)n有關，n越大，測量效率提升越明顯，表2列出了部分階數(shù)序列對應的測量頻點個數(shù)以及測量效率提高的倍數(shù)．

表2 不同階數(shù)L序列測量效率表Table 2 Efficiency of L-sequence measurement with different orders

2 抗工頻干擾參數(shù)選取條件

巖礦石頻譜測量工作中觀測數(shù)據(jù)受自身固有特性、外部因素、溶液條件、儀器和實驗條件的共同作用，從干擾類型與強度分析，工頻干擾最為明顯且難以屏蔽．由于L(t)具有良好的自相關性，通過合理選擇編碼參數(shù)，可以使工頻干擾信號正交于L(t)，從而通過相關辨識的方法壓制工頻干擾．

工頻干擾為典型的周期性干擾信號，在接收端表現(xiàn)為L(t)經(jīng)過巖礦石標本后的響應信號與干擾信號疊加．設周期干擾信號n(t)=Acos(ω0t+φ)，周期為T0=2π/ω0，其自相關函數(shù)為：

(11)

功率譜為：

Sn n(ω)=F[Rn n(τ)]

(12)

(13)

式中，f0=1/T0，其自功率譜也是離散譜，只在±f0處取值．設兩信號的互功率譜為SLn(f)，則有不等式：

(14)

RLn(τ)=F-1[SLn(f)]=0.

(15)

因此，當L(t)信號中不包含工頻頻率分量時，其與工頻干擾的互相關函數(shù)為0，可以有效壓制工頻干擾．考慮到實際測量過程中存在工頻干擾頻率偏移等問題，在選取碼元頻率fc時，可根據(jù)實際工作條件，避開工頻干擾的頻率．

圖3為抑制工頻干擾效果仿真圖，L(t)參數(shù)為：fc=2000 Hz，n=7，在50 Hz附近的頻點為39.37 Hz和55.12 Hz，波形如圖3a．經(jīng)過圖10所示的系統(tǒng)后，輸出響應為圖3b．引入信噪比為-20 dB的50 Hz干擾信號后，輸出響應如圖3c．未被干擾信號和被干擾信號的辨識結(jié)果如圖3d，兩曲線基本吻合，50 Hz附近無明顯波動，按照此要求選擇頻組，具有良好的工頻干擾抑制效果．

圖3 工頻干擾壓制仿真(a) 系統(tǒng)輸入信號; (b) 系統(tǒng)輸出響應; (c) 被50 Hz干擾的輸出響應; (d) 測得系統(tǒng)頻率響應對比.Fig.3 Suppression simulation of power frequency interference(a) System input signal; (b) System output response; (c) Output response disturbed by 50 Hz; (d) Comparison of measured system frequency responses.

3 觀測系統(tǒng)精細測量的解決方法

在實現(xiàn)精準的測量過程中，面臨電阻率分布跨度大、頻率響應存在觀測誤差大、頻率分布范圍寬等問題，在開展被測對象特征的研究基礎上，確定觀測系統(tǒng)精細測量的解決方法．

針對巖礦石標本電阻率分布跨度大(n×10～106Ωm)，采用高驅(qū)動能力設計及阻抗匹配技術(shù)；針對頻率響應存在觀測誤差，采用系統(tǒng)頻率響應校正技術(shù)；針對頻率分布范圍寬(10-3～104Hz)以及非線性效應問題，采用高頻微弱信號產(chǎn)生與拾取的相關技術(shù)．圖4為觀測系統(tǒng)示意圖，由發(fā)送機產(chǎn)生逆重復m序列電流信號，作為待測標本的激勵源，使用數(shù)據(jù)采集卡采集標本響應．

圖4 觀測系統(tǒng)示意圖Fig.4 Schematic diagram of the observation system

圖5 激發(fā)極化效應等效電路(a)及恒流源供電示意圖(b)Fig.5 Equivalent circuit of IP effect (a) and power supply of constant current source (b)

3.1 發(fā)送端驅(qū)動能力及接收端阻抗匹配問題

巖礦石標本電阻率可達n×106Ωm，以長l為5 cm，半徑r為2 cm的標準圓柱體標本為例，其直流阻值約為n×39.8 MΩ.W.H.Pelton經(jīng)過對巖礦石的大量實驗和研究，提出巖礦石激發(fā)極化效應可以等效為圖5a的等效電路．

此模型復阻抗可由式(16)計算得到：

(16)

恒流信號源通常采用壓控電流源的方式實現(xiàn)，其驅(qū)動能力主要受限于最高輸出電壓，原理如圖5b所示．恒流源輸出恒定電流，當負載RL發(fā)生變化時，輸出電壓UO亦隨負載變化．實際器件輸出電壓受供電電源限制，當負載足夠大時，輸出電壓達到最大值，此時恒流源將輸出恒定電壓，反饋回路失效，類似于恒壓源，電流將隨負載變大而減小，輸出信號將嚴重失真．另一方面，在頻率較高時，過大的負載將與線路間的雜散分布電容CP形成RC網(wǎng)絡，使τ值變大，導致輸出波形上升沿平緩，增大測量誤差．

圖6 測量回路阻抗等效電路Fig.6 Measurement circuit impedance equivalent circuit

對于接收端而言，測量回路阻抗等效電路如圖6所示，實際測得阻抗為負載RL與測量系統(tǒng)輸入阻抗Rin的并聯(lián)阻抗，測量時必須加以考慮．

因此在系統(tǒng)設計時，發(fā)送端要保證較高的輸出電壓范圍，選用高帶寬、高壓擺率的運放器件，并盡量降低回路分布電容，以擴大測量負載范圍，提高頻率響應范圍．接收端應采用差分輸入模式，并設置幾組不同的輸入阻抗，根據(jù)負載阻抗大小合理切換，在保證測量精度的前提下，保持較高的信噪比．

3.2 系統(tǒng)頻率響應校正

在實際工作中，盡管測試系統(tǒng)采用高性能器件，但系統(tǒng)頻率響應誤差仍然存在且不可避免．下面根據(jù)相關辨識測量原理，給出系統(tǒng)頻率響應誤差校正方法．相關辨識測量原理如圖7所示，L(t) 經(jīng)過測試系統(tǒng)H′(ω)后，引入干擾信號e(t)，最終測得輸出信號y(t)．然后將y(t)與L(t)通過前文第一小節(jié)的方法，分別計算出互功率譜和自功率譜并相除，即可求得系統(tǒng)響應H′(f)．

圖7 相關辨識原理圖Fig.7 Schematic diagram of correlation identification

測試系統(tǒng)H′(ω)簡化框圖如圖8所示．信號L(t)經(jīng)過發(fā)送通道的V-I轉(zhuǎn)換后變?yōu)殡娏餍图钚盘?，激發(fā)標本產(chǎn)生響應，后經(jīng)接收端測量通道進行采集，記錄為時間序列y[n]．

圖8 測量系統(tǒng)信號傳遞框圖Fig.8 Signal transmission block diagram of the measuring system

對測得數(shù)據(jù)處理后，得到的系統(tǒng)響應為發(fā)送通道響應h0(t)、標本響應h(t)和測量通道響應h1(t)共同作用的結(jié)果，最終測得系統(tǒng)頻率響應H′(ω)可用式(17)表示：

(17)

其中，H0(ω)、H(ω)、H1(ω)分別對應圖7中各部分的頻率響應，X(ω)、Y(ω)為系統(tǒng)輸入輸出采樣信號x[n]、y[n]的傅里葉變換．由式(18)可得標本真實頻率響應的修正公式：

(18)

根據(jù)式(18)可知，采用標準電阻作為負載，測得系統(tǒng)頻率響應即可修正誤差．

3.3 避免非線性效應的高頻弱信號產(chǎn)生及拾取

對于巖礦石的非線性效應，何繼善等(1995)提出不同礦物出現(xiàn)非線性特征所需的電流密度不同，最小僅1 μA/cm2，因此發(fā)送信號的電流幅值需實現(xiàn)μA級穩(wěn)定輸出．本文提出的觀測方法，所用信號采用雙極性形式，且碼元頻率fc是測量帶寬上限頻率的2.23倍，逆變頻率至少應達到22.3 kHz．微弱電流的高速逆變，不僅要保證電流的直流輸出精度，還要保證高頻時電流不會因寄生感抗而衰減，實現(xiàn)較為困難．經(jīng)過研究與多次試驗，采用圖9所示的運放級聯(lián)技術(shù)，實現(xiàn)了設計要求．輸入級采用帶寬盡量高的精密運放(對應圖9中的A1、A3)，輸出級采用高速運放(對應圖9中的A2)，保證精度與帶寬．輸入信號Uin為電壓信號L(t)，經(jīng)過運放電路實現(xiàn)電壓-電流轉(zhuǎn)換，將電壓信號轉(zhuǎn)換為電流信號輸出．輸出電流信號幅值由輸入信號幅值決定，通過調(diào)節(jié)輸入端電壓信號幅值獲得不同輸出電流幅值．

圖9 恒流源級聯(lián)原理示意圖Fig.9 Schematic diagram of the flow source cascade

接收端采集的信號是電壓信號，在測量低阻標本時，電壓通常只有n×100 μV，因此測量端至少需要保證0.1 μV測量分辨率．根據(jù)采樣定理，采樣頻率至少為44.6 kHz才能保證頻率不失真．在室內(nèi)環(huán)境中，工頻干擾耦合在測量回路中的幅值可達幾十毫伏，為了提高測量信噪比，在分辨率足夠的前提下，可在測量回路中串聯(lián)精密電阻提高輸出響應幅值，再根據(jù)基爾霍夫電壓定理(KVL)校正為標本響應．

4 實驗測試

為驗證方法的可行性與設計儀器的性能，實驗分別對所研發(fā)電流源樣機的直流恒流性能和交流測量精度進行測試，根據(jù)本文提出的測量方法，研制了信號發(fā)送樣機，電路板實物圖如圖10所示，采用NI公司的NI-PXIe4497數(shù)據(jù)采集卡采集數(shù)據(jù)，最大采樣率為204.8 kS/s，具有24位A/D轉(zhuǎn)換精度．

圖10 樣機電路板實物圖Fig.10 Real object diagram of prototype circuit board

直流恒流性能選取輸出電流幅值為2 μA，依次使用阻值為100 Ω、1 kΩ、10 kΩ、100 kΩ、1 MΩ的電阻作為負載，采用電流-電壓(I-V)轉(zhuǎn)換增益為104的I-V轉(zhuǎn)換電路，測量I-V轉(zhuǎn)換電路的輸出電壓，測量結(jié)果如圖11所示．在不同負載下，電流輸出穩(wěn)定度較高，整體波動范圍小于0.5%．從圖11中可以看到，在高阻負載情況下輸出電流偏低，是因為電阻存在制造公差，導致恒流源V-I轉(zhuǎn)換電路的反饋回路電阻阻值不嚴格匹配，可經(jīng)過調(diào)試進一步降低誤差，獲取更高的輸出穩(wěn)定度．

圖11 不同負載的穩(wěn)定性Fig.11 Stability under different loads

交流精度測試采用阻容模型模擬巖礦石標本，圖12給出了測量所用阻容模型及其幅頻特性曲線．

圖12 阻容模型及其幅頻特性Fig.12 Resistive and capacitive model and amplitude-frequency characteristics

圖13 轉(zhuǎn)折頻率附近頻組響應波形及功率譜(a) (e) fc=62 Hz； (b) (f) fc=625 Hz； (c) (g) fc=6.25 kHz； (d) (h) fc=62 kHz.Fig.13 Response waveforms of the frequency group near the turning frequency and power spectrum

圖14 頻率響應對比圖Fig.14 Comparison of frequency response

發(fā)送信號詳細參數(shù)為I=50 μA，n=5，fc=0.062 Hz、0.62 Hz、6.2 Hz、62 Hz、625 Hz、6.25 kHz、62 kHz．測量頻帶范圍為0.001～11 kHz，按十倍頻設置頻組，兩相鄰頻組間有重復頻段，共包含90個頻點如表4所示，實現(xiàn)單次單周期覆蓋測量頻帶用時約18.5 min，采用變頻法單次單周期測量以上頻點需要耗時約42.6 min，與變頻法相比測量效率提高2.3倍．圖13展示了轉(zhuǎn)折頻率附近頻組的波形及其功率譜，圖13a—d依次為fc=62 Hz、625 Hz、6.25 kHz、62 kHz的時域響應波形圖，圖13e—h為對應的功率譜圖．

圖14為測得頻率響應與理論值對比圖，從圖中可知測得幅頻曲線在工頻附近無明顯變化．頻率較大時存在系統(tǒng)誤差，按照式(18)的修正方法，計算出系統(tǒng)頻率響應并標定后，另取fc=25 kHz的頻組重新測量，測得響應呈現(xiàn)于圖14中，測量誤差如表3所示，標定后的測量誤差小于0.5%．

表3 測量誤差表Table 3 Measurement errors

5 結(jié)論與展望

本文基于逆重復m序列的編碼原理，論述了逆重復m信號的功率譜特點，推導了抗工頻干擾的參數(shù)選擇方法．從觀測裝置角度詳細分析了驅(qū)動能力及阻抗匹配、測試系統(tǒng)的頻率響應、以及高頻微弱信號的產(chǎn)生與拾取等關鍵問題．以5階序列為例的實驗結(jié)果表明：

(1)該方法可在10-3～104kHz頻帶內(nèi)測量90個有效頻點，實現(xiàn)寬頻帶、多頻點的觀測目的；同時，可根據(jù)需要改變逆重復m信號參數(shù)加密頻點，獲得更高的頻率分辨率．

(2)抗工頻干擾能力強，通過合理選擇信號參數(shù)，使發(fā)送信號不包含工頻分量，可以有效壓制工頻干擾．

(3)測量精度高，通過系統(tǒng)響應修正公式對系統(tǒng)進行標定后，測量誤差小于0.5%．

表4 5階信號在測量頻帶內(nèi)包含的頻點Table 4 Frequency points of 5th-order signal contained in the measurement frequency band

(4)測量效率高，測量效率由逆重復m序列階數(shù)決定，5階序列的測量效率相比變頻法測量提高2.3倍，在實際觀測工作中可根據(jù)需要適當提高階數(shù)以獲得更高測量效率．

(5)恒流源輸出采用分檔連續(xù)可調(diào)方案，在μA 級精度上可實現(xiàn)以1 μA為單位的輸出電流調(diào)整.

(6)觀測數(shù)據(jù)中有效頻點增加，為巖礦石電性參數(shù)計算與分析提供了基礎數(shù)據(jù)保證.

通過理論分析與實驗結(jié)果，可以展望下一步工作重點：

(1)接收系統(tǒng)通過MN拾取大于3 dB有效信號后利用虛擬儀器技術(shù)上傳至上位機(PC機)對時間波形保存，以便于進行后續(xù)計算與分析.

(2)通過時間波形的時頻計算，可在一次測量中，分別求取到時間域與頻率域電性參數(shù)的各類曲線，同時可利用導電模型，得到模型的各二級譜參數(shù)，進一步增強觀測數(shù)據(jù)的通用性.

感謝在文章撰寫過程，得到中國科學院地質(zhì)與地球物理研究所底青云研究員的指導，在此表示感謝.