鬼佬大哥大
  • / 14
  • 下載費用:30 金幣  

一種地面核磁共振二維反演方法.pdf

摘要
申請專利號:

CN201410252243.4

申請日:

2014.06.09

公開號:

CN103984033A

公開日:

2014.08.13

當前法律狀態:

授權

有效性:

有權

法律詳情: 授權|||實質審查的生效IPC(主分類):G01V 3/14申請日:20140609|||公開
IPC分類號: G01V3/14; G01V3/38 主分類號: G01V3/14
申請人: 桂林電子科技大學
發明人: 王國富; 張法全; 葉金才; 張海如; 韋秦明; 龐成; 陳俊婷
地址: 541004 廣西壯族自治區桂林市金雞路1號
優先權:
專利代理機構: 桂林市持衡專利商標事務所有限公司 45107 代理人: 陳躍琳
PDF完整版下載: PDF下載
法律狀態
申請(專利)號:

CN201410252243.4

授權公告號:

||||||

法律狀態公告日:

2017.01.11|||2014.09.10|||2014.08.13

法律狀態類型:

授權|||實質審查的生效|||公開

摘要

本發明公開了一種地面核磁共振二維反演方法,其用拉直變換方法將二維正演模型進行降維處理,將其抽象為矩陣方程求解模型,并用最小二乘奇異值分解(LS-SVD)與改進的隨機梯度下降法(ISGD)相結合的方法進行反演求解,采用LS-SVD求取矩陣方程的粗略解,在該粗略解的基礎上,用ISGD求取其精細解。在不同信噪比的條件下,本發明的反演結果均與模型中含水構造分布相吻合,即使在信噪比為0dB時,其反演結果仍能分辨出地下水文地質構造,其反演得到的含水量值的方均根為8.26%,而此時LS-SVD和ISGD兩種方法的反演結果均無效,其方均根值分別為30.14%和15.35%。

權利要求書

權利要求書
1.  一種地面核磁共振二維反演方法,其特征是包括如下步驟:
步驟1、用拉直變換方法將二維正演模型進行降維處理,將其抽象為矩陣方程求解模型;
步驟2、采用最小二乘奇異值分解法求取步驟1所抽象出的矩陣方程求解模型的粗略解nLS-SVD;
步驟2.1、對步驟1所抽象出的矩陣方程求解模型中的核函數矩陣K做奇異值分解,以獲得核函數矩陣K的奇異值σ;
步驟2.2、根據步驟2.1所獲得的核函數矩陣K的奇異值σ計算核函數矩陣K的有效秩r*;
步驟2.3、根據步驟2.1所獲得的核函數矩陣K的奇異值σ和步驟2.3所計算出的核函數矩陣K的有效秩r*,求取步驟1所抽象出的矩陣方程求解模型的最小二乘解;
步驟3、采用隨機梯度下降法求取步驟1所抽象出的矩陣方程求解模型的精細解n#;
步驟3.1、采用吉洪諾正則化方法構建模型的適應度函數,并將該適應度函數的當前一次迭代的解初始化為nLS-SVD,初始迭代次數初始化為0;
步驟3.2、根據搜索路徑更新公式計算模型的適應度函數的下一次迭代的搜索路徑,并將更新的搜索路徑視為當前一次迭代的搜索路徑;
步驟3.3、將當前一次迭代的搜索路徑中的每一個個體nh+1作為當前值ncur帶入步驟3.1所構建的模型的適應度函數中,計算適應度函數的當前一次迭代的最優解n#;
步驟3.4、如果迭代次數達到設定的最大迭代次數Nmax或者當前最優適應度函數值小于設定的反演精度閾值φ,則停止迭代,反演結果為當前最優解n#;否則,返回步驟3.2。

2.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟1所抽象出的矩陣方程求解模型為:
  ①
式中,E為初始振幅,K為核函數矩陣,n為所求向量,M為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數。

3.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟2.1,采用式②對矩陣方程求解模型中的核函數矩陣K做奇異值分解
KM×L=UM×MΛVL×LH]]>  ②
式中,UM×M和是正交矩陣,為VL×L的復共軛轉置,M為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數,Λ為對角陣,
  ③
式中,Λ為對角陣,M為剖面方向總共發射的激發脈沖矩的個數,L為矩陣K的列向量的個數,r為核函數矩陣K奇異值的個數,σ為核函數矩陣K的奇異值。

4.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟2.2具體為:

ψ(s)=σ1+σ2+...+σsσ1+σ2+...+σr,s=1,2,...,r]]>  ④
則有效秩r*為第一個滿足ψ(s)>θ的s值;
式中,ψ(s)為求取有效秩r*的表達式,σ為核函數矩陣K的奇異值,r為核函數矩陣K奇異值的個數,s為核函數矩陣K的第s個奇異值,θ為設定值。

5.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟2.3具體為:
nLS-SVD=Σf=1r*ufTEσfvf]]>  ⑤
式中,nLS-SVD為用最小二乘奇異值分解求取的解,E為初始振幅,σ為核函數矩陣K的奇異值,uf和vf分別為UM×M和VL×L的第f個列向量(f=1,2,…r*),ufT為uf的轉置,M為剖面方向總共發射的激發脈沖矩的個數,r*為核函數矩陣K的有效秩。

6.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟3.1構建出的模型的適應度函數為:
  ⑥
式中,Η(ncur)為模型的適應度函數,K為核函數矩陣,ncur為當前值,E為初始振幅,R為正則化因子,函數為用于求解ncur的1次偏導數。

7.  根據權利要求1所述的一種地面核磁共振二維反演方法,其特征是步驟3.2中,搜索路徑更新公式為:
nh+1*=nh+wΔvhnh+1=0,nh+1*<0nh+1*,0nh+1*11,nh+1*>1]]>  ⑦
式中,nh+1為當前一次迭代的適應度函數中ncur的值,nh為上一次迭代的適應度函數中ncur的值,w為移動權重,△vh為移動方向向量,△vh為移動方向矩陣△v中的第h列即△vh=(△v1h,△v2h,...,△vLh)T,其中
  ⑧
式中,△v為移動方向矩陣,△v的子個體△vij∈{x|-1,0,1},且每次迭代△v均隨機生成,M為每一次迭代的子個體數,亦為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數。

說明書

說明書一種地面核磁共振二維反演方法
技術領域
本發明涉及地面核磁共振領域,具體涉及一種地面核磁共振二維反演方法。
背景技術
地面核磁共振(Surface Nuclear Magnetic Resonance,簡稱SNMR)技術是目前世界上唯一的一種直接找水的物探方法,該項技術已在探測地下水、考古、地下水污染檢測等領域得到了一定的應用。近年來,隨著專家和學者們的逐漸深入研究,SNMR技術得到了進一步的完善。反演計算含水率是該技術研究過程中的關鍵環節,而反演準確度和分辨率是衡量反演算法性能的關鍵指標。其中,一維正反演理論較為成熟,已經相繼刊登出多種有效算法,如:文獻1[DAI Miao,HU Xiangyun,WU Haibo,et al.“Inversion of surface nuclear magnetic resonance for groundwater exploration,”Chinese Journal of Geophysics,2009,52(5):1166-1173.]提出了改進的模擬退火算法反演,提高了現有反演算法的穩定度和收斂速度;文獻2[Mueller-Petke M.,Yaramanci U..QT inversion-comprehensive use of the complete surface NMR data set[J].Geophysics,2010,75:199–209.]提出了QT反演算法,利用各個激發脈沖矩對應的全部采樣點數據進行反演,充分挖掘了接收信號信息,在一定程度上提高了反演精度;文獻3[Ahmad A.B ehroozmand,Esben Auken,Gianluca Fiandaca,et al.Efficient full decay inversion of MRS data with a stretched-exponential approximation of the T2*distribution[J].Geophysical Journal International,2012,190:900–912.]采用了積分門技術接收信號,提高各個采樣點數據的精度,并進行全衰減反演,是對QT反演的一種改進。在二維反演方面,Boucher、Girard和Legchenko等研究了在二維剖面方向上E0-q曲線隨地下含水構造的變化趨勢,但他們只對二維反演做了定性研究,沒有給出具體的二維反演公式。Legchenko等對三維反演做了一定的研究,雖然能在三維空間反演出模型的含水構造,但是由于在三維空間設定的網格尺寸較大,只能粗略的估計出地下含水構造,其反演分辨率有待提高。由于二維、三維反演算法存在運算量大、待求解變量數多、非線性等問題,目前世界上唯一商業版反演軟件NUMISPLUS仍采用一維反演,而二維、三維正反演研究仍處于起步階段。
發明內容
本發明所要解決的技術問題是現有地面核磁共振技術的實用性不強,一 維反演算法橫向分辨率低的不足,提出一種地面核磁共振二維反演方法。
為解決上述問題,本發明是通過以下技術方案實現的:
一種地面核磁共振二維反演方法,包括如下步驟:
步驟1、用拉直變換方法將二維正演模型進行降維處理,將其抽象為矩陣方程求解模型;
步驟2、采用最小二乘奇異值分解法求取步驟1所抽象出的矩陣方程求解模型的粗略解nLS-SVD;
步驟2.1、對步驟1所抽象出的矩陣方程求解模型中的核函數矩陣K做奇異值分解,以獲得核函數矩陣K的奇異值σ;
步驟2.2、根據步驟2.1所獲得的核函數矩陣K的奇異值σ計算核函數矩陣K的有效秩r*;
步驟2.3、根據步驟2.1所獲得的核函數矩陣K的奇異值σ和步驟2.3所計算出的核函數矩陣K的有效秩r*,求取步驟1所抽象出的矩陣方程求解模型的最小二乘解;
步驟3、采用隨機梯度下降法求取步驟1所抽象出的矩陣方程求解模型的精細解n#;
步驟3.1、采用吉洪諾正則化方法構建模型的適應度函數,并將該適應度函數的當前一次迭代的解初始化為nLS-SVD,初始迭代次數初始化為0;
步驟3.2、根據搜索路徑更新公式計算模型的適應度函數的下一次迭代的搜索路徑,并將更新的搜索路徑視為當前一次迭代的搜索路徑;
步驟3.3、將當前一次迭代的搜索路徑中的每一個個體nh+1作為當前值ncur帶入步驟3.1所構建的模型的適應度函數中,計算適應度函數的當前一次迭代的最優解n#;
步驟3.4、如果迭代次數達到設定的最大迭代次數Nmax或者當前最優適應度函數值小于設定的反演精度閾值φ,則停止迭代,反演結果為當前最優解n#;否則,返回步驟3.2。
上述步驟1所抽象出的矩陣方程求解模型為:
  ①
式中,E為初始振幅,K為核函數矩陣,n為所求向量,M為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數。
上述步驟2.1,采用式②對矩陣方程求解模型中的核函數矩陣K做奇異值分解
KM×L=UM×MΛVL×LH]]>  ②
式中,UM×M和是正交矩陣,為VL×L的復共軛轉置,M為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數,Λ為對角陣,
  ③
式中,Λ為對角陣,M為剖面方向總共發射的激發脈沖矩的個數,L為矩陣K的列向量的個數,r為核函數矩陣K奇異值的個數,σ為核函數矩陣K的奇異值。
上述步驟2.2具體為:

ψ(s)=σ1+σ2+...+σsσ1+σ2+...+σr,s=1,2,...,r]]>  ④
則有效秩r*為第一個滿足ψ(s)>θ的s值;
式中,ψ(s)為求取有效秩r*的表達式,σ為核函數矩陣K的奇異值,r為核函數矩陣K奇異值的個數,s為核函數矩陣K的第s個奇異值,θ為設定值。
上述步驟2.3具體為:
nLS-SVD=Σf=1r*ufTEσfvf]]>  ⑤
式中,nLS-SVD為用最小二乘奇異值分解求取的解,E為初始振幅,σ為核函數矩陣K的奇異值,uf和vf分別為UM×M和VL×L的第f個列向量(f=1,2,…r*),ufT為uf的轉置,M為剖面方向總共發射的激發脈沖矩的個數,r*為核函數矩陣K的有效秩。
上述步驟3.1構建出的模型的適應度函數為:
  ⑥
式中,Η(ncur)為模型的適應度函數,K為核函數矩陣,ncur為當前值,E為初始振幅,R為正則化因子,函數為用于求解ncur的1次偏導數。
上述步驟3.2中,搜索路徑更新公式為:
nh+1*=nh+wΔvhnh+1=0,nh+1*<0nh+1*,0nh+1*11,nh+1*>1]]>  ⑦
式中,nh+1為當前一次迭代的適應度函數中ncur的值,nh為上一次迭代的適應度函數中ncur的值,w為移動權重,△vh為移動方向向量,△vh為移動方向矩陣△v中的第h列即△vh=(△v1h,△v2h,...,△vLh)T,其中
  ⑧
式中,△v為移動方向矩陣,△v的子個體△vij∈{x|-1,0,1},且每次迭代△v均隨機生成,M為每一次迭代的子個體數,亦為剖面方向總共發射的激發脈沖矩的個數,L為核函數矩陣K的列向量的個數。
與現有技術相比,本發明的有益效果是:
1)通過拉直變換將二維正演模型進行降維處理,將二維反演問題抽象為矩陣方程求解模型,可參考現有的一維反演算法對其求解,降低了反演求解問題的復雜度。
2)提出了將成熟的一維反演算法LS-SVD與ISGD相結合的二維反演算法。首先,采用LS-SVD求取矩陣方程的粗略解nLS-SVD;然后,以nLS-SVD作為模型的初始化值,以Tikhonov正則化方法構建模型的適應度函數,用ISGD求取矩陣方程的精細解。
3)對現有的隨機梯度下降法進行改進,采用了變步長搜索,這既能加快算法的收斂速度,又能保證計算精度。
4)在不同信噪比下性能均優于LS-SVD和ISGD,其反演結果均與模型中含水構造分布比較吻合,即使在信噪比為0dB時,其反演結果仍能分辨出地下水文地質構造,其反演得到的含水量值的方均根為8.26%,而此時LS-SVD和ISGD反演結果無效,其方均根值分別為30.14%和15.35%。
附圖說明
圖1為二維剖面含水量模型;
圖2為fSNR=15dB時,S-SVD反演(圖2a)、ISGD反演(圖2b)、本發明反演(圖2c)三種算法反演結果對比圖。
圖3為fSNR=10dB時,LS-SVD反演(圖3a)、ISGD反演(圖3b)、本發明反演(圖3c)三種算法反演結果對比圖。
圖4為fSNR=5dB時,LS-SVD反演(圖4a)、ISGD反演(圖4b)、本發明 反演(圖4c)三種算法反演結果對比圖。
圖5為fSNR=0dB時,LS-SVD反演(圖5a)、ISGD反演(圖5b)、本發明反演(圖5c)三種算法反演結果對比圖。
具體實施方式
本發明所設計的地面核磁共振二維反演方法,采用最小二乘奇異值分解與改進的隨機梯度下降法相結合的方法進行反演求解。
在收/發天線共圈模式地面核磁共振(簡稱SNMR)二維正反演研究中,含水量沿著深度和剖面方向變化,測點p處接收信號的初始振幅E0(p,q)公式為:
E0(p,q)=&Integral;-&Integral;-K2D(p,q;x,z)&CenterDot;n(p;x,z)dxdzK2D(p,q;x,z)=&Integral;-K3D(p,q;x,y,z)dy&PartialD;n(p;x,y,z)/&PartialD;y=0---(1)]]>
其中,q為激發脈沖矩;K2D(p,q;x,z)和K3D(p,q;x,y,z)分別為p點處二維和三維核函數;n(p;x,z)和n(p;x,y,z)分別為p點處二維和三維含水量分布函數。
將探測剖面區域二維含水量n(x,z)沿著深度和剖面方向做離散化處理為:
zk-1zzk,k=1,2,...,Nzxj-1xxj,j=1,2,...,Nx]]>
得到L=Nz×Nx個含水微元。在剖面方向上有Np個測點,每個測點處發射Nq個激發脈沖矩,則剖面方向上總共發射M=Np×Nq個激發脈沖矩。測點pi處的第qh個激發脈沖矩對應的初始振幅為:
E0(pi,qh)=&Integral;-&Integral;-K2D(pi,qh;x,z)&CenterDot;n(pi;x,z)dxdz=Σj=1NxΣkNzK2D(pi,qh;xj,zk)&CenterDot;nkjΔzkΔxj---(2)]]>
對二維矩陣和進行拉直變換,將其轉化為一維列向量形式,即:
K2D&RightArrow;Km=[Km1,Km2,...KmLn&RightArrow;n*=[n1,n2,...,nL]T]]>
則(2)式轉化為:
Em=E0(pi,qh)=Km·n*  (3)
因此,SNMR二維反演問題可以抽象為求解矩陣方程E=K·n*中的n*值問題,即:

SNMR二維反演是一個帶約束的高維欠定非線性方程求解問題,且核函數具有病態性,為了高精度快速求解方程(4),本發明提出了基于LS-SVD的隨機梯度下降法,首先用LS-SVD求取粗略解,之后用該解作為ISGD的初始解,以在該粗略解周圍搜索精細解。
步驟1:采用LS-SVD求取方程(4)的粗略解nLS-SVD。
對方程(4)中的核函數矩陣K做奇異值分解
KM×L=UM×MΛVL×LH---(5)]]>
其中,UM×M和是正交矩陣,Λ為對角陣,即:

且σ1≥σ2≥…≥σr。
計算K的有效秩r*,令
ψ(s)=σ1+σ2+...+σsσ1+σ2+...+σr,s=1,2,...,r---(6)]]>
則r*為第一個滿足ψ(s)>θ的s值,通常θ取接近于1的值,在本發明優選實施例中θ=0.95。
方程(4)的最小二乘解為:
nLS-SVD=Σf=1r*ufTEσfvf---(7)]]>
其中,uf和vf分別為UM×M和VL×L的第f個列向量(f=1,2,…r*)。
步驟2:采用ISGD求取方程(4)的精細解n#。
2.1)算法初始化。解向量初始化為n0=nLS-SVD;最大迭代次數為Nmax;初始迭代次數為Nnum=0;每一代子個體數為M;反演精度閾值為φ,在本發明 優選實施例中,φ的取值為10-8;采用Tikhonov正則化方法構建模型的適應度函數為:

其中,R為正則化因子,在本發明優選實施例中,R的取值為10-8;函數用于求解ncur的1次偏導數。
2.2)計算下一代搜索路徑。下一代搜索路徑更新公式為:
nh+1*=nh+wΔvhnh+1=0,nh+1*<0nh+1*,0nh+1*11,nh+1*>1---(9)]]>
其中,w為移動權重,△vh為移動方向向量。
移動方向矩陣為其中,△vij∈{x|-1,0,1},且每次迭代△v均隨機生成。將△v向量化為△v=(△v1,△v2,…,△vM),其中,△vh=(△v1h,△v2h,...,△vLh)T。
為了加快收斂速度,同時保證計算精度,本發明對現有的隨機梯度下降法進行改進,采用變步長搜索,即:w=(1-Nnum/Nmax)△w,其中,△w為固定步長,在本發明優選實施例中,△w的取值為0.05。
w=(1-Nnum/Nmax)△w,
2.3)計算適應度函數,更新最優解。將步驟2.2)中所得當前一代中每一個體nh+1分別帶入公式(8),計算當前一代最優解。
2.4)判斷迭代是否停止。如果迭代次數達到Nmax,或者當前最優適應度函數值小于φ,則停止迭代,反演結果為當前最優解n#;否則,返回步驟2.2),重復執行步驟2.2)-2.4)。
本發明提出了采用拉直變換方法將二維正演模型進行降維處理,將二維反演模型轉化為矩陣方程求解的問題;由于轉化后的矩陣方程具有待求解變量數多、非線性、病態等特性,為了提高對其的求解速度和精度,并提出了用最小二乘奇異值分解(the least-squares singular value decomposition,簡稱LS-SVD)與改進的隨機梯度下降法(the improved stochastic gradient descent,簡稱ISGD)相結合的方法進行反演求解,采用LS-SVD求取矩陣方程的粗略解,在該粗略解的基礎上,用ISGD求取其精細解。最后,通過仿真實驗在不同信噪比下對本文研究內容進行了驗證。
為了驗證本發明所提出的SNMR二維反演算法的有效性,首先建立二維水文地質模型;然后,對其進行正演計算,求出各個測點的E0-q曲線;之后,對E0加入一定信噪比的噪聲,得到含噪信號E;最后,分別用LS-SVD、ISGD和本文算法對上述含噪信號進行二維反演,并用反演得到的含水量值的方均根fRMS評價反演結果性能,fRMS公式為:
fRMS=1LΣi=1L(ni-ni#)2×100%---(10)]]>
其中,ni和分別為理論模型和反演結果中各地質微元的含水量值。信噪比fSNR公式為:
fSNR=10lg(Σi=1MEi2Σi=1M(Ei-E0i)2)---(11)]]>
其中,Ei為信號序列,E0i為噪聲序列。
為了保證實驗結果的實用性,仿真模型盡可能的模擬實測環境,模型在X-Z剖面的尺寸為600m×90m,剖分的微元大小為20m×5m,剖面內共有30×18個含水微元,剖面270m至380m之間垂深25m至40m的區域含水量為50%,其他區域含水量為5%,其二維模型如圖1所示。共鋪設21個測點,測點坐標分別為p1(100m,0)、p2(120m,0)、……、p21(500m,0)。由于地面核磁共振設備接收靈敏度為10nV,每個測點pi(x,0)處只能接收X-Z剖面上x-100m至x+100m之間垂深90m以內區域的含水量信息,區域以外的核函數微元Kml為0。每個測點發送20個激發脈沖矩,其最大值為10As。地磁場強度為44630nT,磁偏角為24°,收/發線圈邊長為100m,背景電阻率為100Ω·m。在不同信噪比下,三種算法反演結果對比圖如圖2-圖5所示,表1中列出了其對應的方均根值,其中,ISGD和本發明的迭代次數都是1000次。
表1不同信噪比下,三種算法反演結果對比表
fSNR/dBfRMS(LS-SVD)/%fRMS(ISGD)/%fRMS(the new algorithm)/%156.1318.585.19109.2418.715.85516.1717.287.06030.1415.358.26
從圖2-圖5和表1的反演結果對比中,可以看出LS-SVD反演在高信噪比下,其反演精度較高,但是隨著信噪比的降低其精度快速降低,在0dB時,已經不能從反演結果中分辨出目標含水構造;2D SNMR反演是一個高維欠定非線性方程求解問題,由于求解過程中待求解變量數較多,ISGD在不同信噪比下反演精度均較差,雖然其反演結果能分辨出目標含水構造,但是剖面圖 的邊界處和垂深70m以下區域存在虛假含水信息;本發明的反演算法在不同信噪比下反演精度均優于LS-SVD和ISGD,算法穩定性好,抗噪能力強,隨著信噪比的降低反演結果的方均根值略有增加,在fSNR=0dB時,其反演精度仍然很高。

關 鍵 詞:
一種 地面 核磁共振 二維 反演 方法
  專利查詢網所有資源均是用戶自行上傳分享,僅供網友學習交流,未經上傳用戶書面授權,請勿作他用。
關于本文
本文標題:一種地面核磁共振二維反演方法.pdf
鏈接地址:http://www.wwszu.club/p-6140517.html
關于我們 - 網站聲明 - 網站地圖 - 資源地圖 - 友情鏈接 - 網站客服客服 - 聯系我們

[email protected] 2017-2018 zhuanlichaxun.net網站版權所有
經營許可證編號:粵ICP備17046363號-1 
 


收起
展開
鬼佬大哥大