蔣川東 王琦 杜官峰 易曉峰 田寶鳳

1)(地球信息探測儀器教育部重點實驗室(吉林大學),長春 130061)

2)(吉林大學儀器科學與電氣工程學院,長春 130061)

3)(吉林大學通信工程學院,長春 130012)

1 引 言

地面核磁共振(surface nuclear magnetic resonance,SNMR)是一種以非侵入方式直接定量探測地下水水文參數(shù)的地球物理方法[1].與常見的核磁共振成像和核磁共振測井[2]等核磁共振技術(shù)的基本原理相同,SNMR具有對地下水含量分布成像和評估地層孔隙尺寸的能力[3].在地磁場(B0)中,地下水中的氫質(zhì)子(1H)的核磁矩繞B0方向旋進,并沿B0方向產(chǎn)生宏觀磁化矢量M,其角頻率ωL=|γB0|稱為Larmor角頻率,其中γ為1H磁旋比.通過在地表鋪設線圈并發(fā)射Larmor頻率(fL=ωL/2π)的交變磁場B1,M繞B1方向傾倒,產(chǎn)生橫向分量M⊥.當發(fā)射停止時,B1消失,M⊥繞B0做Larmor進動,同時在接收線圈中產(chǎn)生感應信號,稱為自由感應信號(free induction decay,FID)[4].通常發(fā)射頻率(fT)與fL相等,稱為共振激發(fā)(on-resonance excitation);當fT與fL不相等時,稱為偏共振激發(fā)(off-resonance excitation).在實際應用中,不可避免地會出現(xiàn)偏共振現(xiàn)象,如果近似為共振則會導致地下水成像結(jié)果不準確.因此本文以地面核磁偏共振(surface nuclear magnetic off-resonance,SNMOR)為主要研究目標,分析和討論偏共振激發(fā)時各種因素對SNMOR信號的影響.

Larmor頻率fL取決于B0的大小,影響B(tài)0的因素有以下四個方面：1)B0的不均勻性,由于地層介質(zhì)和水的磁化率不同[5],B0在地層內(nèi)部存在差異,近似為柯西-洛倫茲分布[6,7];2)B0的時變性,在同一個地點,B0的日平均變化量為70—120 nT,對應的fL變化在±5 Hz范圍內(nèi)[8];3)磁性介質(zhì)干擾.當測量區(qū)域存在鐵磁性礦產(chǎn)或大量金屬干擾時,B0會出現(xiàn)異常變化(＞100 nT),使得fT與fL的偏離量達到5 Hz以上[9];4)電磁噪聲,在強電磁噪聲環(huán)境中,無法準確獲得B0,利用頻率估計值fT發(fā)射,與fL存在偏差[10].因此,B0隨時間和空間變化、或受噪聲干擾都會出現(xiàn)測量不準確的情況,fT與fL不同將會發(fā)生SNMOR現(xiàn)象,利用儀器獲得的信號幅度和相位均不同于SNMR.如果采用傳統(tǒng)的共振模型對含水量和弛豫時間(T?2)等水文參數(shù)進行反演,將會嚴重影響結(jié)果的準確度.

針對SNMOR,Walbrecker等[11]首先分析了偏共振激發(fā)時頻率偏量?f=fT?fL對FID信號幅度和相位的影響.Grombacher等[12]利用偏共振現(xiàn)象,提出組合激發(fā)序列方式,通過優(yōu)化發(fā)射頻率提高FID信號的幅度.Irons和Li[13]提出基于洛倫茲變換的正反演模型,解決由B0不均勻性引起的FID信號失相(dephasing)問題.Grombacher等[14]討論了多種B0不均勻情況下的仿真模型,以及由此引起的含水量反演誤差;并提出了頻率環(huán)(frequency-cycling)激發(fā)序列,解決在fL未知時的探測問題[15].但是,以上研究主要集中在?f對SNMOR信號和反演結(jié)果的影響,沒有涉及利用?f修正SNMOR模型,進而提高含水量和T?2反演準確度的問題.目前,QT反演(pulse to time inversion,QTI)是最穩(wěn)定的反演方法[16],該方法只利用信號包絡的幅度信息,忽略其相位參數(shù).Roy和Lubczynski[17]提出信號初始幅度和初始相位同時參與的反演方法,在低噪聲環(huán)境情況下可提高對深部含水層的解釋準確度.Chen等[18]也提出了利用初始信號實部和虛部參與反演提高含水量解釋準確度的方法.但是,這些反演方法只利用FID信號的初始值,存在不穩(wěn)定的問題,而且只考慮了電阻率信息引起的相位偏量,沒有考慮?f引起的相位偏量.因此,基于?f建立SNMOR模型,同時利用FID信號包絡的實部和虛部進行復包絡反演,將會提高地下水反演準確度和穩(wěn)定性.

本文針對偏共振激發(fā)產(chǎn)生的SNMOR現(xiàn)象,分析SNMOR響應信號特征,提出自動搜索相位偏量的復包絡反演方法(complex envelope inversion,CEI).首先,從SNMR的基本原理入手,推導偏共振下的磁矩和靈敏度核函數(shù)的表達式,研究了?f對FID信號的影響;其次,提出相位偏量的自動搜索方法,實現(xiàn)CEI反演算法;第三,通過仿真模型分析SNMOR信號的特征和?f對反演結(jié)果的影響;最后,利用野外實測數(shù)據(jù)證明SNMOR模型和CEI的有效性和準確性.

2 SNMOR原理

根據(jù)NMR原理,FID信號表達式為[1]

其中,q=I0τp稱為激發(fā)脈沖矩,是發(fā)射電流強度I0和持續(xù)時間τp的乘積;V0是初始振幅,與含水量w(含水體積與總體積之比)成正比;T?2是橫向弛豫時間,與孔隙大小相關;?0是初始相位,與地下電阻率(ρ)、線圈空間位置、頻率偏差(?f)和儀器系統(tǒng)等多種因素有關.根據(jù)SNMR正演理論[19],信號初始振幅V0表示為

其中,K稱為靈敏度核函數(shù),即單位含水量產(chǎn)生的FID信號大小;r為探測體積vol內(nèi)的坐標.

2.1 激發(fā)磁化矢量

地面上的發(fā)射線圈產(chǎn)生的激發(fā)磁場B1(ωL,I0,ρ,r)是發(fā)射頻率ωL、電流強度I0、地下電阻率ρ和空間位置r的函數(shù),可利用有限元等方法計算[20]. 由于存在地磁傾角I和地磁偏角D,只有垂直于地磁場B0方向的分量B⊥對FID信號有貢獻,B⊥(ωL,I0,ρ,r,I,D)=RIRDB1(ωL,I0,ρ,r)(其 中RI和RD為 旋 轉(zhuǎn) 矩陣[10],以下簡記為B⊥).在地磁場B0和激發(fā)磁場B⊥作用下,氫質(zhì)子自旋產(chǎn)生的宏觀磁化矢量M用Bloch方程[21]表示:

其中,Beff代表有效激發(fā)磁場;?x,?y和?z代表以fT為頻率的旋轉(zhuǎn)坐標系的單位向量,?z的方向與B0相同;T1是縱向弛豫時間.Beff具體表示為

1)共振激發(fā)(fT=fL)

當?f=0 Hz時,Beff=B⊥即發(fā)生SNMR現(xiàn)象,B1在導電介質(zhì)中發(fā)生橢圓極化[22].由于平行于B0方向的分量B⊥z非常微弱,B⊥在旋轉(zhuǎn)坐標系下可以表示為

其中,?=?T+?σ,即發(fā)射電流相位?T與電導率引起的極化相位?σ之和.

在理想條件下,忽略激發(fā)期間的弛豫現(xiàn)象[23]和假設發(fā)射電流包絡為矩形脈沖,求解表達式(3),

此時,B⊥和τp決定了共振激發(fā)的扳倒角度θon=|γB⊥|τp. 通過控制?T=0, 得到在共振時單位磁化矢量Mon在垂直于B0方向(?z方向)的平面上(?x-?y平面)的分量M⊥,on為

因此,對于共振激發(fā),M⊥,on在?x方向上為0,只在?y方向上存在分量.

2)偏共振激發(fā)(fT?=fL)

當?f=0 Hz,發(fā)生SNMOR現(xiàn)象,由表達式(4),Beff不但與B⊥有關,還在方向上存在與?f相關的分量?B0=2π?f/γ.由表達式(3)描述的感應扭矩導致了在Beff方向上的章動,與平面的方向角度α=?B0/|B⊥|.此時有效的扳倒角度為單位磁化矢量Moff的垂直分量M⊥,off表示為

圖1(a)解釋了當?f=0 Hz(虛線)和4 Hz(實線)時,M⊥的單位方向矢量(即對M⊥幅度歸一化)在旋轉(zhuǎn)坐標系上的運動軌跡.在共振(?f=0 Hz)時,扳倒角度為θon,M⊥,on只在?y方向上存在分量.而對于偏共振(?f=0 Hz)時,扳倒角度為θoff,M⊥,off在?x和?y方向上均存在分量,因此偏共振激發(fā)產(chǎn)生了額外的相位偏差??f,見圖1(b)所示.

圖1 旋轉(zhuǎn)坐標系下共振激發(fā)和偏共振激發(fā)的磁化矢量軌跡圖(a)M在三維坐標上的運動軌跡;(b)M在?x-?y平面運動軌跡的二維投影Fig.1. Trajectory of magnetization in a rotating reference frame during on-resonance and off-resonance excitation:(a)Trajectory of M in three dimensional coordinates;(b)trajectory projection of M in two dimensional coordinates.

2.2 靈敏度核函數(shù)

根據(jù)表達式(1)和(2),當發(fā)射線圈和接收線圈共線時,接收磁場BR與激發(fā)磁場B1相等,K表示為[18]

其中,M⊥在共振時取M⊥,on,在偏共振時取M⊥,off;M0是凈磁化矢量;ξ是發(fā)射和接收磁場由于橢圓極化產(chǎn)生的相位分量.對于一維層狀大地模型,SNMOR的正演表達式可以化簡為

其中,V為FID信號的復包絡,z為深度.圖2分別計算了在SNMR和SNMOR模型中,K的實部和虛部隨z和q變化的結(jié)果.仿真正演模型為：B0為54721 nT,fT為2330 Hz,I為60°,D為0°,ρ為100 ?m均勻分布,線圈為邊長100 m的正方形,q設置為0.01—12 As,按對數(shù)分布,數(shù)量為24個,τp為40 ms.如圖2(a)—(c)所示,K的實部圖像在?f=?3,0和3 Hz時基本一致,主要信號(黃色部分)隨深度增加,所需q也相應增加,且幅度逐漸減小,說明含水層深度增加,探測難度逐漸增大.圖2(d)—(f)中K的虛部圖像在?f=?3,0和3 Hz時呈現(xiàn)完全不同的結(jié)果.對于共振,因為m⊥,on無額外的相位,K的虛部幾乎為0,其中較小的數(shù)值是由ρ引起的?ρ,可以忽略.而對于偏共振,M⊥,off在Mx和My均有分量,因此產(chǎn)生了額外的相位??f,K的虛部信號明顯增加,而且?f=?3 Hz(藍色部分)和3 Hz(黃色部分)的符號相反.因此,相比于SNMR,SNMOR能夠獲得額外的虛部信息,更強的信號幅度和更多的相位信息有利于反演地下的含水層分布.

圖2 在SNMR和SNMOR模型下,K隨q和z變化的圖像 (a)?f=?3 Hz,(b)?f=0 Hz,(c)?f=3 Hz時K的實部;(d)?f=?3 Hz,(e)?f=0 Hz,(f)?f=3 Hz時K的虛部Fig.2.Complex kernel function K with q and z in SNMR and SNMOR models:(a)?f=?3 Hz,(b)?f=0 Hz,(c)?f=3 Hz,the real part of K;(d)?f=?3 Hz,(e)?f=0 Hz,(f)?f=3 Hz,the imagery part of K.

3 復包絡反演方法

由SNMOR的原理和正演模型得到接收線圈獲得的V是復數(shù),其相位與發(fā)射電流相位?T(可控制為0 rad)、電阻率相位?ρ、頻率偏量相位??f相關,還與發(fā)射線圈和接收線圈的位置變化引起的幾何相位?g、儀器系統(tǒng)的接收電路引入的相位?c和數(shù)據(jù)處理時產(chǎn)生的相位?d相關.對于一維地下水探測,發(fā)射和接收線圈共線,?g可以忽略.而且,采用相同的儀器和數(shù)據(jù)處理方法,?c和?d為固定的常數(shù)?s=?c+?d,統(tǒng)稱為系統(tǒng)相位.因此,接收信號V的相位可以表示為

復包絡反演的核心是利用表達式(10)中接收信號V(q,t)的實部和虛部,根據(jù)已知的ρ和?f,正演建模計算K,然后通過自動搜索匹配未知的?s,多次迭代實現(xiàn)反演模型計算數(shù)據(jù)與實測數(shù)據(jù)相互吻合.

復包絡反演算法的示意圖見圖3,V的實部和虛部組成測量數(shù)據(jù)d矩陣：

其中,dq,p是第q個脈沖矩、第p個時刻的接收信號V(q,tp),R代表信號實部,I代表信號虛部,Q代表q的個數(shù),P代表時間序列個數(shù),H代表轉(zhuǎn)置.

CEI的變量為每層地下水的含水量w和(簡寫為T),組成變量矩陣m:

其中,wn代表第n層的w,Tn代表第n層的N代表地下分層數(shù).

根據(jù)表達式(9)和(10),假設地下ρ分布和?f已知,正演獲得當不存在?s時,d=Kw;當存在?s且未知時,可以先通過搜索(G(·)為模型正演算子,m0為變量初值)獲得初步的相位估計然后,更新的核函數(shù)

采用迭代方法計算反演變量mk+1=mk+?mk.利用表達式(10)分別計算V對w和兩個變量的導數(shù)(Jacobian)矩陣：

和

為了將數(shù)據(jù)的實部和虛部同時參與反演,將Jacobian矩陣整理成：

為了將變量m限制在合理范圍內(nèi)(0＜w＜100%,0＜T?2＜1 s),利用tan函數(shù)將約束空間轉(zhuǎn)換成非約束空間,同時更新J[16].本文采用基于約束加權(quán)和先驗的共軛梯度最小二乘(CGLSCDP)算法[23,24]進行模型增量?m求解：

其中,Dε是觀測數(shù)據(jù)的噪聲估計加權(quán)矩陣,C是模型平滑度矩陣,λ是正則化參數(shù),可利用χ2分布檢驗或L曲線優(yōu)化選取[25].

在每次迭代過程中,需要選擇最優(yōu)的迭代步長ηk以防止模型發(fā)散：

當確定最優(yōu)步長?ηk后,重新搜索匹配相位??k+1,用于確定下一次迭代的相位偏差：

通過多次迭代后,m和??均趨于穩(wěn)定,由反演模型計算數(shù)據(jù)與測量數(shù)據(jù)和噪聲數(shù)據(jù)的加權(quán)均方差χ2趨近于1時,最終獲得反演結(jié)果.

圖3 復包絡反演算法示意圖Fig.3.Schematic of complex envelop inversion algorithm.

4 仿真模型結(jié)果

4.1 SNMOR信號

為了分析?f對SNMOR信號的影響,建立包含三個含水層的仿真模型,厚度均為25 m,如圖4(a)和圖4(b)所示;第一層深度為0—25 m,含水量為5%,T2?為100 ms,對應細砂層;第二層深度為25—50 m,含水量為25%,T2?為200 ms,對應粗砂層;第三層深度為50—75 m,含水量為10%,T2?為50 ms,對應黏土層.SNMOR正演參數(shù)與圖2中計算K的參數(shù)一致.

利用正演表達式(9)和(10)分別計算?f=±10 Hz范圍內(nèi)的信號V0,即t=0 s時刻的信號幅度,如圖4(c)—(e)所示.圖4(c)是V0的幅度,當q＜2 As時,?f對V0的幅度值影響較小(＜7.5%);當q＞2 As時,?f對V0的幅度值影響逐漸增大.具體表現(xiàn)為SNMR時,V0最小(綠色曲線);SNMOR情況下,隨著?f變化到±10 Hz,V0也逐漸增大(藍色系和紅色系曲線),最大增加了65.9%.說明對于較大的q,SNMOR與SNMR信號幅度明顯不同,因此不能忽略?f的影響.由圖4(d)和圖4(e)得到,當?f由0 Hz變化到±10 Hz時,V0的實部逐漸減小,而V0的虛部逐漸增大,但是對于正?f(紅色系曲線)和負?f(藍色系曲線),V0虛部的符號相反.因此,?f既影響FID信號的幅度,也影響其相位,與圖2中K的結(jié)果一致.

圖4 SNMOR仿真模型和正演計算信號的初始振幅V0(a)仿真模型含水量;(b)仿真模型T?2;(c)正演計算信號V0的幅度;(d)V0的實部;(e)V0的虛部Fig.4.Simulated SNMOR model and forward calculation of the initial values V0of the signal:(a)Water content of model;(b)T?2of model;(c)amplitude of V0;(d)real part of V0;(e)imagery part of V0.

4.2 SNMOR反演結(jié)果

對于圖4中的仿真模型和SNMOR信號(模型參數(shù)見2.2小節(jié)),利用第3節(jié)的CEI,獲得含水量和T?2分布結(jié)果如圖5所示.為了仿真野外環(huán)境,在SNMOR信號中加入了50 nV高斯分布的隨機噪聲,并將FID信號進行了抽道,抽道數(shù)為30.在反演過程中,假設初始模型為同一含水層,含水量和T?2分別為10%和100 ms,最大值分別設定為50%和350 ms,正則化參數(shù)選擇為1000,經(jīng)過50次迭代后,所有的反演數(shù)據(jù)與仿真數(shù)據(jù)的誤差統(tǒng)計χ2≈1.

在復包絡反演中,利用了SNMR核函數(shù)和SNMOR核函數(shù)進行了對比.通常,在地磁場測量不準確或者分布不均勻時,只能假設為SNMR情況,即忽略?f的影響.由圖5(a)和圖5(c)中的含水量和T?2反演結(jié)果可知,假設SNMR情況時,?f的變化對于第一層含水層影響較小,這是因為小脈沖矩時FID信號隨?f變化較小,可以忽略.但是對于第二層和第三層含水層,?f的影響較大.當?f＜0 Hz時(藍色系曲線),第二層含水層的反演結(jié)果出現(xiàn)先偏小后偏大的趨勢,含水量最大偏差18.4%,T?2最大偏差56.1 ms.當?f＞0 Hz時(紅色系曲線),第二層和第三層含水層的反演結(jié)果出現(xiàn)先偏大后偏小再偏大的振蕩趨勢,含水量最大偏差36.1%,T?2最大偏差275.6 ms.由此可以得到,對于較大脈沖矩,FID信號隨?f變化大,對中部和深部含水層的反演結(jié)果影響也較大.因此,為了準確探測中部和深部的含水層,不能忽略偏共振的影響.

當準確獲得?f時,可以計算SNMOR核函數(shù)應用到復包絡反演中,見圖5(b)和圖5(d).當?f從?10 Hz變化到10 Hz時,對于第一層和第二層含水層,反演獲得的含水量和T?2分布均與仿真模型基本一致,只在層邊界處由于平滑過渡的假設出現(xiàn)緩慢變化;對于第三層含水層,隨著?f不同,反演結(jié)果與仿真模型稍有不同,其中含水量最大偏差4.2%,T?2最大偏差39.3 ms.第三層出現(xiàn)的偏差是因為隨深度增加,反演分辨率下降,尤其對于大于50 m的反演結(jié)果,分辨率較低.對比SNMR核函數(shù)和SNMOR核函數(shù)的復包絡反演結(jié)果得到結(jié)論,當頻率偏量?f存在時,利用?f計算的SNMOR核函數(shù)進行反演,可以顯著減小?f對含水量和T?2反演結(jié)果的影響.

圖5 不同?f情況下的復包絡反演結(jié)果 (a)SNMR核函數(shù)反演的含水量;(b)SNMOR核函數(shù)反演的含水量;(c)SNMR核函數(shù)反演的T?2;(d)SNMOR核函數(shù)反演的T?2Fig.5.Inversion results with different ?f:(a)Water content using SNMR kernel;(b)water content using SNMOR kernel;(c)T?2using SNMR kernel;(d)T?2using SNMOR kernel.

4.3 CEI與QTI對比

為了對比CEI與QTI的結(jié)果,利用圖5中的仿真模型和當?f=3 Hz時的正演FID信號,分別加入噪聲水平為50—500 nV高斯分布的隨機噪聲,進行反演對比實驗,反演結(jié)果如圖6所示.

兩種反演方法都利用?f=3 Hz時的SNMOR核函數(shù)K,但是QTI方法計算核函數(shù)的模|K|,只有信號幅度參與計算,隨著噪聲水平增加,含水量和T?2出現(xiàn)偏差,如圖6(a)和圖6(c)所示.對于第一層含水層,隨著噪聲水平增加,QTI結(jié)果在淺層出現(xiàn)了較大偏差,含水量和T?2最大偏差分別為19.9%和177.8 ms.第二層含水層,含水量隨著噪聲增大先減小后增大,T?2隨噪聲增大而增大.對于第三層含水層,含水量和T?2均隨著噪聲增大而增大,最大偏差分別為14.1%和205.0 ms.QTI結(jié)果說明,隨著噪聲增加,模型分辨率降低,尤其是噪聲較大時(＞300 nV),QTI不能區(qū)分第二層和第三層含水層.相比之下,CEI同時利用核函數(shù)和信號的實部和虛部,結(jié)果更加穩(wěn)定和可靠.對于不同噪聲情況,在第一層含水層中均未出現(xiàn)假含水層,含水量和T?2最大偏差分別為1.8%和53.3 ms;第二層和第三層的含水量和T?2分布均比QTI結(jié)果的偏差小,說明CEI的模型分辨率優(yōu)于QTI.因此可以證明,在噪聲較小時,CEI和QTI結(jié)果無明顯差別;隨著噪聲增加,CEI的優(yōu)勢逐漸明顯.但是當噪聲很大時,兩者的反演結(jié)果均不能反映仿真模型.

圖6 不同噪聲水平的CEI和QTI結(jié)果對比(?f=3 Hz)(a)QTI的含水量;(b)CEI的含水量;(c)QTI的T2?;(d)CEI的T2?Fig.6. Results of CEI and QTI with different noise levels(?f=3 Hz):(a)Water content using QTI;(b)water content using CEI;(c)T2?using QTI;(d)T2?using CEI.

5 野外實測數(shù)據(jù)

為了驗證SNMOR探測的有效性,利用JLMRS型地下水探測儀器在長春市郊的太平池水庫旁進行了野外實驗.該地點的環(huán)境噪聲較低,實驗環(huán)境較好,在此進行過多次SNMR和其他方法測量,并且有鉆孔資料進行對比驗證.綜合已知資料得到地下水分布為:1)0—13 m,黏土層,平均含水量為2%;2)13—25 m,細砂和砂礫互層,平均含水量12%;3)25—45 m,黏土層,平均含水量為4%;4)45—62 m,細砂和粘土互層,平均含水量7%.測量地點的平均B0為54720 nT,在一天內(nèi)變化小于12 nT,對應fT為2330 Hz±0.5 Hz,I為62.11°,測量地下ρ電阻率平均為30 ?m.SNMOR測量采用邊長為50 m的單匝正方形線圈,q為20個,在0.5—10 As之間按照對數(shù)分布,每組q下重復測量32次,獲得噪聲水平平均為11 nV.圖7給出了在q=2.4 As時,實測FID信號實部和虛部隨?f變化的曲線,其中藍色曲線是實測數(shù)據(jù),紅色曲線是擬合結(jié)果,表1總結(jié)了實測數(shù)據(jù)的擬合結(jié)果和不確定度.

由圖7和表1可以得到:在SNMR情況(?f=0 Hz)下,實測信號的實部和虛部均按照指數(shù)規(guī)律單調(diào)變化,T?2為107.2 ms,?f為0.13 Hz,可以忽略其影響. 在SNMOR情況下(?f=0 Hz),隨著?f的增加,實測信號的實部和虛部曲線均出現(xiàn)振蕩波動,其波動頻率與?f基本一致.在表1中,當?f從±1 Hz變化到±7 Hz時,V0逐漸增大,最大值比SNMR時增加了52.6%,T?2基本保持不變(95.5—102.6 ms),環(huán)境噪聲水平εn不變(9.3—13.6 nV).雖然曲線擬合的初始相位?0結(jié)果變化不規(guī)律,但是經(jīng)過后續(xù)復CEI后可以獲得對應的系統(tǒng)相位?s,兩者的相位差??呈現(xiàn)明顯的規(guī)律(如表1所列),即隨著?f的增加,??也逐漸增加.由以上分析可以證明,當fT與fL出現(xiàn)偏差時,發(fā)生SNMOR現(xiàn)象,現(xiàn)有的JLMRS型儀器即可采集到有效的SNMOR信號.而且,SNMOR信號的幅度大于SNMR信號幅度,與理論計算結(jié)果相符.

圖7 太平池水庫野外實測FID信號 (a)共振激發(fā),?f=0 Hz;(b)偏共振激發(fā),?f=1 Hz;(c)?f=?1 Hz;(d)?f=3 Hz;(e)?f=?3 Hz;(f)?f=5 Hz;(g)?f=?5 Hz;(h)?f=7 Hz;(i)?f=?7 HzFig.7. FID signal measured in Taipingchi Reservoir:(a)On-resonance excitation,?f=0 Hz;(b)off-resonance excitation,?f=1 Hz;(c)?f= ?1 Hz;(d)?f=3 Hz;(e)?f= ?3 Hz;(f)?f=5 Hz;(g)?f=?5 Hz;(h)?f=7 Hz;(i)?f=?7 Hz.

表1 實測FID信號數(shù)據(jù)擬合結(jié)果和不確定度(q=2.4 As)Table 1.Fitted parameters and uncertainty of the fi eld measured FID signals(q=2.4 As).

對于20組q下不同?f的SNMOR探測結(jié)果如圖8所示.其中圖8(a)是實測FID信號V0的幅度,圖8(b)和圖8(c)分別是經(jīng)過系統(tǒng)相位?s修正后V0的實部和虛部.圖8中的曲線隨?f變化的變化規(guī)律與圖4仿真結(jié)果的規(guī)律類似.對于圖8(a)中V0的幅度,SNMR時(綠色曲線)相對較小;而SNMOR時,隨著?f增加(藍色系和紅色系曲線),大脈沖矩(q＞2 As)對應的V0幅度逐漸增大.對于圖8(b)中的V0實部,SNMR時(綠色曲線)最大,而SNMOR時(藍色系和紅色系曲線)稍有下降.但是,對于圖8(c)中的V0虛部,SNMR時(綠色曲線)在0 nV左右,只有q＞7 As時,由于ρ的影響,略有增大.而在SNMOR時(藍色系和紅色系曲線),隨著正?f和負?f變大,V0虛部的數(shù)值均逐漸增加,且符號相反.因此,由實測結(jié)果可以證明本文建立的SNMOR正演模型是正確的.

圖8 實測FID信號初始振幅V0隨?f的變化曲線 (a)V0的幅度;(b)V0的實部;(c)V0的虛部Fig.8.Initial value V0of the FID signals:(a)The amplitude of V0;(b)the real part of V0;(c)the imagery part of V0.

利用本文提出的CEI,對不同?f的SNMOR實測數(shù)據(jù)進行反演得到結(jié)果如圖9所示.在反演過程中,假設初始模型為同一含水層,含水量為5%,T?2為100 ms,設定兩者的最大值分別為15%和300 ms,正則化參數(shù)根據(jù)χ分布檢驗選擇,經(jīng)過50次迭代,得到最終的數(shù)據(jù)誤差統(tǒng)計χ2≈1.圖9(a)和圖9(c)是不考慮?f的情況下,假設均采用SNMR核函數(shù)時的反演結(jié)果,可以得到只有當fT和fL近似相等時(綠色曲線),反演結(jié)果的含水量和T?2與已知資料相近,當頻率不相等時,反演誤差較大,如在淺層(＜5 m處)出現(xiàn)假含水層(藍色系曲線),對于第二層和第四層的含水量偏低,第三層含水量偏高(紅色系曲線)等.隨著?f增大,反演結(jié)果的誤差也隨之增大,當?f=±7 Hz時已經(jīng)不能反映地下含水層的真實分布情況.圖9(b)和圖9(d)是采用?f計算的SNMOR核函數(shù)的反演結(jié)果,與共振情況相比,含水量和T?2分布能夠正確反演地下含水層分布情況.但是,隨著頻率偏量的增加,反演結(jié)果的誤差也逐漸增大,尤其是對于?f=±7 Hz的藍色系和紅色系曲線,第二層和第三層的含水量和T?2的誤差較大,只能近似反映真實情況.因此,通過實測數(shù)據(jù)反演結(jié)果可以證明,SNMOR現(xiàn)象發(fā)生時,需要準確獲得?f,通過CEI才能獲得更準確的地下含水層分布.

圖9 不同?f時實測FID信號的CEI結(jié)果 (a)SNMR核函數(shù)反演的含水量;(b)SNMOR核函數(shù)反演的含水量;(c)SNMR核函數(shù)反演的T?2;(d)SNMOR核函數(shù)反演的T?2Fig.9.Inversion results of the fi eld measured FID signals with different?f:(a)Water content using SNMR kernel;(b)water content using SNMOR kernel;(c)T?2using SNMR kernel;(d)T?2using SNMOR kernel.

最后,利用不同噪聲水平下的實測數(shù)據(jù)對比反演結(jié)果,驗證CEI的優(yōu)勢.噪聲水平與疊加次數(shù)相關,32次、24次、16次和4次疊加數(shù)據(jù)分別獲得的噪聲水平為11,15,21和26 nV,由此可以獲得不同噪聲水平實測數(shù)據(jù)的反演結(jié)果,如圖10所示.圖10(a)和(c)是利用?f=3 Hz的核函數(shù)進行QTI的結(jié)果,可見隨著噪聲增大,反演結(jié)果逐漸偏離真實模型.當噪聲為11 nV時,除了T2?顯示出淺層的假含水層以外,其他結(jié)果與真實模型類似,但是當噪聲增大到20 nV以上后,含水量和T2?大幅度偏離了真實模型.圖10(b)和圖10(d)利用?f=3 Hz的核函數(shù)進行CEI的結(jié)果,當噪聲增大時,反演結(jié)果基本與真實模型一致,只有第四個含水層由于噪聲增加,分辨率下降,因此含水量和T2?逐漸偏小.由CEI和QTI的結(jié)果對比可以得到結(jié)論,在較大噪聲環(huán)境下,CEI充分利用數(shù)據(jù)的有用信息,能夠獲得比QTI更穩(wěn)定和可靠的結(jié)果.

6 結(jié) 論

SNMR方法是目前唯一一種直接探測地下水的地球物理方法.由于地磁場不均勻性,時變性和易受干擾等問題影響,在實際測量中經(jīng)常進行SNMOR探測.本文針對SNMOR中偏共振激發(fā)引起FID信號的幅度和相位變化,提出包絡信號實部和虛部同時參與的CEI,通過自動搜索方案估計未知的系統(tǒng)相位,實現(xiàn)了SNMOR探測中含水量和弛豫時間(T?2)的準確反演.最后,通過仿真模型和實測數(shù)據(jù)驗證了SNMOR正演模型和CEI的有效性和準確性.

圖10 不同噪聲水平的CEI和QTI結(jié)果對比(?f=3 Hz) (a)QTI的含水量;(b)CEI的含水量;(c)QTI的T?2;(d)CEI的T?2Fig.10. Results of CEI and QTI with different noise levels(?f=3 Hz):(a)Water content using QTI;(b)water content using CEI;(c)T?2using QTI;(d)T?2using CEI.

在SNMOR探測中,隨著頻率偏量的增加,FID信號初始振幅的實部稍有減小,而虛部的數(shù)值逐漸增加,其符號與頻率偏量的符號一致,說明信號相位受頻率偏量的影響較大.在本文的仿真模型中,信號幅度隨頻率偏量增加而增大,最大增幅達到65.9%.常規(guī)反演不考慮頻率偏量的影響,只利用SNMR核函數(shù),因此含水量和T?2的反演結(jié)果誤差較大.本文根據(jù)頻率偏量準確計算了SNMOR核函數(shù),利用CEI方法提高了含水量和T?2反演結(jié)果的準確度,在仿真模型中最大偏差為4.2%和39.3 ms.通過與QTI方法進行對比得到,CEI充分利用測量數(shù)據(jù)的全部信息,反演結(jié)果優(yōu)于QTI,當噪聲較大時,仍然能夠獲得更穩(wěn)定和可靠的反演結(jié)果.

本文通過SNMOR探測中實測FID信號的擬合結(jié)果獲得頻率偏量,可以修改發(fā)射頻率重新測量以保證共振激發(fā),也可以利用本文提出的偏共振模型和CEI反演方法直接獲得準確的地下水文信息.但是,如果頻率偏量較大,超過了交變磁場引起的章動頻率,將會導致無法有效激發(fā)氫質(zhì)子,從而造成較大的測量誤差.此時,必須調(diào)整發(fā)射頻率使其等于或接近Larmor頻率.另外,在測量過程中,有時頻率偏量也存在隨時間和空間變化,或者強噪聲影響測量結(jié)果.因此,需要進一步研究頻率偏量變化時的SNMOR反演問題,以及利用頻率環(huán)等測量序列解決Larmor頻率未知時的探測問題.

[1]Legtchenko A 2013Magnetic Resonance Imaging for Groundwater(Hoboken:John Wiley&Sons)p10

[2]Li X,Xiao L Z,Liu H B,Zhang Z F,Guo B X,Yu H J,Zong F R 2013Acta Phys.Sin.62 147602(in Chinese)[李新,肖立志,劉化冰,張宗富,郭葆鑫,于慧俊,宗芳榮2013物理學報62 147602]

[3]Behroozmand A A,Keating K,Auken E 2015Surv.Geophys.36 27

[4]Hertrich M 2008Prog.Nucl.Mag.Res.Sp.53 227

[5]Hürlimann M D 1998J.Mag.Res.131 232

[6]Chen Q,Marble A E,Colpitts B G,Balcom B J 2005J.Mag.Res.175 300

[7]Grunewald E,Knight R 2012Geophysics77 EN1

[8]Legchenko A,Vouillamoz J M,Lawson F M A,Alle C,Descloitres M,Boucher M 2016Geophysics81 WB23

[9]Legchenko A,Vouillamoz J M,Roy J 2010Geophysics75 L91

[10]Lin J,Jiang C D,Lin T T,Duan Q M,Wang Y J,Shang X L,Fan T H,Sun S Q,Tian B F,Zhao J,Qin S W 2013Chin.J.Geophys.56 3619(in Chinese)[林君, 蔣川東,林婷婷,段清明,王應吉,尚新磊,范鐵虎,孫淑琴,田寶鳳,趙靜,秦勝武2013地球物理學報56 3619]

[11]Walbrecker J O,Hertrich M,Green A G 2011Geophysics76 G1

[12]Grombacher D,Walbrecker J O,Knight R 2014Geophysics79 E329

[13]Irons T P,Li Y 2014Geophys.J.Int.199 1372

[14]Grombacher D,Knight R 2015Geophysics80 E329

[15]Grombacher D,Müller-Petke M,Knight R 2016Geophysics81 WB33

[16]Mueller-Petke M,Yaramanci U 2010Geophysics75 WA199

[17]Roy J,Lubczynski M W 2014Near Surf.Geophys.12 309

[18]Chen B,Hu X,Li J,Liu Y 2016Groundwater55 171

[19]Weichman P B,Lavely E M,Ritzwoller M H 2000Phys.Rev.E62 1290

[20]Lehmann-Horn J A,Hertrich M,Greenhalgh S A,Green A G 2011IEEE Trans.Geosci.Remote.Sens.49 3878

[21]Bloch F 1946Phys.Rev.70 460

[22]Jiang C D,Lin J,Duan Q M,Tian B F,Hao H C 2011Chin.J.Geophys.54 2973(in Chinese)[蔣川東,林君,段清明,田寶鳳,郝薈萃2011地球物理學報54 2973]

[23]Walbrecker J O,Hertrich M,Green A G 2009Geophysics74 G27

[24]Günther T,Rücker C,Spitzer K 2006Geophys.J.Int.166 506

[25]Jiang C,Liu J,Tian B,Sun S,Lin J,Müller-Petke M 2016Geophysics81 E363