基于數(shù)字體圖像相關(guān)法的TATB 基PBX 材料內(nèi)部變形測(cè)量

王延珺，鄒 翔，潘 兵，郭昱辰，劉 晨，韓 超，呂珂臻

（1. 中國(guó)工程物理研究院化工材料研究所，四川 綿陽(yáng) 621999；2. 北京航空航天大學(xué)固體力學(xué)所，北京 100191）

0 引言

高聚物黏結(jié)炸藥（Polymer Bonded Explosives，PBX）作為一種高能炸藥晶體和高聚物黏結(jié)劑組成的復(fù)合材料，因其具有能量高、鈍感、耐熱的優(yōu)點(diǎn)，在國(guó)防工業(yè)中得到了廣泛應(yīng)用。PBX 材料中的炸藥晶體和黏結(jié)劑在含量、熱物理性能方面均存在較大差異，因此組分之間存在性質(zhì)各異的界面。同時(shí)，PBX 材料是一種典型的脆性材料，其抗拉強(qiáng)度較低，力學(xué)性能極易受振動(dòng)、沖擊、熱加載等環(huán)境因素影響［1-3］。將PBX 材料壓制成型為部件時(shí)，內(nèi)應(yīng)力釋放機(jī)制會(huì)使成型部件發(fā)生彈性回彈，影響其穩(wěn)定性。此外，PBX 材料在熱處理、機(jī)械加工、運(yùn)輸及存儲(chǔ)過(guò)程中產(chǎn)生的內(nèi)應(yīng)力會(huì)使晶體斷裂、界面脫粘、黏結(jié)劑基體開(kāi)裂，造成PBX 材料力學(xué)性能的下降［4-5］，因此極需了解PBX 內(nèi)應(yīng)力演化規(guī)律及微細(xì)觀結(jié)構(gòu)演化機(jī)理的識(shí)別方法。

現(xiàn)有PBX 的原位力學(xué)性能測(cè)試方法主要有接觸式與非接觸式，然而都存在一定局限性。如采用表面貼片應(yīng)變計(jì)測(cè)量炸藥部件的力學(xué)響應(yīng)時(shí)，由于PBX 材料屬于非均勻材料，所獲得的測(cè)量值僅能反映小范圍被測(cè)表面的平均效應(yīng)；若采用預(yù)埋傳感器測(cè)量，介入式傳感器會(huì)改變材料內(nèi)部結(jié)構(gòu)，引入較大的測(cè)量誤差。非接觸方法相比于以上接觸式方法，能夠極大地減小測(cè)量誤差，主要有射線(xiàn)法、光測(cè)力學(xué)法。基于布拉格晶格衍射原理的射線(xiàn)法（X/中子衍射）僅能獲得淺表層小體積內(nèi)的晶格畸變，從而可以反映出束斑區(qū)域的晶體變形信息［6］。然而，對(duì)PBX 材料的裂紋檢測(cè)多見(jiàn)于晶體與高分子黏結(jié)劑的界面失效［2］。因此，射線(xiàn)法無(wú)法反映高分子粘結(jié)劑對(duì)PBX 復(fù)合材料的三維應(yīng)變應(yīng)力響應(yīng)的影響。近10 余年來(lái)，利用光測(cè)力學(xué)手段，如數(shù)字圖像相關(guān)方法（Digital Image Correlation，DIC）測(cè)量PBX 材料在熱力作用下的表面位移場(chǎng)及應(yīng)變場(chǎng)開(kāi)始興起［7］。然而2D-DIC，3D-DIC 僅可用于表面變形測(cè)量［8］，獲取的是材料表面位移的2 個(gè)或3 個(gè)自由度分量。非均勻材料的表面變形場(chǎng)分布與內(nèi)部變形場(chǎng)分布往往有顯著的差異，僅依賴(lài)表面的位移分析并不能得到非均勻材料體系內(nèi)部變形的準(zhǔn)確信息。故而，常規(guī)測(cè)試手段無(wú)法實(shí)現(xiàn)PBX 這種多晶、多界面、非均質(zhì)的復(fù)雜材料體系在外載荷作用下的的三維力學(xué)響應(yīng)表征。

基于X 射線(xiàn)微焦點(diǎn)計(jì)算機(jī)斷層掃描（μ-CT）技術(shù)和數(shù)字體圖像相關(guān)（Digital Volume Correlation，DVC）算法的物體內(nèi)部變形測(cè)量技術(shù)是實(shí)驗(yàn)固體力學(xué)的前沿研究方法。作為可用于不透明固體材料內(nèi)部變形場(chǎng)測(cè)量的實(shí)驗(yàn)力學(xué)技術(shù)，DVC 方法在實(shí)驗(yàn)力學(xué)、巖土力學(xué)、材料科學(xué)、生物醫(yī)學(xué)等領(lǐng)域受到關(guān)注［9］，在含能材料領(lǐng)域內(nèi)研究剛剛起步。美國(guó)Los Alamos 國(guó)家實(shí)驗(yàn)室的研究團(tuán)隊(duì)利用IMPULSE 同步輻射先進(jìn)光源APS，對(duì)奧克托今（HMX）基PBX 炸藥在動(dòng)態(tài)壓縮實(shí)驗(yàn)進(jìn)行了表征和追蹤，測(cè)量了沖擊過(guò)程的晶體破碎形貌，但未能對(duì)應(yīng)力應(yīng)變演化進(jìn)行深入分析［10］。法國(guó)巴黎薩克雷大學(xué)的Hild［11］與法國(guó)軍用核能中心CEA-Gramat 研究團(tuán)隊(duì)將DVC 方法應(yīng)用于ESRF 加速器CT 原位觀測(cè)下的HTPB 推進(jìn)劑材料單軸拉伸實(shí)驗(yàn)，通過(guò)計(jì)算全局體積應(yīng)變獲得材料在不同載荷狀態(tài)下的泊松比，通過(guò)局部分析各計(jì)算單元內(nèi)的平均應(yīng)變獲得模擬材料晶體與黏結(jié)劑2 種主要組分的非均勻變形狀態(tài)，并且DVC 法的灰度殘余分布圖能夠精確地表明材料的損傷程度與開(kāi)裂機(jī)理，該方法對(duì)于揭示含能材料內(nèi)部變形機(jī)制具有重要意義［11］。英國(guó)牛津大學(xué)的Siviour 等［12］借助原位CT 實(shí)驗(yàn)，利用DVC 法獲得了含能材料造型粉顆粒在壓制過(guò)程中的三維位移場(chǎng)、應(yīng)變場(chǎng)等高精度信息，DVC 方法獲得的材料單軸壓縮力學(xué)響應(yīng)與有限元仿真模型數(shù)值預(yù)測(cè)結(jié)果符合度良好，為進(jìn)一步推導(dǎo)應(yīng)力場(chǎng)分布提供了必要條件。該團(tuán)隊(duì)與Los Alamos 國(guó)家驗(yàn)室將DVC 方法應(yīng)用于PBX 試樣在單軸壓縮過(guò)程中的變形測(cè)量，揭示了內(nèi)部不均勻應(yīng)變場(chǎng)分布與細(xì)觀結(jié)構(gòu)演變的關(guān)聯(lián)關(guān)系，挖掘出豐富的力學(xué)信息：內(nèi)部單個(gè)造型粉顆粒的運(yùn)動(dòng)，內(nèi)部裂紋的形成與拓展，黏結(jié)劑與晶體的脫粘以及試樣的斷裂機(jī)制。美國(guó)德州大學(xué)達(dá)拉斯分校的Hu 等［13］采用DVC 方法測(cè)量了曼森砂在圍壓實(shí)驗(yàn)中的內(nèi)部應(yīng)變場(chǎng)，并通過(guò)分析各個(gè)顆粒的最小本征應(yīng)變值評(píng)估了壓縮過(guò)程中的力鏈形成與演化過(guò)程［14］，對(duì)于揭示顆粒材料體系壓制成型的機(jī)理機(jī)制具有重要意義。北京航空航天大學(xué)的潘兵等［15］利用自己開(kāi)發(fā)的增量DVC 算法系統(tǒng)研究了推進(jìn)劑代用材料在單軸壓縮實(shí)驗(yàn)過(guò)程中壓縮應(yīng)變從10%增加到70%過(guò)程中其內(nèi)部應(yīng)變場(chǎng)變化情況。由此可見(jiàn)，該方法能夠定性定量地分析含能材料在原位力學(xué)實(shí)驗(yàn)中的變形狀態(tài)，對(duì)于開(kāi)展相關(guān)研究具有一定的參考意義。當(dāng)晶體顆粒填充比例較低時(shí)，所獲得的數(shù)字圖像灰度對(duì)比度大，對(duì)變形狀態(tài)的測(cè)量相對(duì)容易實(shí)現(xiàn)。

為此，本研究將基于μ-CT 和DVC 結(jié)合的方法應(yīng)用于TATB 基PBX 材料內(nèi)部變形測(cè)量的探索。由于PBX 內(nèi)部組分密度接近，數(shù)字圖像灰度值對(duì)比度相對(duì)較低，在小變形狀態(tài)下即會(huì)發(fā)生斷裂，因而研究在采用μ-CT 對(duì)單軸壓縮下的PBX 圓柱試件進(jìn)行了體成像掃描后，為提高內(nèi)部組分紋理對(duì)比度，將圖像進(jìn)行預(yù)處理，采用自編程的基于C++語(yǔ)言的局部DVC 算法，對(duì)所得PBX 的CT 成像的灰度圖原始數(shù)據(jù)進(jìn)行計(jì)算，揭示了TATB 基PBX 材料在單軸壓縮實(shí)驗(yàn)中的三維內(nèi)部變形狀態(tài)，并提出了通過(guò)DVC 計(jì)算獲得的零均值歸一化互相關(guān)（zero-mean normalized cross-correlatio，ZNCC）系數(shù)分布圖來(lái)判定亞體素裂紋的空間位置的方法。

1 原位實(shí)驗(yàn)與圖像分析方法

1.1 TATB 基PBX 材料單軸壓縮實(shí)驗(yàn)

樣品試件：TATB 基PBX 試件為120 ℃、150 MPa壓制的直徑為10 mm、高度為10.2 mm 藥柱（圖1），藥柱密度為1.876 g·cm-3。

圖1 TATB 基PBX 材料試件在單軸壓縮下的力-位移曲線(xiàn)Fig.1 Loading curve of TATB-based PBX sample under uniaxial compression test

儀器：?jiǎn)屋S壓縮實(shí)驗(yàn)采用配置Deben Microtest CT5000-TEC 原位加載裝置的原位μ-CT 進(jìn)行。

實(shí)驗(yàn)方法：室溫下對(duì)TATB 基PBX 試件進(jìn)行同步的單軸壓縮實(shí)驗(yàn)和CT 觀測(cè)。為采集壓縮過(guò)程中的變化圖像，實(shí)驗(yàn)對(duì)樣件加載后，再進(jìn)行CT 掃描。由于應(yīng)力松弛，載荷均會(huì)出現(xiàn)下降，因此不同的時(shí)刻進(jìn)行CT 掃描后繼續(xù)進(jìn)行加載，如此循環(huán)加載掃描步驟，直至試件發(fā)生破壞后停止實(shí)驗(yàn)，本研究為此連續(xù)進(jìn)行5 次加載。加載模式為位移加載模式，加載速率為0.1 mm·min-1，實(shí)驗(yàn)力-位移載荷如圖1 所示。CT 掃描電壓為100 kV，電流為100 μA，曝光時(shí)間為0.6 ms，試件旋轉(zhuǎn)過(guò)程中每隔0.33°采集1 張圖像，共采集1080 張圖像，獲取的體 圖 像 尺 寸 為550 voxel×550 voxel×500 voxel。實(shí)驗(yàn)過(guò)程中保持CT 掃描參數(shù)相同，以減小測(cè)量參數(shù)對(duì)結(jié)果的影響，保證在相同條件下對(duì)試件內(nèi)部形變帶來(lái)的灰度變化進(jìn)行追蹤，進(jìn)而實(shí)現(xiàn)對(duì)內(nèi)部變形及損傷演化的定量分析。實(shí)驗(yàn)中調(diào)節(jié)樣品與探測(cè)器之間的距離，使TATB 基PBX 材料的CT 體圖像分辨率為19.8 μm。

1.2 數(shù)字體圖像相關(guān)法(DVC)

數(shù)字體圖像相關(guān)法（DVC）可通過(guò)對(duì)比分析數(shù)字體成像設(shè)備得到的試件變形前后的2 組數(shù)字體圖像，經(jīng)由計(jì)算感興趣區(qū)域內(nèi)的灰度變化獲得物體內(nèi)部亞體素精度的三維位移場(chǎng)和全場(chǎng)應(yīng)變［9］。DVC 法分析主要包括以下3 個(gè)步驟：

（1）數(shù)字體圖像的采集：采用體成像設(shè)備對(duì)不同加載狀態(tài)下的試件進(jìn)行掃描，得到試件的數(shù)字體圖像；（2）三維位移場(chǎng)的測(cè)量：通過(guò)體圖像配準(zhǔn)算法［9］確定參考體圖像中各離散計(jì)算點(diǎn)在目標(biāo)體圖像中的位置，獲得全場(chǎng)三維位移；（3）三維應(yīng)變場(chǎng)計(jì)算：利用數(shù)值微分方法處理得到三維位移場(chǎng)，獲得三維應(yīng)變場(chǎng)。

為此，本研究通過(guò)單軸壓縮實(shí)驗(yàn)獲取了4 種載荷狀態(tài)下的試件體圖像。為準(zhǔn)確研究試件變形引起的部分計(jì)算點(diǎn)灰度變化，同時(shí)對(duì)圖像進(jìn)行了預(yù)處理，以提升體圖像紋理對(duì)比度，進(jìn)而保證DVC 計(jì)算精度。用DVC 對(duì)原始圖片直接進(jìn)行處理，在4 種載荷下，基于C++編程的DVC 算法對(duì)原始圖片的計(jì)算速度分別為8.26 點(diǎn)/秒，7.53 點(diǎn)/秒，6.22 點(diǎn)/秒 和5.11 點(diǎn)/秒，從 計(jì)算時(shí)間看，隨著載荷的加大，試樣變形越大，故而所需計(jì)算時(shí)間越長(zhǎng)。

三維應(yīng)變場(chǎng)獲取的關(guān)鍵是三維位移場(chǎng)的測(cè)量精準(zhǔn)度。為確定某計(jì)算點(diǎn)的三維位移矢量，需要圍繞該點(diǎn)選擇包含足夠灰度對(duì)比度的立方體參考圖像子體塊，并建立一個(gè)可定量評(píng)價(jià)參考圖像子體塊和待尋找的目標(biāo)圖像子體塊相似程度的目標(biāo)函數(shù)。本研究采用零均值歸一化互相關(guān)（zero-mean normalized sum-of-square difference，ZNSSD）函數(shù)來(lái)評(píng)價(jià)2 子體塊的相似程度［18］。由于目標(biāo)圖像子體塊可能出現(xiàn)包括平移、轉(zhuǎn)動(dòng)、均勻應(yīng)變甚至更為復(fù)雜的非均勻變形，因而需要采用包括12 個(gè)變量的一階形函數(shù)［16］來(lái)近似其內(nèi)部變形，再通過(guò)數(shù)值優(yōu)化算法和灰度插值算法相結(jié)合來(lái)優(yōu)化目標(biāo)函數(shù)，以在變形后的體圖像中尋找與圖像子體塊最為相似的目標(biāo)圖像子體塊，從而確定其當(dāng)前體素點(diǎn)（參考圖像子體塊中心點(diǎn)）的位移矢量。通過(guò)簡(jiǎn)單的整體素搜索，DVC 能快速確定各計(jì)算點(diǎn)的整體素位移。隨后，為了提高測(cè)量精度，需要采用亞體素配準(zhǔn)算法對(duì)整像素結(jié)果進(jìn)行優(yōu)化以獲得具有亞體素精度的位移結(jié)果，這一過(guò)程通常被認(rèn)為是提高DVC 測(cè)量精準(zhǔn)度的關(guān)鍵［9］。在眾多亞體素位移算法中，反向組合高斯-牛頓算法（Inverse Compositional Gauss-Newton，IC-GN）因?yàn)榫哂杏?jì)算精度和計(jì)算效率高等優(yōu)點(diǎn)被廣泛應(yīng)用。本研究采用最新的三維IC-GN 配準(zhǔn)算法，插值方法采用3 次B 樣條插值方法，形函數(shù)采用一階形函數(shù)，應(yīng)變計(jì)算采用逐點(diǎn)最小二乘（pointwise least-square，PLS）算法進(jìn)行運(yùn)算。

為了保證計(jì)算效率和精度，使計(jì)算能夠在普通計(jì)算機(jī)中使用DVC 方法對(duì)高分辨率體圖像進(jìn)行分析，本次分析采用逐層可靠性導(dǎo)向的初值估計(jì)策略為計(jì)算點(diǎn)提供位移初值［17］，由于在加載過(guò)程中，試樣沒(méi)有發(fā)生較大的旋轉(zhuǎn)運(yùn)動(dòng)，因此可靠性導(dǎo)向種子點(diǎn)的位移初值采用簡(jiǎn)單的整體素搜索方法。與現(xiàn)有文獻(xiàn)中提出的常規(guī)三維可靠性導(dǎo)向位移跟蹤方法不同，逐層可靠性導(dǎo)向方法僅僅引導(dǎo)每個(gè)計(jì)算層內(nèi)的計(jì)算路徑并傳遞相鄰計(jì)算層之間的變形矢量以保證初值精度。

2 結(jié)果與討論

2.1 體圖像預(yù)處理對(duì)DVC 結(jié)果的影響

在DVC 計(jì)算中，ZNCC 系數(shù)（變化范圍為［-1，1］）被廣泛應(yīng)用于判斷DVC 計(jì)算點(diǎn)的匹配程度，ZNCC 系數(shù)越高，說(shuō)明計(jì)算點(diǎn)的參考子體塊和變形子體塊匹配越好，此時(shí)獲得的位移場(chǎng)精度越高，相反地，當(dāng)一計(jì)算點(diǎn)的ZNCC 系數(shù)低于某一設(shè)定值時(shí)就認(rèn)為該點(diǎn)匹配失敗，位移場(chǎng)誤差較大［9］。為了獲取較高的ZNCC 系數(shù)，研究采用降噪方法對(duì)試件的體圖像圖像進(jìn)行了預(yù)處理，并對(duì)比分析了預(yù)處理前后的DVC 計(jì)算結(jié)果。

圖2 為原始CT 體圖像預(yù)處理前，DVC 計(jì)算得到的各計(jì)算點(diǎn)處的ZNCC 系數(shù)分布圖。由圖2 可以看出，在4 種載荷下，各計(jì)算點(diǎn)ZNCC 系數(shù)均較低，均值僅分別為0.696、0.672、0.670 和0.668，這主要是由于試件體圖像內(nèi)部體圖像特征不明顯，抗噪聲能力較弱；此外，隨著載荷的增加，計(jì)算點(diǎn)ZNCC 系數(shù)均值呈現(xiàn)小幅度下降趨勢(shì)，分析認(rèn)為這主要是因?yàn)樵嚰冃卧斐闪瞬糠钟?jì)算點(diǎn)的灰度變化?？梢?jiàn)，直接對(duì)試件CT 的原始圖片進(jìn)行DVC 計(jì)算，所得的結(jié)果較差。因此接下來(lái)將會(huì)對(duì)原始圖像進(jìn)行預(yù)處理以提高DVC 結(jié)果的準(zhǔn)確性。

圖2 采用原始CT 體圖像不同載荷下各計(jì)算點(diǎn)ZNCC 系數(shù)分布情況Fig.2 Distribution of ZNCC coefficient using original CT volume images under different loading states

在DIC 測(cè)量中，操作者可在被測(cè)試樣表面制作人工散斑圖案以提供高對(duì)比度的變形信息載體。在絕大多數(shù)情況下，DVC 方法只能依賴(lài)被測(cè)物體內(nèi)部的復(fù)雜成分或微結(jié)構(gòu)體系在體圖像中形成的灰度差異作為變形信息載體，因此體圖像的質(zhì)量對(duì)DVC 的計(jì)算精度影響較大。相比于具有良好對(duì)比度的材料（如巖石、泡沫等）以及具有大量?jī)?nèi)含物的材料（如球墨鑄鐵、鋁銅合金等），PBX 材料組分密度接近，內(nèi)部紋理天然較弱。為了提高DVC 的測(cè)量精準(zhǔn)度對(duì)圖像進(jìn)行處理，圖3 為預(yù)處理前后灰度分布直方圖，圖4 為預(yù)處理前后試件體圖像和切片圖。由圖3a 可以看出，原始圖像中大量體素點(diǎn)灰度值集中在0～150，導(dǎo)致圖像整體偏暗，難以辨別（圖4a）。為了解決這個(gè)問(wèn)題，將灰度值大于130的體素點(diǎn)直接設(shè)為0（這些點(diǎn)屬于噪聲點(diǎn)），然后將剩余體素點(diǎn)灰度值拉伸到0～255，拉伸公式見(jiàn)式（1）。由圖4b 可知，實(shí)驗(yàn)試件內(nèi)的灰度進(jìn)行優(yōu)化處理后，對(duì)比度有所提升，可以辨別試件內(nèi)的顆粒邊界。剩下4組變形體圖像均進(jìn)行類(lèi)似操作。

圖3 實(shí)驗(yàn)試件體圖像優(yōu)化前后灰度分布直方圖Fig.3 Histograms of original CT image and optimized CT image

圖4 優(yōu)化前后的試件體圖像和切片圖Fig.4 Comparison between original CT image and optimized CT image

圖5 為實(shí)驗(yàn)試件內(nèi)的灰度進(jìn)行優(yōu)化處理后，4 種不同載荷下各計(jì)算點(diǎn)ZNCC 系數(shù)分布情況。由于PBX 材料是脆性材料，在外載荷作用下，發(fā)生較小的變形即會(huì)斷裂，故而圖5 中ZNCC 系數(shù)分布大體一致。由圖5可以看出各計(jì)算點(diǎn)ZNCC 系數(shù)均得到提升，均值高于0.75。這是因?yàn)樽鳛樽冃涡畔⑤d體的內(nèi)部紋理的對(duì)比度提高，體圖像的噪音大大降低，可見(jiàn)，預(yù)處理對(duì)于提升DVC 計(jì)算精度具有較大的影響。

圖5 不同載荷下各計(jì)算點(diǎn)ZNCC 系數(shù)分布情況Fig.5 Distribution of ZNCC coefficient of optimized CT images under different loading states

2.2 單軸壓縮實(shí)驗(yàn)過(guò)程中的位移場(chǎng)

在進(jìn)行DVC 分析前，為了降低圖像噪聲對(duì)測(cè)量精度的影響，利用高斯濾波器（尺寸為5×5×5，標(biāo)準(zhǔn)差為1）對(duì)4 組體圖像進(jìn)行高斯濾波［18］。然后，在濾波后的參考體圖像中心選取尺寸為300 voxel×300 voxel×300 voxel 的感興趣區(qū)域（Volume of Interest，VOI）。在DVC 計(jì)算時(shí)，將每次加載獲得的體圖像（變形體圖像）分別與加載前采集的體圖像（參考體圖像）進(jìn)行匹配，采用3D IC-GN 亞體素配準(zhǔn)算法，最優(yōu)收斂條件為：位移矢量增量的模小于0.001 voxel 或迭代次數(shù)不超過(guò)20 次，子體塊尺寸為41 voxel×41 voxel×41 voxel，計(jì)算步長(zhǎng)為10 voxel，共有30×30×30=27000 個(gè)均勻分布的計(jì)算點(diǎn)。在4 種載荷下，DVC 計(jì)算速度分別為8.31，8.02，7.31 點(diǎn)/秒和6.45 點(diǎn)/秒，可以看出試樣變形越大，所需計(jì)算時(shí)間越長(zhǎng)。

在4 種載荷下，大部分計(jì)算點(diǎn)ZNCC 系數(shù)較高，4種載荷下ZNCC 系數(shù)均值分別為0.798，0.772，0.770和0.768，均在0.75 以上，因此結(jié)果具有較高的可靠性；與此同時(shí)，隨著載荷的增加，計(jì)算點(diǎn)ZNCC 系數(shù)均值呈小幅度下降趨勢(shì)，這主要是由于試件變形引起部分計(jì)算點(diǎn)灰度變化而引起的。為了保證計(jì)算結(jié)果的可靠性，需要將ZNCC 系數(shù)較低的計(jì)算點(diǎn)剔除。根據(jù)計(jì)算 結(jié) 果，將0.75 設(shè) 為ZNCC 閾 值，ZNCC 低 于 該 閾 值的計(jì)算點(diǎn)將被剔除，這些點(diǎn)的位移結(jié)果（包括位移u，v，w）將通過(guò)周?chē)鶽NCC 系數(shù)高于0.75 的計(jì)算點(diǎn)插值得到。在4 種載荷下，ZNCC 系數(shù)高于所設(shè)閾值的點(diǎn)分別為23325、19356、19356 個(gè)和18497 個(gè)，占總共27000 個(gè)點(diǎn)的比例分別為86.3889%、71.6889%、70.0185 %和68.5074 %，均大于65%，因此計(jì)算結(jié)果具有較高的可靠性。

圖6 給出了在不同載荷下各計(jì)算點(diǎn)的u、v、w位移場(chǎng)。由圖6 可知，獲得的位移場(chǎng)表現(xiàn)出極大的非均勻性，主要原因?yàn)樵嚰?nèi)部結(jié)構(gòu)各向異性；造型粉顆粒的尺寸及大小呈現(xiàn)非均勻分布，因此在加載過(guò)程中各點(diǎn)位移差別較大。與此同時(shí)，由圖6 可知，在載荷1 下位移u的最大值和最小值分別為10.9906 voxel 和9.0158 voxel，在載荷2 下位移u的最大值和最小值為11.2879 voxel 和7.8982 voxel，在 載 荷3 下 位 移u的最大值和最小值為12.4122 voxel 和7.8921 voxel，在載荷4 下位移u的最大值和最小值為11.9612 voxel和8.0347 voxel。在4 種載荷下，相對(duì)位移u而言，位移v值較小，其變化范圍依次為［-2.6380 voxel，0.1821 voxel］、［3.0077 voxel，9.4698 voxel］ 、［3.2914 voxel，9.3240 voxel］和［2.4368 voxel，8.9796 voxel］，這說(shuō)明在加載過(guò)程中試件在x方向還存在滑移。相對(duì)于位移u、v，位移w的變化更加均勻，這主要是由于z方向是加載方向，各點(diǎn)變形比較大，結(jié)果信噪比較高。

2.3 應(yīng)變應(yīng)力結(jié)果的分析

研究通過(guò)對(duì)每次加載獲得的變形體圖像與參考體圖像的灰度進(jìn)行亞體素配準(zhǔn)計(jì)算分析，獲得了不同載荷下各計(jì)算點(diǎn)的u、v、w位移場(chǎng)，結(jié)果如圖6 所示。由圖6 可以看出隨著單軸壓縮載荷的增加，感興趣區(qū)域內(nèi)的位移逐步增大，且由于非均勻材料體系特征，位移分布也體現(xiàn)出不均勻性。然而材料的力學(xué)性能由其受外載荷下的應(yīng)力應(yīng)變響應(yīng)所反映，基于獲得的三維位移場(chǎng)采用合適的差分算法即可獲得應(yīng)變場(chǎng)。

圖6 不同載荷下位移場(chǎng)Fig. 6 Displacement fields under different loading states

由于各種因素的影響，該三維位移場(chǎng)不可避免地包含噪聲，因此本研究采用逐點(diǎn)局部最小二乘擬合法計(jì)算各計(jì)算點(diǎn)的位移場(chǎng)，應(yīng)變窗口大小為11×11×11。由于應(yīng)變計(jì)算時(shí)間較短，約為5 s，因此不再給出具體的應(yīng)變計(jì)算時(shí)間。為了更加直觀地展示應(yīng)變結(jié)果，結(jié)果如表1 所示，試件內(nèi)部6 個(gè)應(yīng)變分量（εxx，εyy，εzz，γxy，γxz，γyz）在不同載荷下的平均值。

由表1 可知，最大應(yīng)變?yōu)閱屋S載荷方向εzz，當(dāng)載荷穩(wěn)定以后，從載荷2 開(kāi)始，平面內(nèi)應(yīng)變?chǔ)舩x?εyy，3 個(gè)切應(yīng)變分量的變化較小，且均維持在較低水平。試件沿x、y方向的主應(yīng)變相差不大，這說(shuō)明試件沿壓制面內(nèi)的2個(gè)方向的力學(xué)性能差別較小，應(yīng)變測(cè)量結(jié)果符合模壓?jiǎn)蜗驂嚎s成型的PBX 材料的力學(xué)特征。

表1 不同載荷下試件內(nèi)部各應(yīng)變分量平均值Table 1 Mean value of strain component under different loading states

根據(jù)加載過(guò)程中的載荷力-位移關(guān)系，表2 給出了每個(gè)加載步中應(yīng)力-應(yīng)變曲線(xiàn)的擬合直線(xiàn)，由表2 可知，由載荷0 狀態(tài)到載荷1 過(guò)程中，擬合直線(xiàn)的斜率為532.25 MPa；從載荷1 到載荷2，擬合直線(xiàn)的斜率為649.13 MPa；從載荷2 到載荷3，擬合直線(xiàn)的斜率為1317.00 MPa；從載荷3 到載荷4，擬合直線(xiàn)的斜率為1372.30 MPa。

表2 不同載荷階段應(yīng)力應(yīng)變線(xiàn)性擬合曲線(xiàn)Table 2 Linear fitting curve of stress-strain relationship under different loading states

通過(guò)前期材料力學(xué)實(shí)驗(yàn)測(cè)得該種材料的彈性模量E為729.60 MPa，泊松比ν為0.368，通過(guò)對(duì)比本次原位壓縮實(shí)驗(yàn)結(jié)果，可近似認(rèn)為在載荷2 之前材料均處于彈性階段。為了更加直觀地展示試件在壓縮載荷下的變形情況，在載荷1 和2 下，通過(guò)胡克定律彈性本構(gòu)關(guān)系，即可求得沿加載方向的彈性應(yīng)力場(chǎng)σzz在試件內(nèi)的分布［1］，見(jiàn)公式（2）。由于TATB 基PBX 材料本身是由尺寸、形狀各異的造型粉顆粒壓制而成，其結(jié)構(gòu)的高度非均勻性導(dǎo)致試件在單軸壓縮作用下的變形也是不均勻的。由于TATB 基PBX 材料含有隨機(jī)分布的張開(kāi)或閉合的初始損傷微裂紋，在材料只承受單向壓縮載荷時(shí)，由于其微結(jié)構(gòu)的局部非均勻性，試件中存在局部的拉伸應(yīng)力集中［1］。而TATB 基PBX 材料的力學(xué)性能對(duì)拉伸載荷更為敏感，拉伸開(kāi)裂是微裂紋激活和擴(kuò)展的主要形式［2］。為了揭示TATB 基PBX 材料的斷裂機(jī)制，根據(jù)DVC 算法獲得的應(yīng)變張量的6 個(gè)分量，即可采用胡克定律對(duì)試件的應(yīng)力分布進(jìn)行分析。如圖7 所示，沿載荷方向的應(yīng)力場(chǎng)（σzz）分布中有局部拉應(yīng)力集中現(xiàn)象，材料的損傷和破壞往往是由這些集中區(qū)域開(kāi)始。值得注意的是，對(duì)TATB 基PBX 材料內(nèi)部應(yīng)力場(chǎng)的分析，還需要考慮材料物理性能的各向異性，借助于多軸非線(xiàn)性本構(gòu)關(guān)系的發(fā)展，以更精確地揭示外載荷下TATB 基PBX 材料的變形非均勻性和破壞機(jī)理。

圖7 三維應(yīng)力場(chǎng)σzz（MPa）彈性階段Fig.7 Stress distribution of σzz（MPa）during elastic state

2.4 基于ZNCC 系數(shù)局部化的裂紋分析

受μ-CT 設(shè)備分辨率的限制，在對(duì)內(nèi)部裂紋的研究中，主要是觀測(cè)大于μ-CT 空間分辨率的裂隙。在裂紋萌生階段，由于寬度較小，往往肉眼不可見(jiàn)，容易被忽略。小于μ-CT 尺度的微裂紋難以從圖像中直接觀測(cè)，從而影響了對(duì)裂紋孕育發(fā)展全過(guò)程的揭示。微裂紋的存在會(huì)使得采集到的參考體圖像與變形體圖像中的灰度局部不守恒，通過(guò)DVC 方法分析獲得的全局應(yīng)變場(chǎng)分布中能夠測(cè)得變形局部化區(qū)域，進(jìn)而對(duì)低于μ-CT 分辨率的微裂紋進(jìn)行探測(cè)。通過(guò)應(yīng)用數(shù)字圖像相關(guān)方法的計(jì)算參數(shù)，如利用灰度殘余分布、相關(guān)性系數(shù)、局部應(yīng)變集中等參量，能夠很好地對(duì)材料內(nèi)部的損傷進(jìn)行量化表征，如已有研究中對(duì)顆粒材料的應(yīng)變場(chǎng)局部集中分布［19］，對(duì)奧氏體不銹鋼材料的灰度殘余局部集中分布的量化分析［20］，均觀測(cè)到了低于分辨率的微裂紋萌生，并對(duì)裂紋擴(kuò)展和加寬后的方向與形貌進(jìn)行了準(zhǔn)確探測(cè)。

通過(guò)對(duì)圖片進(jìn)行判別，發(fā)現(xiàn)在載荷4 下試件出現(xiàn)裂紋，裂紋處于試件邊緣位置，為便于理解，圖8a中給出了在slice100上裂紋的大概位置（圖8a 中黃框處）。為了研究DVC 方法在裂紋識(shí)別中的有效性，在裂紋尖端附近選擇一個(gè)大小為40 voxel×80 voxel×60 voxel 的計(jì)算區(qū)域（圖8a中紅框處）。該計(jì)算區(qū)域如圖8b 所示，x坐標(biāo)的范圍為［450，490］，y坐標(biāo)的范圍為［300，380］，z坐標(biāo)的范圍為［70，130］。采用DVC 分析時(shí)子 體 塊 尺 寸 為21 voxel×21 voxel×21 voxel，步 長(zhǎng) 為2 voxel，其他設(shè)置與前文一致。圖8c 給出了該計(jì)算區(qū)域內(nèi)的ZNCC 系數(shù)分布圖，由圖8c 可知，在裂紋附近，ZNCC 系數(shù)普遍偏低，這主要是由于在裂紋附近試件發(fā)生較大的變形，從而引起匹配失效。在CT 體圖像分辨率不足的位置，明顯可見(jiàn)ZNCC 系數(shù)局部化，即通過(guò)DVC 計(jì)算的ZNCC 系數(shù)分布圖，可以簡(jiǎn)單、高效地判定亞體素微裂紋位置在試件內(nèi)部的分布。

圖8 基于DVC 計(jì)算ZNCC 系數(shù)分布的內(nèi)部微裂紋分析：（a）邊緣開(kāi)裂的CT 灰度切片圖，（b）裂紋尖端的紅框區(qū)域體圖像，（c）紅框區(qū)域內(nèi)的ZNCC 系數(shù)分布圖Fig.8 Internal microcrack detection based on the distribution of ZNCC coefficient from DVC calculation：（a）a slice image of sample with a crack initiated from edge ，（b）CT image of crack tip region indicated by red box，（c）distribution of ZNCC coefficient of crack tip region indicated by red box

3 結(jié)論

在TATB 基PBX 材料單軸壓縮實(shí)驗(yàn)中，通過(guò)DVC方法分析不同載荷下試件三維體圖像，獲得了加載過(guò)程試件內(nèi)部變形場(chǎng)的演化情況。研究結(jié)果表明：

（1）在DVC 計(jì)算前對(duì)原始圖像進(jìn)行灰度直方圖拉伸，經(jīng)體圖像優(yōu)化，材料在4 種載荷下的ZNCC 系數(shù)均值在0.75 以上，DVC 計(jì)算結(jié)果具有較高的可靠性。

（2）壓縮過(guò)程中隨著載荷的增加，試件發(fā)生較大的壓縮變形，斷裂前達(dá)到了0.8%壓縮應(yīng)變。但是由于試件的顆粒體系的材料不均勻性，使其內(nèi)部三維變形場(chǎng)也呈現(xiàn)出一定的不均勻性。

（3）將ZNCC 系數(shù)分布圖應(yīng)用于試件內(nèi)部微裂紋的識(shí)別，能夠簡(jiǎn)單、高效地確定小于μ-CT 尺度的微裂紋分布情況。