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

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

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

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

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

1 核磁共振測井反演算法

1.1 反演算法模型

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

(1)

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

(2)

則式(1)可表示為

bm×1=Am×nxn×1

(3)

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

(4)

設(shè)p=ATb，Q=ATA，則

▽f=Qx-P

(5)

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

(6)

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

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

(7)

(8)

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

1.2 APGD算法設(shè)計

APGD算法結(jié)合了Nesterov加速算法和投影梯度下降(Projected Gradient Descend, PGD)法，采用自適應(yīng)重置參數(shù)使算法單調(diào)。

FPGA中實現(xiàn)的APGD算法流程如下面的算法1所示，APGD算法共分為兩部分，主體部分執(zhí)行APGD算法，當殘差值增大時，切換到PGD算法對數(shù)據(jù)和參數(shù)進行更新和重置[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(參數(shù)復(fù)位) endif k=k+1end

2 FPGA算法實現(xiàn)

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

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

圖1 FPGA程序設(shè)計結(jié)構(gòu)圖

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

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

(9)

式中:y(n)為累乘加計算結(jié)果;Ak為常數(shù);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)將標準的乘加轉(zhuǎn)換為查找表操作和累加，輸入變量的位作為地址進行查找表映射，在FPGA中將映射值通過移位實現(xiàn)對應(yīng)的二次冪加權(quán)，得到輸出結(jié)果。

2.2 多部分表并行計算結(jié)構(gòu)

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

(12)

表1 乘加數(shù)據(jù)舉例

表2 查找表構(gòu)造規(guī)則

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

圖2 分割查找表硬件結(jié)構(gòu)圖

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

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

2.3 仿真測試

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

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

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

圖4 PC端與FPGA殘差迭代

3 實驗結(jié)果與分析

3.1 實驗環(huán)境與實驗裝置

正在研制的隨鉆核磁共振測井儀器處于實驗室樣機驗證階段，為驗證回波反演算法的性能，在儀器的回波采集和數(shù)據(jù)處理電路中運行所設(shè)計的回波快速反演算法?；夭ù杉蛿?shù)據(jù)處理電路如圖5所示，電路采用DSP+FPGA的結(jié)構(gòu)作為電路的核心，F(xiàn)PGA選用ProASIC3系列的A3P1000，實現(xiàn)回波串采集和基于APGD方法的回波反演計算功能。DSP采用TMS320F2812，通過RS485接口將數(shù)據(jù)發(fā)送至主控電路。

圖5 回波串采集和數(shù)據(jù)處理電路功能框圖

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

3.2 實驗操作與過程

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

圖6 原始回波串信號

3.3 實驗數(shù)據(jù)分析

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

圖7 PC端與FPGA反演結(jié)果對比

(13)

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

3.4 方法設(shè)計評估

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

表3 IP資源和內(nèi)存消耗

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

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

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

4 結(jié)束語

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