楊玉晶 葉瑞 趙汗青 萬玲2) 林婷婷2)?

1) (吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長春 130061)

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

1 引 言

過去人們探測地下水主要采用的是電法和電磁法[1-4]， 通過探測含水構(gòu)造或者含水結(jié)構(gòu)體， 從而評價地下水的賦存狀態(tài)， 然而這些方法都是間接方法， 無法直接對地下水進行定性、定量估計， 也難以直接獲得賦水層相關(guān)參數(shù).相比于傳統(tǒng)電磁法，地面磁共振(magnetic resonance sounding， MRS)技術(shù)以其無損、原位， 且獨有的直接探測優(yōu)勢， 被廣泛應(yīng)用于淺層地下水勘探、水源性地質(zhì)災(zāi)害預(yù)測等領(lǐng)域[5].通常， 在實施探測時， 僅需施加地磁場( B0)條件下的拉莫爾頻率( f0)交流場， 激發(fā)自由水中氫原子核進動， 即可探知飽和地下水存在、含量及其賦存狀態(tài).激發(fā)結(jié)束后， 收集氫原子核能釋產(chǎn)生的自由感應(yīng)衰減(free induction decay， FID)信號， 其初始振幅、平均弛豫時間分別與地下含水量w、飽和水孔隙分布相關(guān)[6-8].近年來， 地面MRS 相關(guān)理論與方法日趨成熟， 其正反演技術(shù)已經(jīng)由簡單的一維層狀水探測成像逐漸發(fā)展為更為復(fù)雜的二維、三維成像[9-13].然而， 受單脈沖激發(fā)序列限制， 現(xiàn)有技術(shù)仍局限在w 及探測及解釋，在某些特定探測環(huán)境下， FID 信號易受地下鐵磁性介質(zhì)影響而大幅縮短(衰減快)， 造成含水結(jié)構(gòu)賦存狀態(tài)估測誤差極大[14-16].

相比于傳統(tǒng)的FID 探測， 基于自旋回波(spin echo， SE)信號的橫向弛豫時間 T2探測適應(yīng)性更強， 不僅能夠更為準(zhǔn)確地估計地下水孔隙分布、滲透率[17-20]， 且在水力傳導(dǎo)率等更多水文地質(zhì)信息預(yù)測方面表現(xiàn)出一定潛力[21，22]， 是目前地面MRS 領(lǐng)域的國際熱點問題.與測井技術(shù)中應(yīng)用的自旋回波探測序列[23]相似， 地面磁共振 T2探測通常由多組雙脈沖或多脈沖CPMG (Carr-Purcell-Meiboom-Gill)序列[24，25]激發(fā)， 通過擬合自旋回波峰值并代入反演， 得到地下水w 和 T2分布.隨著相關(guān)探測需求的提出， 兼具雙脈沖及CPMG 序列發(fā)射功能的地面MRS 儀器應(yīng)運而生[26].但由于早期工作更多出于檢測短弛豫信號的基本需要， 相關(guān)數(shù)據(jù)建模及解釋方法并無明顯發(fā)展， 仍采用與普通FID 信號相似策略， 即將SE 信號包絡(luò)峰值擬合線作為原始信號， 直接代入反演求解， 獲取地層含水信息[17].

為了更好地預(yù)測地下水賦存狀態(tài)， 充分利用SE 信號攜帶信息， 2014 年， Grunewald 等[19]提出了基于線性空間變換-非線性擬合的 T2反演框架.即首先通過奇異值分解(singular value decomposition， SVD)法[27]將時域自旋回波信號與脈沖矩間的函數(shù)關(guān)系變換至與層深相關(guān)， 得到地下水SE 分布; 隨后通過非線性擬合地下水SE 結(jié)果， 求得w 和 T2信息.該方法的提出， 直接建立了w 和T2與深度間的空間關(guān)系， 有效避免了原信號峰值取樣-反演擬合流程產(chǎn)生的二次誤差.然而， 雖然SVD 方法計算步驟簡單、操作方便， 但易受數(shù)據(jù)摻雜噪聲的影響， 難以避免病態(tài)化矩陣問題.雖然施加閾值約束等策略可有效降低分解矩陣的病態(tài)化程度[28]， 但大部分情況下， 其解矩陣依舊受噪聲干擾嚴(yán)重， 且存在部分不具備水文地質(zhì)意義的畸變元素， 嚴(yán)重影響后續(xù)w 及 T2擬合精度.

針對以上問題， 本文首先介紹了地面MRS 自旋回波正演建模過程， 推導(dǎo)了其靈敏度核函數(shù)及正演公式.再基于Grunewald 等[19]提出的 T2反演框架， 總結(jié)了地面SE 反演策略.針對SVD 法病態(tài)化嚴(yán)重、精度低等問題， 提出將高效、自帶非負(fù)約束的同時迭代重建技術(shù)(simultaneous iterative reconstruction technology， SIRT)作為SVD 的后續(xù)處理步驟[29]， 提高地面MRS 數(shù)據(jù)解釋精度.本文研究成果將進一步推動SE 探測在地面MRS 領(lǐng)域的應(yīng)用， 并促進相關(guān)正演建模及反演解釋方法的進步.

2 地面磁共振橫向弛豫探測原理

2.1 傳統(tǒng)地面磁共振探測原理

地面MRS 探測利用地下水中氫原子核的自旋現(xiàn)象， 其旋進角頻率(拉莫爾頻率)由 B0決定[9]:

其中， ω0=2πf0， γ 為地下水中氫原子核的磁旋比.所有原子核的自旋行為形成微量的宏觀磁矩M，其幅度表示為 M0.為探測地下水存在， 應(yīng)用布設(shè)于地表的線圈發(fā)射激發(fā)脈沖， 氫原子核吸收脈沖能量而進入激發(fā)態(tài)， M 的旋轉(zhuǎn)軸逐漸偏離原方向.此時， M 垂直于 B0方向分量被稱為激發(fā)橫向磁化強度:

激發(fā)場停止后， M 以角頻率 ω0旋進并回到激發(fā)前狀態(tài)(平衡狀態(tài))， 釋放對應(yīng)頻率交流場， 即為MRS 信號.在該過程中， 通常用 T1(縱向弛豫時間)、 T2等參數(shù)描述氫原子核的恢復(fù)平衡態(tài)過程(弛豫過程).T1與 T2分別為M 平行(縱向)、垂直(橫向)于 B0分量， 其恢復(fù)初始平衡狀態(tài)時間， 通常被認(rèn)為與地下水的孔隙度、滲透率有關(guān).

2.2 地面磁共振橫向弛豫探測發(fā)展

其中， Δ B0為地磁場不均勻程度.在均勻地磁場條件下(地下不存在鐵磁性介質(zhì))， (3)式第二項可忽略不計，可近似替代 T2.但當(dāng)?shù)卮艌霾痪鶆蚯闆r下， (3)式第二項的存在會導(dǎo)致遠小于 T2.此時直接應(yīng)用探測值代替 T2， 不僅使探測難度進一步增加(受儀器死區(qū)時間影響， 信號越短接收越困難)， 而且估測的地下孔隙半徑將遠小于真實情況.

近十年來， MRS 領(lǐng)域高效充放電及多脈沖發(fā)射技術(shù)逐漸發(fā)展成熟， 為基于SE 信號的 T2探測提供了條件.通常， 地面SE 信號探測有兩種方式: 多組不同時間間隔的雙脈沖探測序列組合和多脈沖連續(xù)發(fā)射的CPMG 序列.對于前者， 要求針對同一激發(fā)脈沖， 在其關(guān)斷后設(shè)置不同時間間隔的重聚脈沖.接收并擬合多個SE 信號包絡(luò)峰值， 即可通過數(shù)據(jù)解釋獲得地下 T2分布.由于每次測量均需要對同一激發(fā)脈沖矩求取其不同時間間隔下的SE信號， 所以該種探測方式效率較低(探測時間是傳統(tǒng)FID 實驗的數(shù)倍).相比于雙脈沖組合序列，CPMG 序列(圖1)由單個激發(fā)脈沖與相同時間間隔的多個重聚脈沖組成， 通過單次間斷發(fā)射， 即可獲得多重SE 信號.其過程為: t =0 時發(fā)射激發(fā)脈沖q， 激發(fā)地下水中氫原子核產(chǎn)生磁共振現(xiàn)象， 形成衰減時間為的FID 信號; t =τ，3τ，5τ，··· 時刻分別發(fā)射相同大小的重聚脈沖2q (大小為激發(fā)脈沖2 倍)， 即可接收在 t =2τ，4τ，6τ，··· 時刻達到峰值的多重SE 信號.針對不同探測層深， 分別發(fā)射對應(yīng)該深度的CPMG 序列(激發(fā)脈沖不同)， 最終完成整個區(qū)域內(nèi)的地下水探測.相比于雙脈沖激發(fā)方式， CPMG序列極大提高了探測效率， 進一步推動了相關(guān)理論及解釋算法的發(fā)展.

圖1 地面磁共振CPMG 序列激發(fā)示意圖(為獲取包含T2 衰減信息的完整SE 信號， 發(fā)射激發(fā)脈沖q， 并間隔時間τ 重復(fù)發(fā)射多個重聚脈沖2q， 重聚脈沖間的間隔為2τ)Fig.1.Typical time diagram of CPMG apply in surface MRS.Measurement for excitation pulse moment q requires a set of refocusing pulse 2q repeating for a defined 2τ delay after the first τ delay and obtains a complete SE response to determine the T2 decay curve.

需要注意的是， 受地下水?dāng)U散現(xiàn)象影響， 通過MRS 自旋回波探測得到的橫向弛豫時間也不是完全等同于真實橫向弛豫時間 T2.具體來說， 將前者定義為 T2MRS， 則[25]:

其中， D 為孔隙內(nèi)填充流體(地下水)的擴散系數(shù)， G為空間磁場(地磁場)梯度.為盡量縮小忽略(4)式第二項對探測的影響， 宜選取較小的 τ 值， 以保證探測到的 T2MRS更接近 T2.在本文中， 由于該關(guān)系與主體研究內(nèi)容無關(guān)， 故默認(rèn) T2MRS=T2， 且下文統(tǒng)一用 T2代替 T2MRS表述地面MRS 探測得到的橫向弛豫時間.

3 地面磁共振橫向弛豫正反演理論

3.1 靈敏度核函數(shù)及正演計算

如前文所述， 氫原子核受激發(fā)后章動能夠釋放交流弛豫信號.對應(yīng)于不同激發(fā)脈沖q， 其MRS 信號響應(yīng)可描述為[9]

式中， w (r) 為對應(yīng)于地下位置r 處的含水量;K(q，r)為MRS 的正演靈敏度核函數(shù)，

與FID 建模方式不同， SE 信號的橫向磁化強度計算還需充分考慮多個重聚度脈沖帶來的影響.參考磁共振測井的直接建模理論[32]， 地面MRS 自旋回波的激發(fā)橫向磁化強度可以表示為[17]

其中， θ′為重聚脈沖形成的扳倒角.考慮到通常CPMG 的重聚脈沖設(shè)定為對應(yīng)激發(fā)脈沖的兩倍，則SE 信號扳倒角 θ′=2θ ， (7)式可簡化為

為對比FID 與SE 探測區(qū)別， 驗證重聚脈沖對激發(fā)過程的影響， 計算了同一配置下， 兩種激發(fā)方式的橫向磁化強度分布(發(fā)射、接收線圈均為邊長100 m 方形、單匝)， 結(jié)果見圖2.由于采用同一激發(fā)脈沖(2 A·s)， 二者的橫向磁化強度分布趨勢相似， 但在效率上， 由于SE 信號接收發(fā)生在FID 信號之后， 其理論幅值明顯低于FID 信號.在上述結(jié)果基礎(chǔ)上， 仿真二者的一維靈敏度核函數(shù)， 即分別將橫向磁化強度(2)式、(8)式代入靈敏度核函數(shù)表達式(6)式， 并按層狀大地結(jié)構(gòu)積分， 如圖3(a)和圖3(b)所示.由圖可知， 對應(yīng)于橫向磁化強度分布， SE 信號的靈敏度核函數(shù)幅值也相對降低.

由于SE 信號形式較為復(fù)雜(包含多種不同狀態(tài)氫原子核自旋散相過程)， 難以通過固定公式精確定義其包絡(luò)， 通常僅定義其包絡(luò)峰值.根據(jù)(6)式計算SE 信號靈敏度核函數(shù)， 則層狀含水結(jié)構(gòu)(一維)下， t =0 時刻(即激發(fā)脈沖關(guān)斷瞬間)開始接收的第n 個自旋回波包絡(luò)峰值為

其中， z 對應(yīng)于一維大地層深， tn為第n 個SE 信號達到峰值的時刻.圖4 為典型含水結(jié)構(gòu)下的FID 與SE 磁共振響應(yīng)(地下8—20 m 存在單一含水層， 含水量最大30%， 對應(yīng)或 T2為400 ms).由于對應(yīng)SE 信號(圖4(b))的靈敏度核函數(shù)整體幅度偏低， 且隨著接收時間增大幅度逐漸衰減， 所以對應(yīng)于同一含水模型其整體幅度明顯低于FID信號(圖4(a))， 且伴隨每一次重聚脈沖發(fā)射， 呈現(xiàn)多峰值分布.值得注意的是， 圖4 仿真中認(rèn)為含水層= T2， 根據(jù)前文的(3)式， 實際情況下的應(yīng)略小于 T2， 但由于此處僅用于對比FID 與SE響應(yīng)形式、量級上的差別， 故并沒有嚴(yán)格區(qū)分二者間的大小.

為對SE 信號進行線性空間反演， 建立相應(yīng)正演形式如下:

圖2 2 A·s 激發(fā)脈沖下的地面磁共振(a) FID 與(b) SE 激發(fā)橫向磁化強度分布剖面(實驗配置為100 m 方形發(fā)射/接收線圈， 單匝; 填充顏色表示激發(fā)橫向磁化強度相對于總磁化強度比例)Fig.2.Excitation profile (2 A·s) for (a) FID and (b) SE responses in a surface MRS case with 100 m square transmitting/receiving loop， 1-turn configuration.Color indicates the amplitude ratio of the excited transverse magnetization to the total magnetization.

圖3 單匝100 m 方形發(fā)射/接收線圈探測配置下的地面磁共振(a) FID 與(b) SE 一維靈敏度核函數(shù)(填充顏色反映在對應(yīng)于激發(fā)脈沖矩， 某一固定深度下存在地下水能夠誘發(fā)得到的磁共振信號幅度)Fig.3.Comparison of kernel function with 100 m transmitting/receiving square loop， 1-turn configuration for (a) FID and (b) SE excitation.Color reflects the signal amplitude induced by underground water at a given depth layer corresponding to each pulse moment.

圖4 單匝100 m 方形發(fā)射/接收線圈探測配置下， 同一含水情況下的地面磁共振響應(yīng)(a) FID 與(b) SE 信號(填充顏色反映對應(yīng)激發(fā)脈沖矩， 其FID 或SE 信號隨時間衰減的磁共振信號幅度)Fig.4.Comparison of forward responses for the same aquifer distribution， with 100 m transmitting/receiving square loop， 1-turn configuration for (a) FID and (b) SE excitation.Color reflects the signal amplitude induced by underground water decays with receiving time corresponding to each pulse moment.

式中， m (z，t，tn) 為 對應(yīng)于 Ve(q，t，tn) 的自旋回波變換結(jié)果， 也是反演要求解的量， 此處稱為地下水SE信號， 包含了最終要求取的w 與 T2信息，

其中， s =1， 2， 3，···， g 為多種弛豫分量對應(yīng)序號，ws(z) 與 T2s(z) 分別代表單一巖性孔隙物質(zhì)局部含水量及橫向弛豫時間.

3.2 反演策略

將(10)式離散為

即為反演求解的矩陣形式.其中， Iq， Iz與 It分別為脈沖矩、大地層數(shù)及信號時間采樣點.由于K 與Ve均為復(fù)數(shù)域矩陣， 為使通過m 矩陣求取的w， T2信息具有實際意義， 通常對K 陣進行雅克比變換[30].通過SVD 方法， 將變換后的靈敏度核函數(shù)矩陣分解為

其 中， U 與D 為 正 交 矩 陣，S =diag[s1，s2，··· ，sI，0，··· ，0]， I =min(Iq，Iz)為對角矩陣， 包含了奇異值信息.則矩陣(12)式的解為

由于奇異值分解過程中， 奇異值矩陣S 存在不同程度的病態(tài)化， 所以通常在上述求解過程中，還需引入合理的截斷方案.本文中， 引入Müller-Petke 與Yaramanci 在地面MRS 分辨率計算中采用的濾波法[33]， 即應(yīng)用奇異值結(jié)合正則化參數(shù)構(gòu)建對角矩陣 F =diag[f1，f2，···，fI] ， 作為濾波系數(shù)矩陣， 對S 進行預(yù)處理， 則第i 個濾波因子 fi可表示為

其中， λ 為正則化參數(shù).則此時通過SVD 分解求得的矩陣(12)式解為

在正則化參數(shù)選擇過程中， 通常引入偏差度準(zhǔn)則等方案， 采用環(huán)境噪聲水平等因素判斷當(dāng)前正則化參數(shù) λ 是否合適.即便如此， 由于實測數(shù)據(jù)中環(huán)境噪聲難以被完全壓制， 殘留的噪聲將進一步加劇SVD 分解病態(tài)程度， 并直接影響反演精度.圖5顯示了相同含水模型(圖5(a)， 地下8—20 m 存在單一含水層， 含水量最大30%， 對應(yīng)T2為400 ms)，在不同殘留環(huán)境噪聲情況下的SVD 反演結(jié)果.可以得出結(jié)論， 當(dāng)原始數(shù)據(jù)中無噪聲殘留時(圖5(b))，SVD 方法的整體反演效果較好， 除末期數(shù)據(jù)由于SE 信號衰減幅度較低， 出現(xiàn)輕微誤差外， 基本還原了真實的模型.但當(dāng)原始數(shù)據(jù)中存在噪聲殘留(圖5(c)， 3 nV 高斯白噪聲)， 盡管噪聲殘留較小，依舊會明顯影響反演結(jié)果， 且隨著噪聲殘留增大，其反演結(jié)果受干擾情況更加明顯(圖5(d)， 6 nV高斯白噪聲).

圖5 (a)模型與其在(b)無噪聲， (c) 3 nV 及(d) 6 nV 高斯白噪聲情況下的SVD 線性空間反演結(jié)果Fig.5.(a) Simulated modeling and its linear spatial inversion results employing SVD with (b) no noise， (c) 3 nV and (d) 6 nV Gaussian noise.

由于SVD 反演精度有限， 為進一步提高求解質(zhì)量， 引入SIRT 方法作為SVD 的后續(xù)處理步驟.SIRT 作為代數(shù)重建技術(shù)[34]的改善算法， 被廣泛應(yīng)用于磁共振測井?dāng)?shù)據(jù)的 T1， T2譜解釋[29]， 其不僅操作簡單、自帶非負(fù)約束， 且具有優(yōu)越的收斂特性.但由于該方法對初始模型要求較高， 所以單獨使用時并不穩(wěn)定.而將SVD 分解得到的反演解， 作為初始模型 m(1)代入SIRT 運算， 即可為SIRT 運算提供穩(wěn)定的初始解， 又可優(yōu)化SVD 反演結(jié)果， 能夠進一步提升地面MRS 解釋精度， 其具體過程如下.

首先， 認(rèn)為由初始模型 m(1)得到的模擬數(shù)據(jù)V(1)與實測數(shù)據(jù)Vobs差值， 完全由 m(1)與真實地下自旋回波間誤差造成， 則根據(jù)(12)式可得

將上述矩陣按照時間采樣點t=tit(it=1， 2，··· ，It)分割成若干方程組:

為求解每層地下水自旋回波誤差 Δ miz，it， 根據(jù)靈敏度核函數(shù)確定每層的分配系數(shù):

并 將 信 號 誤差 Δ viq，it按照分配 系 數(shù) ρiq， 分配給每層地下水誤差 Δ miz，it， 則

通過多次迭代， 不斷更新SIRT 求解結(jié)果:

其中， j 是當(dāng)前的迭代次數(shù).隨著迭代進行， 擬合與實測數(shù)據(jù)誤差不斷縮小， 最終可求得地下水SE 分布矩陣m.

將前文中SVD 反演結(jié)果(圖5)代入SIRT 進一步迭代擬合， 如圖6 所示.可以看到， 對于無噪聲情況下的數(shù)據(jù)反演， SIRT 方法可有效彌補SVD的誤差， 其反演結(jié)果幾乎與原始模型一致(圖6(a)).在分別加入3 nV (圖6(b))與6 nV (圖6(c))高斯白噪聲情況下， 單一SVD 方法僅能反演得到前兩個地下水自旋回波， 其邊緣分辨率極低， 且非含水區(qū)域受噪聲殘留影響， 形成多處明顯的波動.對于反演得到的數(shù)據(jù)尾部， 已幾乎難以判定第4、第5自旋回波是噪聲殘留還是信號.應(yīng)用SIRT 方法進一步處理上述反演結(jié)果， 能夠明顯提升前兩個自旋回波的邊緣分辨率， 并識別最末的自旋回波信號，壓制非含水區(qū)域殘留噪聲引起的波動.

圖6 (a)模型與其在(b) 3 nV 和(c) 6 nV 高斯白噪聲情況下的SVD 與SIRT 聯(lián)合線性空間反演結(jié)果Fig.6.(a) Simulated modeling and its linear spatial inversion results employing SVD and SIRT with (b) 3 nV and (c) 6 nV Gaussian noise.

為驗證與SVD 聯(lián)合處理過程中， SIRT 策略的收斂特性， 基于圖6 數(shù)據(jù)計算過程， 同時比較了不同噪聲情況下的迭代次數(shù)(圖7).受噪聲波動影響， 含噪聲幅度越大， 反演得到的最終數(shù)據(jù)誤差也相應(yīng)增大， 但收斂所需的迭代次數(shù)基本維持不變(約100 次)， 這說明了本文提出的策略用于地面MRS 自旋回波反演時， 始終維持著快速收斂的特性.

圖7 不同噪聲情況下， SIRT 線性空間反演迭代次數(shù)與擬合誤差間的關(guān)系Fig.7.Relationship between iteration and fitting errors for SIRT linear spatial inversion with different noise cases.

為得到最終的w 與 T2分布， 單獨提取對應(yīng)于每層的地下水SE 信號包絡(luò)， 并結(jié)合其峰值信息，進行單指數(shù)或多指數(shù)擬合.為了降低原始數(shù)據(jù)中殘留噪聲及求解誤差對擬合的影響， 在上述擬合過程中， 通常還需考慮加入窄時間窗， 即分別以t=2τ，4τ，6τ，···時間點為中心， 取一定范圍內(nèi)包絡(luò)幅值平均量作為擬合峰值， 其過程如圖8 所示[19].

4 實驗驗證

為評價本文方法的有效性， 本節(jié)通過模擬實際野外的CPMG 探測實驗， 獲取對應(yīng)SE 信號響應(yīng)，并對該信號進行反演， 最終求得地下w 及 T2分布.

圖8 地面磁共振SE 信號線性空間反演及含水量-橫向弛豫時間分布擬合示意圖Fig.8.Schematic diagram of linear spatial inversion and non-linear fitting of water content and transverse relaxation time for surface MRS spin-echo responses.

4.1 仿真實驗設(shè)置

雖然SE 探測主要應(yīng)用在非均勻地磁場分布情況下， 其地下電阻率與拉莫爾頻率分布并不均勻.但本節(jié)的仿真實驗?zāi)康臑轵炞C本文反演方法有效性， 故可忽略以上非均勻性帶來的影響.設(shè)定地磁場強度為54721 nT， 拉莫爾頻率(即設(shè)定的激發(fā)頻率)為2330 Hz， 地磁傾角與地磁偏角分別為60°， 0°.探測線圈采用邊長為100 m 的方形線圈，發(fā)射與接收均為單匝， 探測空間為線圈正下方(地下)電阻率為100 Ω·m 的均勻半空間.由于SE 信號的探測激發(fā)效率略低于FID， 由此設(shè)定反演最大深度為60 m， 激發(fā)脈沖為[0.1 8]A·s 間按指數(shù)取樣的40 組交流脈沖， 重聚脈沖為其倍數(shù)， 即分布在[0.2 16]A·s 間， 對應(yīng)每個激發(fā)脈沖發(fā)射5 次重聚脈沖， 其間隔參數(shù) τ =83 ms.

結(jié)合選取的探測深度及激發(fā)脈沖矩， 設(shè)定分布在—5 m 至—15 m、—30 m 至—50 m 的兩個含水層，背景含水量與橫向弛豫時間分別設(shè)定為1%和0.02 s， 如圖9(a).根據(jù)該含水模型擬合得到的每層最大w 分別為30%和50%， 對應(yīng) T2為0.4 s 和0.5 s.考慮到選取的參數(shù) τ 及重聚脈沖個數(shù)， 設(shè)定模擬信號的采樣時間為0.8 s.由于野外MRS 探測數(shù)據(jù)信噪比通常較低， 即使經(jīng)過硬件及算法多重噪聲壓制策略， 依舊殘留部分噪聲.為給實際的野外勘探提供參考經(jīng)驗， 同時驗證本文反演方法的抗干擾特性， 向模擬數(shù)據(jù)中加入3 nV 的高斯白噪聲， 則根據(jù)上述含水模型計算得到的加噪MRS 信號包絡(luò)，如圖9(b)所示.由于本文為仿真實驗， 且其探究重點是SE 信號的建模過程及反演解釋， 故信號響應(yīng)圖中沒有展示激發(fā)脈沖所引起的FID 信號.但在實際野外的SE 探測實驗中， 不僅能夠在MRS數(shù)據(jù)初始段觀測到FID 信號， 且在后續(xù)SE 信號解釋前， 還需首先排除接收FID 信號對SE 包絡(luò)提取的干擾.

圖9 (a)地面磁共振CPMG 探測實驗?zāi)Ｐ团c(b)正演數(shù)據(jù)(加入3 nV 高斯白噪聲)， 假設(shè)地下存在兩個含水層， 分別分布在—5 m 至—15 m 及—30 m 至—50 m.Fig.9.(a) Simulated modeling and (b) dataset (adding 3 nV Gaussian noise) for CPMG sequence assuming a surface MRS experiment with two aquifers， which distributed from —5 m to —15 m and —30 m to —50 m， respectively.

4.2 反演結(jié)果

作為對比， 仿真數(shù)據(jù)首先由SVD 方法單獨處理.前文提及的吉洪諾夫正則化及濾波器截斷方案被應(yīng)用于該反演過程.為保證反演結(jié)果具有實際水文地質(zhì)意義， 將病態(tài)矩陣導(dǎo)致的m 矩陣負(fù)值元素，均用零替代.如圖10(a)所示， 經(jīng)過反演過程， SE響應(yīng)由時間-脈沖間的函數(shù)關(guān)系轉(zhuǎn)換為時間-深度間函數(shù)關(guān)系.在噪聲影響下， 地下水SE 信號成像質(zhì)量較差， 甚至在地下12 m 左右出現(xiàn)“斷層”， 盡管能夠通過該反演結(jié)果大致判斷兩含水層存在位置，但各個地下水SE 信號邊緣分辨率較低.尤其對于深層含水層， 其各個SE 的幅度遠低于真實情況，且?guī)缀鯚o法分辨含水區(qū)域與非含水區(qū)域的界限.此外， 在模型中并不存在地下水區(qū)域， 均分布著明顯噪點， 且在地下20 m 處， 出現(xiàn)2 m 左右的“假”含水層.

由于原模型并不存在多弛豫分量情況， 故僅通過單指數(shù)對圖10(a)中地下水SE 信號進行擬合(15 ms 時間窗)， 得到w 及 T2結(jié)果如圖10(b)以及圖10(c)所示.對比原始模型擬合結(jié)果， 僅通過SVD 反演能夠相對精準(zhǔn)地求得第一個含水層的含水量信息， 但由于淺層m 矩陣中大量噪點的存在，由其擬合得的 T2幾乎毫無規(guī)律， 第一個含水層對應(yīng)的 T2信息幾乎被淹沒.在地下20 m 處， 由于m 矩陣中“假”含水層的存在， 對應(yīng)含水量與 T2均呈現(xiàn)小范圍突變.對于第二個含水層， 雖然通過其w 擬合結(jié)果基本能夠識別該含水層位置， 但含量遠低于真實值， 對應(yīng) T2信息也并不準(zhǔn)確.在50 m 深度后， 其w 與 T2信息更呈現(xiàn)出完全相反的趨勢(含水量隨深度降低， T2隨深度升高)， 與真實情況偏差極大.

應(yīng)用本文提出的SVD 與SIRT 聯(lián)合策略處理上述數(shù)據(jù)， 在經(jīng)過98 次迭代后， 反演結(jié)果基本達到穩(wěn)定， 如圖11(a)所示.正如預(yù)期， 聯(lián)合策略對淺層含水層幾乎實現(xiàn)了完全還原， 在深部含水層成像上， 雖然其邊緣分辨率稍有降低， 但依舊明顯高于SVD 方法單獨反演的結(jié)果.采用與前文相同手段提取該m 矩陣的w 及 T2信息， 對應(yīng)結(jié)果如圖11(b)和圖11(c)所示.由于第一個含水層的m 矩陣成像質(zhì)量更高， 其w 及 T2信息也擬合較好， 幾乎與原始模型完全重合， 它們最大誤差僅為1.5%和0.02 s.對于深層含水層， 由于m 矩陣中該部分的邊緣分辨率略低， 所以對應(yīng)含水層的邊緣部分?jǐn)M合結(jié)果稍有振蕩， 但w 與 T2信息基本還原了真實情況， 對應(yīng)峰值處誤差僅為5%和0.08 s.

圖10 加入3 nV 高斯白噪聲的地面磁共振仿真數(shù)據(jù)SVD 線性空間反演結(jié)果 (a) 地下水SE 響應(yīng)隨深度及接收時間的變化;(b)與(c)分別為應(yīng)用15 ms 時間窗(單指數(shù))擬合得到的w 與T2 分布結(jié)果Fig.10.Linear spatial inversion results of SVD method for synthetic surface MRS experiment data with 3 nV Gaussian noise polluted: (a) SE responses of underground water separated as a function of depth z and decayed over time， while the amplitude scale to water content; (b) and (c) are the subsurface w and T2 distribution fitted (mono-exponential) from (a) with a time window of 15 ms.

圖11 加入3 nV 高斯白噪聲的地面磁共振仿真數(shù)據(jù)SVD 與SIRT 線性空間反演結(jié)果 (a) 地下水SE 響應(yīng)隨深度及接收時間的變化; (b)與(c)分別為應(yīng)用15 ms 時間窗(單指數(shù))擬合得到的w 與T2 分布結(jié)果Fig.11.Linear spatial inversion results of SVD and SIRT method for synthetic surface MRS experiment data with 3 nV Gaussian noise polluted: (a) SE responses of underground water separated as a function of depth z and decayed over time， while the amplitude scale to water content; (b) and (c) are the subsurface w and T2 distribution fitted (mono-exponential) from (a) with a time window of 15 ms.

5 結(jié) 論

本文從傳統(tǒng)地面MRS 探測理論出發(fā)， 系統(tǒng)性介紹并總結(jié)了SE 信號正演理論.同時針對SE 信號反演問題， 提出了SIRT 與SVD 算法結(jié)合的數(shù)據(jù)解釋策略.該方法不僅能夠有效提升解釋精度，而且能夠避免單一SVD 方法存在的病態(tài)化矩陣問題.為了驗證本文方法的實際意義， 我們應(yīng)用加入噪聲的模擬數(shù)據(jù)作為實測數(shù)據(jù)， 分別用SVD 方法及SVD 與SIRT 結(jié)合策略處理上述數(shù)據(jù)， 取得的反演結(jié)果能夠精確反映地下真實的含水情況.結(jié)合全文理論分析及模擬實驗結(jié)果， 最終得出以下結(jié)論:

1)在地磁場不均勻情況下(地下存在鐵磁性介質(zhì))， 簡單的FID 探測并不適用; 應(yīng)用多脈沖探測序列獲得SE 信號， 求解地下 T2分布， 能更準(zhǔn)確地判斷地下含水層賦水信息; 受探測原理限制，SE 信號激發(fā)效率及幅度略低于同種情況下的FID信號.

2)在SE 實測信號解釋過程中， 引入磁共振測井相關(guān)理論， 能夠有效解決其靈敏度核函數(shù)計算及正演建模問題; 結(jié)合線性空間反演-非線性擬合策略， 可獲取其地下w 及 T2分布.

3)雖然直接將SVD 作為線性空間反演策略，也能實現(xiàn)w 及 T2求解， 但其對殘留噪聲的實測數(shù)據(jù)適應(yīng)性較差; 通過SVD 與SIRT 結(jié)合， 不僅能夠壓制前者病態(tài)化問題帶來的誤差， 而且能夠有效降低信號殘留噪聲影響， 提升數(shù)據(jù)解釋精度.

本文的研究成果， 將促進地面磁共振雙脈沖及多脈沖儀器， 以及相關(guān)正演理論、 T1及 T2反演方法的完善.

6 展 望

本文的正演理論， 主要基于測井領(lǐng)域引入的經(jīng)驗?zāi)Ｐ?在實際探測中， 自旋回波的激發(fā)過程更為復(fù)雜， 目前的正演模型僅能算作近似計算.且正如前文提及， 考慮到自旋回波探測的實際應(yīng)用場合，探測參數(shù)的設(shè)定及反演還需充分結(jié)合地面探測原理， 深入研究脈沖設(shè)定、多序列發(fā)射時間間隔等問題.需重點探究的是如何在現(xiàn)有儀器能夠?qū)崿F(xiàn)的發(fā)射功率下， 盡可能地提高發(fā)射電流、縮短發(fā)射時間，以降低地磁場不均勻造成的模型誤差影響.此外，如何避免激發(fā)FID 對SE 信號包絡(luò)提取影響等問題也需結(jié)合具體探測環(huán)境、探測參數(shù)設(shè)置進行解決.綜上所述， 在未來的研究中， 我們將進一步考慮模型的優(yōu)化及野外現(xiàn)場的實用性問題， 為地面磁共振 T1和 T2探測提供方法理論及儀器支撐.