基于動態輪廓模型的地面核磁共振反演方法技術

本發明專利技術公開一種基于動態輪廓模型的地面核磁共振反演方法，其在SNMR層狀反演求解過程中，首先設定探測區域內的含水層數，并初始化各層的厚度和含水量值；然后，將其劃分為MN個微元，使其滿足動態模型的矩陣方程；最后，對矩陣方程迭代求解，迭代過程中對各個含水層的厚度和含水量值進行動態調整，以搜索滿足矩陣方程的最優解。整個求解過程，對探測區域內含水量值的垂直分布圖的輪廓不斷地進行動態調整，故該方法稱為“動態輪廓模型”，并且設計了隨機梯度下降法(SGD)求解該模型；本發明專利技術具有收斂速度快、反演結果精度高而且穩定，其性能優于正則化方法，能夠實現了對SNMR反演問題的高精度求解。

【技術實現步驟摘要】
基于動態輪廓模型的地面核磁共振反演方法
本專利技術涉及地面核磁共振(SurfaceNuclearMagneticResonance，簡稱SNMR)領域，具體涉及一種基于動態輪廓模型的地面核磁共振反演方法。
技術介紹
地面核磁共振(SurfaceNuclearMagneticResonance，簡稱SNMR)技術是目前世界上唯一的一種直接找水的物探方法，該項技術已在探測地下水、考古、地下水污染檢測等領域得到了一定的應用。近年來，隨著專家和學者們的逐漸深入研究，SNMR技術得到了進一步的完善。反演計算含水率是該技術研究過程中的關鍵環節，而反演結果的準確度和反演過程的運算效率是衡量反演算法性能的關鍵指標。其中，一維正反演理論較為成熟，已經相繼刊登出多種有效算法，如：改進的模擬退火算法反演，提高了現有反演算法的穩定度和收斂速度；QT反演算法，利用各個激發脈沖矩對應的全部采樣點數據進行反演，充分挖掘了接收信號信息，在一定程度上提高了反演精度，但是，由于接收信號呈現近似指數衰減，晚期信號信噪比很低，該方法只適用于高信噪比環境；后來也有學者采用了積分門技術接收信號，提高了各個采樣點數據的精度，并進行全衰減反演，是對QT反演的一種改進。SNMR技術在二維反演方面，Boucher、Girard和Legchenko等研究了在二維剖面方向上E0-q曲線隨地下含水構造的變化趨勢，但他們只對二維反演做了定性研究，沒有給出具體的二維反演公式。Legchenko等對三維反演做了一定的研究，雖然能在三維空間反演出模型的含水構造，但是由于在三維空間設定的網格尺寸較大，只能粗略的估計出地下含水構造，其反演分辨率有待提高。由于二維、三維反演算法存在運算量大、待求解變量數多、非線性等問題，目前世界上唯一商業版反演軟件NUMISPLUS仍采用一維反演。地面核磁共振反演能抽象為一個求解病態矩陣方程的問題，現有的文獻中多采用正則化方法對其進行求解，但是，現有的標準正則化方法以及在其基礎上的改進方法存在兩方面的不足，分別是：正則化參數選取困難；將反演問題的近似解強制添加光滑性約束，弱化了相鄰地質層之間含水量的差異。
技術實現思路
本專利技術所要解決的技術問題是提供一種基于動態輪廓模型的地面核磁共振反演方法，其能夠實現對地面核磁共振反演問題的高精度求解。為解決上述問題，本專利技術是通過以下技術方案實現的：一種基于動態輪廓模型的地面核磁共振反演方法，包括如下步驟：步驟1，建立的地面核磁共振反演的動態輪廓模型，即An＝E①式中，A為核函數矩陣；n為待求解的含水量向量，MN為模型中含水層的微元數；E為地面核磁共振信號初始振幅值向量，IN為激發脈沖矩個數；步驟2，初始化探測區域內的含水層數為L，各層的厚度和含水量值分別為HL＝(h1,h2,...,hL)T和NL＝(n1,n2,...,nL)T；最大迭代次數為Nmax；初始迭代次數為NI＝0；反演精度閾值設定為θ；步驟3，將當前一次迭代的各個子含水量垂直分布圖劃分為MN個微元，每個微元厚度為△h，并生成動態輪廓模型的核函數矩陣A；步驟4，計算動態輪廓模型的下一代搜索路徑，并將更新的搜索路徑視為當前一次迭代的搜索路徑；其中下一代搜索路徑更新公式為20nk+1*=nk+wnΔvknk+1=0,nk+1*<0nk+1*,0≤nk+1*≤11,nk+1*>1]]>②式中，hk和nk為本次迭代各層的厚度和含水量，hk+1和nk+1為下一次迭代各層的厚度和含水量；wh和wn分別為各層的厚度和含水量值的移動權重；△uk和△vk分別為各層的厚度和含水量值的移動方向向量；上述△uk為各層厚度移動方向矩陣△u中的第k列，即△vk為各層含水量移動方向矩陣△v中的第k列，其中式中，各層厚度移動方向矩陣△u中的元素△uij∈{x|-1,0,1}，各層含水量移動方向矩陣△v中的元素△vij∈{x|-1,0,1}，且每次迭代△u和△v均隨機生成。步驟5，計算動態輪廓模型的適應度函數，更新最優解；適應度函數定義為：fmin＝||Ank+1-E0||2⑤式中，A為核函數矩陣，nk+1為當前一代一個子含水量向量，E0為實測地面核磁共振信號中各個激發脈沖矩對應的初始振幅值；計算當前一代中各個子含水量向量中的最優解；當前一代最優解中若含有層厚度為△h米的含水層時，則刪除該層，并將其上一層的層厚度增加△h米，更新探測區域內的含水層數L＝L-1；若是第一層層厚度為△h米時，則將其與第二層合并；步驟6：判斷迭代是否停止；如果迭代次數達到最大迭代次數Nmax，或者當前最優適應度函數值小于反演精度閾值θ，則停止迭代，反演結果為當前最優解；否則，返回步驟4，進行下一次迭代。所述步驟4中，各層的厚度和含水量值的移動權重wh和wn采用變步長搜索，即各層的厚度的移動權重wh為wh＝(1-NI/Nmax)△wh⑥式中，Nmax為最大迭代次數，NI為當前迭代次數，△wh為厚度的固定步長。各層的含水量的移動權重wn為wn＝(1-NI/Nmax)△wn⑦式中，Nmax為最大迭代次數，NI為當前迭代次數，△wn為含水量的固定步長。本專利技術在SNMR層狀反演求解過程中，首先設定探測區域內的含水層數，并初始化各層的厚度和含水量值；然后，將其劃分為MN個微元，使其滿足動態模型的矩陣方程；最后，對矩陣方程迭代求解，迭代過程中對各個含水層的厚度和含水量值進行動態調整，以搜索滿足矩陣方程的最優解。整個求解過程，對探測區域內含水量值的垂直分布圖的輪廓不斷地進行動態調整，故該方法稱為“動態輪廓模型”，并且設計了隨機梯度下降法(SGD)求解該模型；本專利技術具有收斂速度快、反演結果精度高而且穩定，其性能優于正則化方法，能夠實現了對SNMR反演問題的高精度求解。附圖說明圖1為基于動態輪廓模型的地面核磁共振反演方法與基于正則化模型的地面核磁共振反演方法在求解過程中擬合誤差收斂曲線圖。圖2為基于動態輪廓模型的地面核磁共振反演方法中初始振幅值E0和激發脈沖矩q的實際值、測量值和反演結果值的對比曲線。圖3為基于動態輪廓模型的地面核磁共振反演方法中實際模型和反演結果的含水量垂直分布對照圖。具體實施方式為了驗證本專利技術動態輪廓模型用于SNMR反演的有效性，首先建立含水層垂直分布模型，然后通過正演模型計算一系列激發脈沖矩對應的初始振幅的理論值，在初始振幅理論值的基礎上加入一定信噪比的噪聲作為觀測信號，最后分別采用ACM模型(AdvancedComputerModeling，先進的計算機模型)和TRM模型(theTikhonovregularizationmethod,正則化模型)進行反演。同時，為了實現模型對比，采用隨機梯度下降法分別求取上述兩個模型的最優解。定義信噪比SNR公式為：...

【技術保護點】
基于動態輪廓模型的地面核磁共振反演方法，其特征是包括如下步驟：步驟1，建立的地面核磁共振反演的動態輪廓模型，即An＝E????①式中，A為核函數矩陣；n為待求解的含水量向量，MN為模型中含水層的微元數；E為地面核磁共振信號初始振幅值向量，E=(E(q1),E(q2),...,E(qIN))T,]]>IN為激發脈沖矩個數；步驟2，初始化探測區域內的含水層數為L，各層的厚度和含水量值分別為HL＝(h1,h2,...,hL)T和NL＝(n1,n2,...,nL)T；最大迭代次數為Nmax；初始迭代次數為NI＝0；反演精度閾值設定為θ；步驟3，將當前一次迭代的各個子含水量垂直分布圖劃分為MN個微元，每個微元厚度為△h，并生成動態輪廓模型的核函數矩陣A；步驟4，計算動態輪廓模型的下一代搜索路徑，并將更新的搜索路徑視為當前一次迭代的搜索路徑；其中下一代搜索路徑更新公式為hk+1*=hk+whΔukhk+1=1,hk+1*<1hk+1*,1≤hk+1*≤2020,hk+1*>20nk+1*=nk+wnΔvknk+1=0,nk+1*<0nk+1*,0&le;nk+1*≤11,nk+1*>1]]>????②式中，hk和nk為本次迭代各層的厚度和含水量，hk+1和nk+1為下一次迭代各層的厚度和含水量；wh和wn分別為各層的厚度和含水量值的移動權重；△uk和△vk分別為各層的厚度和含水量值的移動方向向量；上述△uk為各層厚度移動方向矩陣△u中的第k列，即△vk為各層含水量移動方向矩陣△v中的第k列，Δvk=(Δv1k,Δv2k,...,ΔvMNk)T;]]>其中式中，各層厚度移動方向矩陣△u中的元素△uij∈{x|?1,0,1}，各層含水量移動方向矩陣△v中的元素△vij∈{x|?1,0,1}，且每次迭代△u和△v均隨機生成。步驟5，計算動態輪廓模型的適應度函數，更新最優解；適應度函數定義為：fmin＝||Ank+1?E0||2????⑤式中，A為核函數矩陣，nk+1為當前一代一個子含水量向量，E0為實測地面核磁共振信號中各個激發脈沖矩對應的初始振幅值；計算當前一代中各個子含水量向量中的最優解；當前一代最優解中若含有層厚度為△h米的含水層時，則刪除該層，并將其上一層的層厚度增加△h米，更新探測區域內的含水層數L＝L?1；若是第一層層厚度為△h米時，則將其與第二層合并；步驟6：判斷迭代是否停止；如果迭代次數達到最大迭代次數Nmax，或者當前最優適應度函數值小于反演精度閾值θ，則停止迭代，反演結果為當前最優解；否則，返回步驟4，進行下一次迭代。...

【技術特征摘要】
1.基于動態輪廓模型的地面核磁共振反演方法，其特征是包括如下步驟：步驟1，建立地面核磁共振反演的動態輪廓模型，即An＝E①式中，A為核函數矩陣；n為待求解的含水量向量；E為地面核磁共振信號初始振幅值向量；步驟2，初始化探測區域內的含水層數為L，各層的厚度和含水量值分別為HL＝(h1,h2,...,hL)T和NL＝(n1,n2,...,nL)T；最大迭代次數為δmax；初始迭代次數為0；反演精度閾值設定為θ；步驟3，將當前一次迭代的各個子含水量垂直分布圖劃分為MN個微元，每個微元厚度為△h，并生成動態輪廓模型的核函數矩陣A；步驟4，計算動態輪廓模型的下一代搜索路徑，并將更新的搜索路徑視為當前一次迭代的搜索路徑；其中下一代搜索路徑更新公式為式中，hk和nk為本次迭代各層的厚度和含水量，hk+1和nk+1為下一次迭代各層的厚度和含水量；wh和wn分別為各層的厚度和含水量值的移動權重；△uk和△vk分別為各層的厚度和含水量值的移動方向向量；上述△uk為各層厚度移動方向矩陣△u中的第k列，即△vk為各層含水量移動方向矩陣△v中的第k列，其中式中，各層厚度移動方向矩陣△u中的元素△uij∈{x|-1,0,1}，各層含水量移動方向矩陣△v中的元素△vij∈{x|-1,0,1}，且每次迭代△u和△v均隨機生成...

【專利技術屬性】
技術研發人員：王國富，葉金才，張法全，張海如，韋秦明，龐成，王猛，
申請(專利權)人：桂林電子科技大學，
類型：發明
國別省市：廣西;45