一種基于APGD方法的隨鉆核磁共振測井快速反演算法

李沐堯，朱萬里，滕 朔，程晶晶*，葸春平

(1.華中科技大學 人工智能與自動化學院，湖北 武漢 430074;2.中國石油集團測井有限公司 測井技術研究院，陜西 西安 710077)

隨鉆核磁測井技術是一種非放射性的孔隙度測量技術，可以檢測泥漿濾液侵入前或者侵入很淺時的底層特征。與傳統電纜核磁測井相比，其可以更真實地反映儲層特征，評估精度更高[1]。受水平井或大斜度井限制，隨鉆測井系統無法采用電纜進行數據通信，僅能選取傳輸速率為4～10 bit/s的泥漿脈沖遙測技術[2]。在低傳輸速率下，無法直接上傳大量的原始測量數據，因此需要在井下的嵌入式設備中完成T2譜反演算法，上傳數據量較少的反演結果。

在反演算法方面，傳統算法包括奇異值分解法[3]、正則化方法[4]、非負最小二乘法[5]等。近些年開始對傳統算法進行改進。2009年，Prange等[6]提出了一種高效的Monte Carlo反演算法，并檢查光譜解決方案的統計特性。2015年，Martakusumah等[7]分別使用LM方法和Occam方法進行反演算法設計，對結果進行比較。2018年，岳鑫等[8]提出奇異值分解法的改進算法。趙宏晨等[9]對于Tikhonov正則化方法進行研究，建立了不適定問題方程組。2019年，Song等[10]提出正則化和約束矩陣擬合的反演方法。目前實際應用較多的方法為奇異值分解法和非負最小二乘法。奇異值分解法涉及大量的矩陣運算，移植到井下嵌入式設備中增加硬件邏輯資源，降低算法執行效率，嚴重影響井下快速反演的評價效果。

根據實時性要求并結合原始回波串數據特點，設計一種基于加速投影梯度下降(Accelerated Projected Gradient Descend,APGD)方法的回波T2譜反演算法，并充分利用FPGA硬件資源，采用分布式算法，將APGD算法中的矩陣乘法轉換為查找表和移位累加操作，并在FPGA中對算法進行化簡，減少FPGA硬件電路規模，提高算法執行效率。

1 核磁共振測井反演算法

1.1 反演算法模型

隨鉆核磁測井儀采用CPMG序列測量得到的原始回波串信號是多種橫向弛豫共同作用的結果，信號的波形呈多指數衰減[11]。T2譜信號的離散表達式為

(1)

式中:bi為i時刻的原始回波幅度；n為橫向弛豫分量數目，也表示布點數；xj為對應于T2j分量零時刻信號幅度的貢獻值，也就是反演所要得到的T2譜幅值，此參數必須滿足非負約束條件；ti為i×TE，TE為回波間隔；T2j為第j中弛豫分量的橫向弛豫時間；εi為實際測量產生的隨機噪聲信號；m為原始回波串的回波個數；e-ti/T2j用矩陣表示為

(2)

則式(1)可表示為

bm×1=Am×nxn×1

(3)

式中:A為系數矩陣，與回波設置的參數和布點數有關;b為原始回波串幅度，兩項為已知矩陣，且m>n時此公式為超定方程，這種方程通常是嚴重病態的，其結果非常不穩定[12]。求解x這類不適定問題可以采用非負約束最小二乘(Non-Negative Least Square,NNLS)的思想進行求解[13]。NNLS模型為

(4)

設p=ATb，Q=ATA，則

▽f=Qx-P

(5)

xk=[xk-1-tk(Qxk-1-p)]+

(6)

式中:tk設為1/L，L為Lipschitz常數，L=‖Q‖2。將式(6)展開整理得：

xk+1=[θ1xk+θ2]+

(7)

(8)

式中:k為迭代次數;xk為每次迭代結果;當滿足停止迭代條件時xk為反演結果。

1.2 APGD算法設計

APGD算法結合了Nesterov加速算法和投影梯度下降(Projected Gradient Descend, PGD)法，采用自適應重置參數使算法單調。

FPGA中實現的APGD算法流程如下面的算法1所示，APGD算法共分為兩部分，主體部分執行APGD算法，當殘差值增大時，切換到PGD算法對數據和參數進行更新和重置[14]。

算法 1:APGD for NNLS初始化:x0=y0=0,θ1=In-ATA‖ATA‖2,θ2=ATb‖ATA‖2,k=1,α0∈(0,1)while(不滿足迭代停止條件) do xk=[θ1yk-1+θ2]+(投影梯度步驟) αk=12(α4k-1+4α2k-1-α2k-1),βk=αk-1(1-αk-1)α2k-1+αk yk=xk+βk(xk-xk-1)(外推) if(殘差增加) do xk+1=[θ1xk+θ2]+(投影梯度步驟) yk+1=xk+1(重啟) αk=α0(參數復位) endif k=k+1end

2 FPGA算法實現

引入分布式算法和查找表操作在FPGA中實現APGD算法，并且通過算法的并行處理方式和部分計算步驟的簡化，使得FPGA硬件電路規?？s減，同時算法執行速度加快，實現反演算法在井下的快速處理。

FPGA總體程序設計結構如圖1所示，共分為5個模塊。SPI模塊和Flash模塊為外部通信驅動；APGD模塊負責算法的主體邏輯流程，將矩陣變量組合成地址以備查表，并在算法完成后將反演結果送入主控系統；查找表模塊存儲了預先計算好的矩陣乘法部分積和加速參數結果β，部分積需要后續計算，加速參數結果則送入APGD模塊；矩陣計算模塊將部分積進行移位累加計算得到矩陣乘法結果，并將其返回至APGD模塊參與算法計算。系統各模塊之間由數據線和使能控制信號連接，嚴格保證算法的執行步驟和時序邏輯。

圖1 FPGA程序設計結構圖

2.1 基于分布式算法的累乘加快速計算

APGD算法中存在大量的矩陣與列向量相乘，采用分布式算法完成多組累乘加的快速計算以提高計算效率。分布式算法則通過查找表獲得部分積結果，再對查找表結果進行累加和移位獲得算法結果[15]。一個標準的乘加運算為

(9)

式中:y(n)為累乘加計算結果;Ak為常數;xk(n)為輸入變量。xk(n)采用B位二進制補碼表示為

xB(n),xb(n)∈[0,1]

(10)

式中：xB(n)為x(n)的符號位；xb(n)為x(n)的第b位。將式(10)代入式(9)，可得：

(11)

式(11)將標準的乘加轉換為查找表操作和累加，輸入變量的位作為地址進行查找表映射，在FPGA中將映射值通過移位實現對應的二次冪加權，得到輸出結果。

2.2 多部分表并行計算結構

將完整表分割成多個部分表，部分表并行操作，可以進一步減小查找表規模[16]，令N=L1+L2+…+Lk，則：

(12)

表1 乘加數據舉例

表2 查找表構造規則

圖2為分割查找表的硬件結構框圖，輸入變量x通過拼接構成多組查找表地址，按序移入查找表中完成映射得到部分積，部分積在狀態機的控制下經過累加器和移位寄存器輸出最終結果。

圖2 分割查找表硬件結構圖

在APGD算法中，有3個不同的常數矩陣需要做矩陣乘法，分別為A、AT、θ1，對應列向量累乘加數為N=12、10、10，所以具體構造3組分割查找表為4×23、5×22、5×22。加速參數α和β與迭代次數有關，以迭代次數作為地址構造查找表映射，當參數需要重置時地址重新賦值為零，參數查找表操作省略了開方、除法、方根等運算過程，減少了FPGA邏輯單元的使用，提高了算法處理速度。

在計算殘差的過程中需要開根號計算，本設計FPGA選用ACTEL公司的A3P1000，其不具有對應功能的IP核可以應用。經分析殘差在該算法中僅需要比較大小，不需輸出具體結果，所以在FPGA中用較大空間來存儲開根號前的數據，再進行數據比較。

2.3 仿真測試

采用ModelSim仿真軟件對所設計的APGD算法進行仿真驗證，系統時鐘SYSCLK設置為30 MHz，迭代次數設置為100次，回波數據采用一組隨機數[17]，并與在MATLAB環境下APGD算法的運行結果作對比說明。

圖3為FPGA下APGD算法的仿真結果，在迭代第60次時取出殘差值c_a[0～11]和反演結果值x[0～9]。殘差值基本保持不變，選取第60次的結果作為仿真測試的結果。

圖4為MATLAB和FPGA的算法殘差歸一化收斂曲線對比，從圖4中看出兩種實現方式對應的殘差衰減基本相同，以驗證APGD算法的收斂性和FPGA實現此算法的可行性。

圖4 PC端與FPGA殘差迭代

3 實驗結果與分析

3.1 實驗環境與實驗裝置

正在研制的隨鉆核磁共振測井儀器處于實驗室樣機驗證階段，為驗證回波反演算法的性能，在儀器的回波采集和數據處理電路中運行所設計的回波快速反演算法?；夭ù杉蛿祿幚黼娐啡鐖D5所示，電路采用DSP+FPGA的結構作為電路的核心，FPGA選用ProASIC3系列的A3P1000，實現回波串采集和基于APGD方法的回波反演計算功能。DSP采用TMS320F2812，通過RS485接口將數據發送至主控電路。

圖5 回波串采集和數據處理電路功能框圖

測試中配置一定濃度的CuSO4溶液加入刻度箱中，使用刻度箱環境模擬真實地層特性，設置100%孔隙度的地層信息作為原始模擬數據，用于算法結果的對比分析。

3.2 實驗操作與過程

實驗中采集電路通過探頭向外界發射CPMG(Garr-Purcell-Meiboom-Gill)脈沖序列，設置回波間隔TE=0.6 ms。電路接收到原始回波串后，設置橫向弛豫時間在區間0.5～5000 ms上均勻布點，布點數n=10。已知以上參數，可根據式(2)計算出Am×n。對原始回波串信號進行預處理校正和信號抽樣，抽樣個數m=12，回波幅度為bm×1。一組典型的抽樣后的回波衰減曲線如圖6所示。

圖6 原始回波串信號

3.3 實驗數據分析

實驗條件下刻度箱設為100%孔隙度的地層特性，地層孔隙度為T2譜幅值累加值。PC端與FPGA反演結果對比如圖7所示，可以看出所得反演結果T2譜呈單峰性，MATLAB和FPGA反演結果的孔隙度分別為97.47%和96.81%。將預設的孔隙度數據Xmod與反演結果Xinv代入式(13)中，計算相對誤差RX。

圖7 PC端與FPGA反演結果對比

(13)

在MATLAB和FPGA中反演相對誤差分別為2.53%和3.19%，可見使用實際回波串數據驗證該算法其結果精度仍然較高，可以滿足井下快速反演的需求。

3.4 方法設計評估

表3為A3P1000中的片上IP資源使用情況和所占總資源的比率。片上共有24567觸發器資源和144 KB的RAM，通過分布式算法簡化分割查找表，僅消耗15.97%的RAM資源。

表3 IP資源和內存消耗

分別在PC端、DSP、FPGA上測試此算法執行時間，對算法耗時進行對比。PC端算法實驗在Win10系統下進行配有2.40 GHz的Intel Core i5-10200H處理器。DSP和FPGA均在30 MHz頻率系統時鐘下工作。分別對APGD算法的單次執行不更新參數、單次執行更新參數和完整迭代流程進行計時，算法耗時對比如表4所示，從耗時來看DSP>PC端>FPGA。

表4 FPGA和DSP算法耗時比較 單位：μs

通過以上實驗算法效率和結果精度對比可以看出，FPGA作為反演算法在井下快速實現的嵌入式設備具有十分重要的意義和研究價值，對比其他現有工作，在FPGA上實現APGD算法在處理效率和性能上均具有較大優勢。

4 結束語

針對隨鉆核磁共振測井儀器的井下快速反演問題，提出了一種基于FPGA的APGD算法實現方法。該方法應用了分布式算法，將矩陣乘法中復雜的累乘加運算轉化為查找表操作，同時在FPGA中簡化原算法，不僅提高了算法的執行效率，而且節約了硬件電路規模。通過自主研發的隨鉆核磁共振測井儀實驗室刻度實驗結果表明，該方法相對誤差較低，算法執行效率較高，實時性和準確率滿足隨鉆核磁測井儀器快速反演的需求，為核磁共振井下反演提供一種可行的快速實現方式。