李沐堯，朱萬(wàn)里，滕 朔，程晶晶*，葸春平

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

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

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

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

1 核磁共振測(cè)井反演算法

1.1 反演算法模型

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

(1)

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

(2)

則式(1)可表示為

bm×1=Am×nxn×1

(3)

式中:A為系數(shù)矩陣，與回波設(shè)置的參數(shù)和布點(diǎn)數(shù)有關(guān);b為原始回波串幅度，兩項(xiàng)為已知矩陣，且m>n時(shí)此公式為超定方程，這種方程通常是嚴(yán)重病態(tài)的，其結(jié)果非常不穩(wěn)定[12]。求解x這類不適定問(wèn)題可以采用非負(fù)約束最小二乘(Non-Negative Least Square,NNLS)的思想進(jìn)行求解[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為L(zhǎng)ipschitz常數(shù)，L=‖Q‖2。將式(6)展開(kāi)整理得：

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

(7)

(8)

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

1.2 APGD算法設(shè)計(jì)

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

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

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

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

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

2.1 基于分布式算法的累乘加快速計(jì)算

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

(9)

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

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

(10)

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

(11)

式(11)將標(biāo)準(zhǔn)的乘加轉(zhuǎn)換為查找表操作和累加，輸入變量的位作為地址進(jìn)行查找表映射，在FPGA中將映射值通過(guò)移位實(shí)現(xiàn)對(duì)應(yīng)的二次冪加權(quán)，得到輸出結(jié)果。

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

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

(12)

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

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

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

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

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

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

2.3 仿真測(cè)試

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

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

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

圖4 PC端與FPGA殘差迭代

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)環(huán)境與實(shí)驗(yàn)裝置

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

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

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

3.2 實(shí)驗(yàn)操作與過(guò)程

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

圖6 原始回波串信號(hào)

3.3 實(shí)驗(yàn)數(shù)據(jù)分析

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

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

(13)

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

3.4 方法設(shè)計(jì)評(píng)估

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

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

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

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

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

4 結(jié)束語(yǔ)

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