井下核磁共振測(cè)井儀微弱信號(hào)處理算法設(shè)計(jì)

宋公仆 張嘉偉 薛志波 王光偉

（中海油田服務(wù)股份有限公司油技事業(yè)部 油田技術(shù)研究院，中國(guó) 北京 101149）

0 引言

核磁共振測(cè)井方法可直接測(cè)量地層孔隙中可動(dòng)流體的信息，可定量確定自由流體、束縛水、滲透率及孔徑分布，其孔隙測(cè)量不受巖石骨架礦物成分的影響，因此目前頗受測(cè)井行業(yè)的廣泛應(yīng)用[1]。與此同時(shí)，通過(guò)核磁共振測(cè)井采集到的地層回波信號(hào)相當(dāng)微弱，其幅度值大概在40nV～2uV之間，幾乎不可能通過(guò)常規(guī)測(cè)試手段直接檢測(cè)到相應(yīng)信號(hào)。本文重點(diǎn)介紹了兩種微弱信號(hào)的提取算法，通過(guò)這兩種算法可以有效地從采集到的微弱回波信號(hào)中快速提取相關(guān)信號(hào)幅度與相位信息，因此目前普遍應(yīng)用于核磁共振微弱信號(hào)的提取處理上。

1 微弱信號(hào)處理算法分析

1.1 相敏檢波算法（DPSD）分析

目前，對(duì)于納伏級(jí)（nV）微弱信號(hào)電壓的測(cè)量?jī)x器主要有鎖定放大器和取樣積分器兩類(lèi)。前者是物質(zhì)表面組份分析和表面電子能態(tài)研究的重要手段；后者使得核磁共振技術(shù)得以真正實(shí)現(xiàn)。這兩類(lèi)微弱信號(hào)檢測(cè)儀器可測(cè)量到淹沒(méi)在強(qiáng)噪聲中的μV～nV量級(jí)的電壓信號(hào)[2]。不過(guò)由于測(cè)量系統(tǒng)以及各種外界干擾引入的各種噪聲使得噪聲強(qiáng)度高出有用信號(hào)幾十倍，一些常規(guī)的電壓測(cè)量方法無(wú)法測(cè)量淹沒(méi)在強(qiáng)噪聲中的電壓信號(hào)。對(duì)于這些信噪比低至-30～-60db范圍的混于噪聲中的微弱信號(hào)的測(cè)量均采用基于最大似然估計(jì)的互相關(guān)方法[3]。下面重點(diǎn)對(duì)互相關(guān)法的一種DPSD相敏檢波算法進(jìn)行詳細(xì)介紹。

圖1 雙通道互相關(guān)器檢測(cè)方法-DPSD算法結(jié)構(gòu)圖

圖1中清晰列出了相敏檢波算法的結(jié)構(gòu)圖，其中設(shè)被測(cè)信號(hào)x（t）為：

式（1）中 K 為測(cè)量電路放大器增益，n（t）為零均值高斯噪聲。圖1中對(duì)應(yīng)的輸出值為對(duì)正弦信號(hào)x（t）的幅值Us以及相位θ的最大似然估計(jì)。其中

式（2）、（3）中參考信號(hào) r1（t）及 r2（t）為相互正交且與 x（t）同頻的正弦信號(hào)。 r1（t）=Urcosωt，r2（t）=Ursinωt。

已經(jīng)證明，由最大似然估計(jì)得到的正弦信號(hào)幅度及相位屬于一致、有效估計(jì)。但是當(dāng)被測(cè)正弦信號(hào)為納伏量級(jí)時(shí)，測(cè)量系統(tǒng)噪聲是信噪比低至-30db～-60db，則在有限測(cè)量時(shí)間T以?xún)?nèi)，U01和U02的測(cè)量結(jié)果會(huì)有較大的起伏波動(dòng)，從而嚴(yán)重妨礙了正弦參量Us以及θ的精確測(cè)量。顯然，如何選擇合適的測(cè)量電路以及參數(shù)是微弱正弦信號(hào)檢測(cè)的關(guān)鍵技術(shù)[4]。

而DPSD算法（數(shù)字相敏檢波算法）是指利用計(jì)算機(jī)或DSP芯片實(shí)現(xiàn)式（2）以及（3）的相敏解調(diào)算法。將以上二式數(shù)字化后，可寫(xiě)成

式中的 U01（n）為第 n 次取樣后的計(jì)算平均值。 r2（n）與 r1（n）相位差90°（如圖 1），故r2（n）可以通過(guò) r1（n）來(lái)產(chǎn)生，但由于 r1（n）是取樣信號(hào)，故要求r1（n）一周期內(nèi)的取樣點(diǎn)數(shù)與原頻率之比為4的整倍數(shù)，即m為4的整倍數(shù)。fs為采樣頻率（r（n）采樣頻率必須等于x（n）采樣頻率）。在預(yù)先存儲(chǔ)中可以先讓Ur為1，最后通過(guò)計(jì)算得出的相位是采樣起始點(diǎn)處的初相[5-6]。最后利用式（4）和（5）求出回波的幅值和相位。

以下是參考信號(hào) r1（n）與 r2（n）和 U（n）02的相關(guān)計(jì)算公式。

為了驗(yàn)證以上算法我們輸入標(biāo)準(zhǔn)正弦信號(hào)并加入噪聲后如下圖2左邊所示，可以看出信號(hào)淹沒(méi)在噪聲中，利用DPSD算法進(jìn)行提取后其提取結(jié)果如圖2右邊所示：

圖2 DPSD算法提取信號(hào)試驗(yàn)測(cè)試結(jié)果

從上圖2中，我們可以清晰地看到原始輸入信號(hào)混入白噪聲后形狀已經(jīng)失真，無(wú)法直觀辨別輸入信號(hào)的實(shí)際形狀。將混入噪聲的輸入信號(hào)進(jìn)行上述DPSD相敏檢波算法處理后得到圖2右邊所示的原始信號(hào)波形。從圖中我們可以直觀看到通過(guò)DPSD算法處理后，利用該算法選頻特性好的特點(diǎn)成功地將信號(hào)從噪聲中提取出來(lái)。

1.2 小波提取算法設(shè)計(jì)

小波變換在當(dāng)前微弱信號(hào)處理中應(yīng)用越來(lái)越多，它是時(shí)間（空間）頻率的局部化分析，通過(guò)伸縮平移運(yùn)算對(duì)信號(hào)（函數(shù)）逐步進(jìn)行多尺度細(xì)化，最終達(dá)到高頻處時(shí)間細(xì)分，低頻處頻率細(xì)分，能自動(dòng)適應(yīng)時(shí)頻信號(hào)分析的要求，從而可聚焦到信號(hào)的任意細(xì)節(jié)，解決了Fourier變換的困難問(wèn)題，成為繼Fourier變換以來(lái)在科學(xué)方法上的重大突破[3]。有人把小波變換稱(chēng)為“數(shù)學(xué)顯微鏡”。

1.2.1 含噪信號(hào)的小波域表述

假定待提取的信號(hào)是一個(gè)有用信號(hào)與噪聲的疊加效果，通過(guò)一種正交小波變換，可以最大程度的消除了信號(hào)的相關(guān)性，將能量集中到少數(shù)的小波系數(shù)上，從而達(dá)到濾除噪聲的目的。首先來(lái)分析其中噪聲小波域的分布規(guī)律。平穩(wěn)白噪聲的正交小波變換仍然是平穩(wěn)的白噪聲，其小波系數(shù)仍是互不相關(guān)的，分布在各個(gè)尺度下的所有時(shí)間軸上。假設(shè) n（t）是一個(gè)方差為 σ2的寬平穩(wěn)白噪聲，ψ（t）是一個(gè)小波函數(shù)，則白噪聲n（t）的小波變換的期望值為：

上式（11）中Ws為小波變換系數(shù)因子。上式（11）表明，隨小波變換尺度S的增加，白噪聲的小波變換幅值的均值E{|Wsn（t）|2}在減少，其衰減正比于1/S。若白噪聲n（t）是高斯白噪聲，在變換尺度S上，其小波變換模的平均密度為：

上式（12）說(shuō)明，高斯白噪聲的小波變換模值的平均密度正比于1/S，隨著分解尺度S的增大，其密度減小。由上述兩式性質(zhì)可知，隨著分解尺度S的增加，噪聲的小波譜因幅度和密度的逐漸減少而將逐漸消失，從而達(dá)到濾出噪聲的目的。下面來(lái)分析信號(hào)的小波域表述。設(shè)信號(hào)f（t）是光滑且P次連續(xù)可微的函數(shù)，ψ（t）是具有p階消失矩的實(shí)正交小波，其支撐為[c，d]，有下式（13）所示關(guān)系：

若上式（13）中 j足夠大，ψjk（t）變得足夠窄，將 f（t+2-jk）用其在 2-jk處的p階泰勒級(jí)數(shù)展開(kāi)近似表示，即：

1.2.2 提升小波算法

小波分析在實(shí)際工程中的應(yīng)用主要通過(guò)以下三種方式：MALLAT算法、多孔算法以及提升算法[5]。前兩種算法均建立在經(jīng)典小波分析基礎(chǔ)之上，其核心是基于二進(jìn)平移和伸縮思想的小波變換和多分辨率分析，而經(jīng)典小波分析則由傅里葉變換發(fā)展而來(lái)，因此MALLAT算法和多孔算法或多或少都會(huì)受傅里葉變換的影響。而提升小波算法相比基于經(jīng)典多分辨率分析的小波（第一代小波）最主要優(yōu)點(diǎn)是：改進(jìn)了第一代小波算法，使其更易于實(shí)現(xiàn)。提升算法不存在卷積運(yùn)算，節(jié)省了大量的計(jì)算時(shí)間和程序空間，計(jì)算量只有MALLAT算法的一半且避免了卷積計(jì)算存在的邊界問(wèn)題；計(jì)算過(guò)程簡(jiǎn)單明確，易于實(shí)現(xiàn)，可以采用原位計(jì)算，只需申請(qǐng)少量的額外存儲(chǔ)單元，節(jié)省了存儲(chǔ)空間；擺脫了傅里葉變換的束縛，可在第一代小波基的基礎(chǔ)上較為方便的設(shè)計(jì)滿足需要的小波基，盡管不能發(fā)現(xiàn)新的小波基，但可以改善現(xiàn)有小波基的特性，如消失矩特性等。小波濾波器組中濾波器h可以表示為式（16）所示的多相位形式，式中 he（z）和 ho（z）如式（17）所示，分別包含了 h（z）的偶系數(shù)和奇系數(shù)[7]。

上式（18）中稱(chēng)P（z）為多相位矩陣。小波分解重構(gòu)的多相位表示如圖3所示。

圖3 小波分解重構(gòu)的多相位表示

Daubechies和Swedens在此基礎(chǔ)上提出了多相位矩陣因子分解的定理：若 P（z）的行列式等于 1，則總存在 Laurent多項(xiàng)式 ui（z）和 pi（z）以及非零常數(shù) K，使得

上式（19）中pm（z）=0。根據(jù)多相位矩陣的因式分解以及圖4所示小波分解重構(gòu)的多相位表示，得到下式（20）和式（21）：

進(jìn)一步得到如圖4所示的小波分解過(guò)程的提升算法流程圖以及圖5所示的小波重構(gòu)過(guò)程的提升算法流程圖如下：

圖4 小波分解過(guò)程的提升算法流程圖

圖5 小波重構(gòu)過(guò)程的提升算法流程圖

上圖4所示為一層提升小波分解流程圖，得到近似系數(shù)a（z）與細(xì)節(jié)系數(shù)d（z），在反復(fù)按照?qǐng)D4所示進(jìn)行多層提升小波分解最終得到最底層的近似系數(shù)與細(xì)節(jié)系數(shù)。圖5為一層小波重構(gòu)過(guò)程的提升算法流程圖，它是提升小波分解的逆運(yùn)算，通過(guò)反復(fù)執(zhí)行提升小波重構(gòu)過(guò)程從而使分解的小波信號(hào)重構(gòu)成去噪后的有用信號(hào)。

2 結(jié)束語(yǔ)

本文論述了在核磁共振測(cè)井中微弱信號(hào)提取的兩種實(shí)現(xiàn)算法（DPSD）相敏檢波算法與小波提取算法，對(duì)這兩種算法進(jìn)行了詳細(xì)的介紹與分析。在核磁共振測(cè)井中，很多都利用這兩種提取算法對(duì)采集到的回波串信號(hào)進(jìn)行信號(hào)提取處理，其有著很好的噪聲抑制和信號(hào)提取的能力[8]。這兩種微弱信號(hào)提取算法已經(jīng)成功應(yīng)用在實(shí)際核磁共振測(cè)井儀器上，并取得了很好的應(yīng)用效果。

［1］George Coates，肖立志，Manfred Prammer.核磁共振測(cè)井原理與應(yīng)用[M].孟繁瑩，譯.北京：石油工業(yè)出版社，2007.

［2］邵維志，莊升，丁娛嬌.一種新型核磁共振測(cè)井儀：MREx[J].石油儀器，2004.

［3］Pollak V L and Slater R R.Input Circuit for PulsedNMR[J].The Review of Scientif ic Instruments，1966.

［4］華中科技大學(xué).微弱信號(hào)檢測(cè)技術(shù)資料[Z].

［5］戴逸松.微弱信號(hào)檢測(cè)方法及儀器[M].北京：國(guó)防工業(yè)出版社，1994.

［6］張新發(fā)，劉富，戴逸松.DPSD算法性能研究及參數(shù)選擇[J].吉林工業(yè)大學(xué)學(xué)報(bào)，1998，3（28）：40-45.

［7］戴逸松.測(cè)量低信噪比電壓的數(shù)字相敏解調(diào)算法及性能分析[J].計(jì)量學(xué)報(bào)，1997，18（2）：126-132.

［8］肖立志.核磁共振成像測(cè)井與巖石核磁共振及其應(yīng)用[M].北京：科學(xué)出版社，1998.