基于CT成像技術(shù)的鹽漬土壤孔隙結(jié)構(gòu)識(shí)別與分析

洪明海，黃介生，曾文治，伍靖?jìng)ィ跹P(yáng)(.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 43007；.江蘇省鹽城市第三人民醫(yī)院，江蘇 鹽城 400)

0 引 言

土壤是植物根系從中獲取營(yíng)養(yǎng)物質(zhì)和水分的基礎(chǔ)，是植物賴以生存的根基，而土壤孔隙結(jié)構(gòu)和孔隙分布是影響土壤質(zhì)地、肥力、含水量等土壤性質(zhì)的重要因素。土壤孔隙是土體內(nèi)部的空隙，是土壤結(jié)構(gòu)的重要組成部分，土壤孔隙的數(shù)量、大小、分布及其連通性對(duì)土壤持水能力有著十分重要的作用[1]，影響著土壤有機(jī)污染物的吸附/解吸過(guò)程[2]。已有的研究多關(guān)注降雨[3-6]、土壤改良措施[7]、耕作方式、種植模式等對(duì)土壤孔隙的影響，很少考慮鹽分對(duì)土壤結(jié)構(gòu)的影響。然而，土壤鹽漬化已經(jīng)成為土地退化和影響作物生長(zhǎng)的主要威脅。位于華北的內(nèi)蒙古河套灌區(qū)，是我國(guó)最大的自流灌區(qū)，由于強(qiáng)烈的蒸發(fā)、較少的降水量和較淺的地下水位，約70%的耕地受到土壤鹽漬化的影響。因此，有必要就鹽分對(duì)土壤結(jié)構(gòu)的影響機(jī)理展開(kāi)研究，這對(duì)進(jìn)一步認(rèn)識(shí)土壤和植物的作用過(guò)程、改良和利用鹽漬土有著重要的意義。

要定量研究鹽分對(duì)土壤孔隙結(jié)構(gòu)的影響機(jī)理，必須獲取高精度的原狀土壤孔隙結(jié)構(gòu)并對(duì)其進(jìn)行預(yù)處理和提取分析。獲取土壤孔隙度最為經(jīng)典的方法是根據(jù)土壤容重和密度進(jìn)行計(jì)算[8]，但是這種方法的工作量大、耗時(shí)長(zhǎng)。另一種較為常用的獲取原狀土壤孔隙結(jié)構(gòu)的方法為序列數(shù)字圖像技術(shù)[9]，但該方法操作復(fù)雜并且獲取的切片圖像精度較低，對(duì)原狀土壤帶有一定的破壞性。近年來(lái)，CT成像技術(shù)[10]以其快速、準(zhǔn)確、便捷等優(yōu)點(diǎn)廣泛應(yīng)用于醫(yī)學(xué)、農(nóng)業(yè)、生物等多個(gè)領(lǐng)域，由于CT成像技術(shù)可以在不破壞土壤樣品的情況下對(duì)研究對(duì)象的內(nèi)部結(jié)構(gòu)進(jìn)行掃描和測(cè)定，且成像及分析速度快，操作過(guò)程相對(duì)簡(jiǎn)便，該技術(shù)已經(jīng)被逐漸應(yīng)用于土壤孔隙結(jié)構(gòu)的可視化和定量化研究[11]。本文通過(guò)CT掃描獲得內(nèi)蒙古河套灌區(qū)6個(gè)不同鹽分水平下原狀土壤的連續(xù)切面圖像，利用數(shù)字圖像處理技術(shù)并結(jié)合Arcgis 10.1軟件以及python程序批量提取連續(xù)圖像的土壤孔隙結(jié)構(gòu)參數(shù)，并對(duì)獲取的原狀土樣圖像和孔隙參數(shù)進(jìn)行對(duì)比分析，探求土壤鹽分對(duì)土壤孔隙的影響機(jī)制，為鹽漬土的改良和合理利用提供借鑒。

1 材料和方法

1.1 供試土壤及CT掃描

供試土樣取至內(nèi)蒙古自治區(qū)五原縣義長(zhǎng)試驗(yàn)站永聯(lián)基地(108°00′16.14″ E，41°04′11.57″ N)，海拔約為1 030 m，地處黃河北岸，為河套灌區(qū)(處106°20′ E～109°19′ E， 40°19′ N～41°18′ N)鹽漬土地帶。

首先利用土壤三參數(shù)傳感器(W.E.T Sensor, DELTA-T Co. Ltd)選擇6塊不同鹽分水平的田塊，并對(duì)每個(gè)田塊取樣帶回實(shí)驗(yàn)室用雷磁DDSJ-318電導(dǎo)率儀準(zhǔn)確測(cè)定土壤水溶液的EC1∶5值，最后通過(guò)公式[12,13]換算，以飽和土壤水浸提液的電導(dǎo)率ECe表示，各田塊鹽分值如下：田塊1(S1:ECe=1.51 dS/m)；田塊2(S2:ECe=4.33 dS/m)；田塊3(S3:ECe=15.33 dS/m)；田塊4(S4：ECe=21.31 dS/m)；田塊5(S5：ECe=28.89 dS/m)；田塊6(S6：ECe=54.92 dS/m)，其余田塊各層ECe值如表 1所示。

然后利用Φ50 mm×50 mm的PVC管在上述6個(gè)田塊分別采集原狀土樣，將取出土樣密封后用冷藏箱運(yùn)至江蘇省鹽城市第三人民醫(yī)院進(jìn)行CT掃描，所用CT機(jī)器型號(hào)為Discovery CT750 HD，CT掃描參數(shù)為：電壓140 kV，電流220 mA，層厚0.63 mm，設(shè)置窗寬(Window width)、窗位(Window level)為1×580，對(duì)每個(gè)土樣連續(xù)掃描80張CT圖片，圖像輸出DICOM格式醫(yī)學(xué)圖片[圖1(a)]。

從表1 可看出，除田塊5(S5)的20～30 mm的17.61 dS/m和田塊6(S6)的0～10 mm的68.97 dS/m兩個(gè)值比較異常外，其余各田塊各層ECe值都比較相近。而按照美國(guó)農(nóng)業(yè)部[14]根據(jù)飽和土壤水浸提液電導(dǎo)率ECe值對(duì)土壤鹽分水平的劃分，所取土樣鹽分水平也列于表 1，2個(gè)非鹽田塊(S1、S2)，1個(gè)中鹽田塊(S3)，3個(gè)高鹽田塊(S4、S5、S6)，無(wú)低鹽田塊。

表1 各田塊不同深度電導(dǎo)率值

1.2 圖像處理

由于考慮擾動(dòng)因素，去除土樣表層的5張和土樣底部的5張CT圖片，選擇余下的70張CT圖像作為此次分析：①需要利用專業(yè)的醫(yī)學(xué)圖像處理軟件Sante DICOM Viewer FREE軟件打開(kāi)CT掃描得到的DICOM醫(yī)學(xué)圖片，并進(jìn)行模式選擇和閾值調(diào)整處理，最終轉(zhuǎn)化為tiff格式保存[如圖1(b)所示]；②為了能使同一批圖片都能用Arcgis10.1地理軟件進(jìn)行處理，先用ArcCatalog 10.1對(duì)tiff格式的圖片構(gòu)建金字塔顯示(Buid Paramids)；③再利用Arcgis 10.1對(duì)圖片定義投影(Define Projection)；④使用裁切工具將tiff格式的柵格圖片剪切(Clip)為我們所要的半徑為20 mm的研究區(qū)域(以圖像中心為原點(diǎn))；⑤用柵格轉(zhuǎn)多邊形(Raster to Polygon)工具將校正后的圖片矢量化，得到如圖1(c)所示的圖片；⑥通過(guò)多邊形的屬性表提取出每一層孔隙的數(shù)量、單個(gè)孔隙的周長(zhǎng)、單個(gè)孔隙的面積等土壤孔隙結(jié)構(gòu)的基礎(chǔ)計(jì)算參數(shù)作為后續(xù)的分析使用。

圖1 土層切片中的孔隙分布情況

然而，在實(shí)際操作過(guò)程中，這一列的操作繁瑣、易出錯(cuò)、處理周期長(zhǎng)，且中間可能出現(xiàn)的問(wèn)題較多，這對(duì)于大量的田間試驗(yàn)的圖像處理是難以接受的；另外在矢量化后的圖片中往往會(huì)產(chǎn)生大量的小多邊形，這主要是因?yàn)镃T得到的土樣圖片中會(huì)有很多微小的小孔隙還有非二值化的一些圖像域，都會(huì)導(dǎo)致這種錯(cuò)誤，比如，某個(gè)土樣剖面的實(shí)際孔隙數(shù)只有52個(gè)，通過(guò)矢量化后得到的多達(dá)613個(gè)，這對(duì)于孔隙結(jié)構(gòu)參數(shù)的提取有很大的影響，精度會(huì)顯著下降。針對(duì)這些問(wèn)題，在圖像處理初期，用python編程統(tǒng)一將420張tiff圖片編號(hào)排序(如：“S1_1.tif”表示的是鹽分水平S1土樣的第一張CT掃描圖片的tiff格式)，并除去那些因?yàn)楹诎组撝刀鴮?dǎo)致的多邊形使計(jì)算結(jié)果更加準(zhǔn)確，由于Arcgis 10.1版本可外置python窗口處理，因此，前面所介紹的整個(gè)過(guò)程均通過(guò)python編程調(diào)用Arcgis 10.1中相應(yīng)的模塊實(shí)現(xiàn)批量處理，且后面的土壤孔隙結(jié)構(gòu)的基礎(chǔ)參數(shù)提取(見(jiàn)1.3)亦通過(guò)python編程實(shí)現(xiàn)并保存為“*.txt”文檔。

1.3 圖像分析

經(jīng)過(guò)圖像處理后提取到的孔隙面積、孔隙周長(zhǎng)，按照以下公式[9]計(jì)算土樣橫斷面孔隙的成圓率：

C=4 πA/P2

(1)

式中：C為成圓率，其值介于0～1之間，當(dāng)C=1時(shí)表示孔隙為圓形；A為孔隙面積，mm2；P為孔隙的周長(zhǎng)，mm。

統(tǒng)計(jì)土樣剖面孔隙面積的成圓率的最大值、最小值、平均值、標(biāo)準(zhǔn)差，分析它們?cè)谕寥乐猩疃鹊淖兓约霸诓煌}分水平之間的差異性。

孔隙數(shù)是土壤某層CT切片中孔隙的數(shù)目，記為n。

孔隙度的定義為單位容積土壤中孔隙容積所占的百分?jǐn)?shù)如下：

f=Vp/V

(2)

式中：f為土壤的孔隙度，無(wú)單位；Vp為土壤孔隙的體積，mm3；V為整個(gè)土壤的體積，mm3。

Vp=Sph

(3)

V=Sh

(4)

式中：Sp、S分別表示某一層孔隙的面積和某一層整個(gè)截面的面積，mm2。

將公式(3)、(4)代入式(2)中，消去相同的層厚h=0.63 mm，則有：

f=Sp/S

(5)

因此，如式(5)所示，計(jì)算土壤的孔隙度則是計(jì)算孔隙的面積占整個(gè)土樣剖面面積百分?jǐn)?shù)；分析的剖面取半徑為20 mm的圓形區(qū)域，因此整個(gè)剖面的面積為1 256.64 mm2。

2 結(jié)果與分析

2.1 土壤孔隙成圓率的變化

土壤孔隙的成圓率(C)是表征土壤孔隙的形態(tài)特征，成圓率值介于0～1之間，其值越接近于1說(shuō)明土壤的剖面的孔隙越接近圓形。一般認(rèn)為，成圓率受到土壤孔隙形狀和孔隙大小的影響，較小的孔隙其成圓率一般更接近1，因此根據(jù)成圓率可以判斷土壤孔隙的大小。因考慮表層和底部的擾動(dòng)影響，從表層和底部各去掉了5張切片，故實(shí)際分析范圍為3.15～47.25 mm。本次試驗(yàn)的6個(gè)土柱不同深度的成圓率分布如圖2所示。總體而言，6個(gè)土柱的成圓率值均介于0.5～0.75之間，并且沒(méi)有表現(xiàn)出隨土壤深度明顯的規(guī)律性變化，都在0.5～0.75之間波動(dòng)。

圖2 土壤孔隙成圓率隨深化度的變化

此外，鹽分影響土壤的成圓率。具體而言，S6土樣的鹽分含量最高(ECe=54.92 dS/m)其各層的成圓率也最大，并且成圓率的最大值0.742 5也在S6土樣中出現(xiàn)；S1土樣的鹽分含量最低(ECe=1.51 dS/m)，其平均成圓率也最小。此外，盡管中間鹽分水平的成圓率沒(méi)有呈現(xiàn)出明顯的規(guī)律性的變化，但總體而言，除S5土樣，土壤成圓率隨鹽分含量的變化(S2～S4)不斷增大，更接近S1。

2.2 土壤孔隙度的變化

圖3給出的是6個(gè)不同鹽分土柱的孔隙度f(wàn)隨土層深度的變化，其對(duì)應(yīng)的描述性統(tǒng)計(jì)分析參數(shù)和各層的平均孔隙度值如表2所示。具體而言，土壤孔隙度在土層深度上的變化是中間(約15～35 mm)比較大，兩頭值較小，而35 mm以后各處理的孔隙度也明顯小于上層土壤的，這在表2中也可以體現(xiàn)。此外，與成圓率規(guī)律類似，鹽分最高的S6處理(ECe=54.92 dS/m)孔隙度值最小(0.007)并且與其他鹽分處理的孔隙度差異明顯。同時(shí)，根據(jù)圖3和表2中給出的標(biāo)準(zhǔn)差0.005 9的值，發(fā)現(xiàn)，S6的波動(dòng)是最小的，S1處理(ECe=1.51 dS/m)的平均孔隙度最大且孔隙度值在整個(gè)深度范圍內(nèi)最不穩(wěn)定，波動(dòng)最大。除S5外，孔隙度隨鹽分的升高而表現(xiàn)出減小的趨勢(shì)，孔隙度的最大值(0.258 9)和最小值(0.000 1)也分別出現(xiàn)在S1處理和S6處理中。

圖3 土壤孔隙度隨深度的變化

表2 土壤孔隙度描述性統(tǒng)計(jì)

2.3 土壤孔隙數(shù)的變化

圖4給出的是土壤孔隙數(shù)隨深度的變化圖，從圖4中可以看出，土壤孔隙數(shù)在土層中的變化存在一定的差異，約在35 mm后都表現(xiàn)為下降的趨勢(shì)，深層(約40 mm以后)的土壤孔隙數(shù)明顯少于土壤表層和中間土層的土壤孔隙數(shù)，因此，從整體趨勢(shì)來(lái)看，隨著土層深度的增加孔隙數(shù)是有先增大后減小的趨勢(shì)，位于土層表面的孔隙數(shù)較土壤深層的孔隙數(shù)多，而中間(約15～35 mm)深度最多，各土層深度范圍內(nèi)的孔隙數(shù)變化都比較集中出現(xiàn)在20～60區(qū)間，并且變化差異較大，除S2(ECe=4.33 dS/m)處理出現(xiàn)81的高孔隙數(shù)值外，其他均未超過(guò)80的孔隙數(shù)值?？紤]不同鹽分的影響時(shí)，圖像最左側(cè)高鹽的S6(ECe=54.92 dS/m)處理的孔隙數(shù)最少，這與孔隙度和成圓率的變化規(guī)律都是一致的，因此可以認(rèn)為高鹽分的土壤孔隙較少，尤其不利于大孔隙的形成。

圖4 土壤孔隙數(shù)隨深度的變化

進(jìn)一步將土壤以層為單位分析每層土壤孔隙的數(shù)量以及整個(gè)土柱的孔隙數(shù)量，其結(jié)果如圖5所示。整體而言，S2(ECe=4.33 dS/m)處理的孔隙數(shù)量最多，為3 458，S4(ECe=21.31 dS/m)處理、S5(ECe=28.89 dS/m)處理的總孔隙數(shù)緊隨其后，S6(ECe=54.92 dS/m)的孔隙數(shù)最少[圖5(a)]。觀察圖5～(f)也可發(fā)現(xiàn)，表層土壤(3.15～10 mm)和深層土壤(40～47.25 mm)的孔隙數(shù)較中間的少。而不同土層的孔隙數(shù)量分析也表明，除了表層3.15～10 mm的S1土樣，其余各層的孔隙數(shù)量隨鹽分的變化規(guī)律均與全剖面的規(guī)律一致。而S1土樣3.15～10 mm土層中孔隙數(shù)較大可能是由于取樣受到擾動(dòng)所致。

3 結(jié) 語(yǔ)

(1)CT掃描技術(shù)能夠非破壞性地獲取土壤結(jié)構(gòu)特征，與其他方法相比，CT掃描方便快速并且能夠準(zhǔn)確描述孔隙的數(shù)目、大小等幾何參數(shù)。

(2)土壤孔隙結(jié)構(gòu)的形成受到鹽分的影響。具體而言，當(dāng)土壤鹽分較高S6(ECe=54.92 dS/m)時(shí)，土壤孔隙的成圓率明顯增加，土壤孔隙度明顯降低并且孔隙度的波動(dòng)隨鹽分增加逐漸減小，而鹽分最低的S1(ECe=1.51 dS/m)處理，其成圓率明顯較小、孔隙度較大。

(3)所取土柱的總孔隙數(shù)以及層間孔隙數(shù)隨著鹽分水平的提高表現(xiàn)為先增加后減小的趨勢(shì)，較高鹽分條件下(S4、S5)的土壤孔隙數(shù)較大，高鹽條件下(S6)的土壤孔隙數(shù)最少。

□

[1] 孫 梅, 黃運(yùn)湘, 孫 楠, 等. 農(nóng)田土壤孔隙及其影響因素研究進(jìn)展[J]. 土壤通報(bào), 2015,(1):233-238.

[2] 王震宇, 于曉冬, 許 穎, 等. 土壤微孔對(duì)有機(jī)物吸附/解吸的影響及其表征[J]. 生態(tài)學(xué)報(bào), 2009,(4):2 087-2 096.

[3] Marshall T J. A relation between permeability and size distribution of pores[J]. Journal of Soil Science, 1958,9(1):1-8.

[4] Belnap Jayne, Büdel B, Lange Otto L. Biological soil crusts: characteristics and distribution[M]. Springer, 2001.

[5] Fox D M, Bryan R B, Fox C A. Changes in pore characteristics with depth for structural crusts[J]. Geoderma, 2004,120(1):109-120.

[6] 郭 磊. 降雨過(guò)程中典型紅壤孔隙時(shí)空變化規(guī)律研究[D]. 武漢：華中農(nóng)業(yè)大學(xué), 2007.

[7] 楊永輝, 武繼承, 毛永萍, 等. 利用計(jì)算機(jī)斷層掃描技術(shù)研究土壤改良措施下土壤孔隙[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2013,(23):99-108.

[8] 程亞南, 劉建立, 張佳寶. 土壤孔隙結(jié)構(gòu)定量化研究進(jìn)展[J]. 土壤通報(bào), 2012,(4):988-994.

[9] 李德成, Velde B, 張?zhí)伊? 利用土壤的序列數(shù)字圖像技術(shù)研究孔隙小尺度特征[J]. 土壤學(xué)報(bào), 2003,(4):524-528.

[10] 周 虎, 李文昭, 張中彬, 等. 利用X射線CT研究多尺度土壤結(jié)構(gòu)[J]. 土壤學(xué)報(bào), 2013,(6):1 226-1 230.

[11] Munkholm Lars J, Heck Richard J, Deen Bill. Soil pore characteristics assessed from X-ray micro-CT derived images and correlations to soil friability[J]. Geoderma, 2012,181:22-29.

[12] Hachicha M, Mansour M, Rejeb S, et al. Applied research for the utilization of Brackish/Saline water in center of Tunisia: water use[C]∥ Salinity evolution and crop response: Proceedings of International Salinity Forum, Riverside, 2005.

[13] Kanzari Sabri, Hachicha Mohamed, Bouhlila Rachida, et al. Simulation of water and salts dynamics in bouhajla (Central Tunisia): Exceptional rainfall effect[J]. Soil and Water Research, 2012,7(1):36-44.

[14] Department Natural Resources Management And. CHAPTER 7-SALTY SOILS[EB/OL]. http:∥www.fao.org/docrep/R4082E/r4082e08.htm.,2016-04-17.