所用土壤选自黑龙江省克山农场, 土壤类型以黏化湿润均腐土为主, 采用自制内径和高分别为10 cm的有机玻璃管于0~40 cm层深的侵蚀沟壁进行原状土取样, 每层重复取样3次。将样本进行饱和水处理后^[7], 共得到原状饱和水(3个)、原状非饱和水(3个)、填充饱和水(3个)、填充非饱和水(3个)等共计12个样本, 以供后续CT扫描所用。

选用黑龙江省中医药大学的Lightspeed16排螺旋CT扫描仪对原状样本进行扫描处理。CT扫描仪选用参数分别为:电压140 kV, 电流60 mA, 扫描间隔1.5 s, 扫描层厚2 mm, 窗宽(显示CT图像是所选用的CT值范围)和窗位(窗宽上、下限CT值的平均数) 均为1 300, 对每个土柱样品进行7次横断面扫描。

本研究的单个扫描样本分别经历0、1、3、6、9次冻融循环, 每次冻融循环后进行1次扫描, 每次扫描可得7张土壤CT图像, 故单个样本冻融循环后共需5次扫描, 得35幅土壤CT图像。本研究合计获得420幅土壤CT断层扫描图像, 用于建立土壤CT图像数据库, 作为测试图像库。

1.3.1 土壤CT图像处理

原始土壤CT图像以医疗图像的格式存储, 并不便于计算机的后续处理。同时, 由于CT机器硬件技术缺陷的存在, 使得图像边界会产生畸变, 而且, 数字图像在获取和传输的过程中也会受到噪声影响, 这就决定了基于原始图像的特征提取可能会影响试验结果的精确性^{[8, 9]}。因此, 须针对土壤CT图像首先进行一些处理。

如图1-a所示的原始图像中包含了与土壤分析无关的CT机器扫描参数的基本信息, 因此, 本试验基于Matlab编程软件, 采用最大内切正方形的方法, 剪裁得到有效土壤CT图像, 并保存为* .bmp格式文件。通过自适应中值滤波算法消除噪声对后续特征提取的干扰, 以此来抑制无用信息, 突出有效信息, 得到用于结果分析的试验用图像, 如图1-b所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 土壤CT图像预处理Fig.1 Preprocessing of soil CT image

1.3.2 灰度共生矩阵法

纹理是土壤CT图像中每个像素的CT值或灰度值在空间中以某种方式变化而产生的, 是图像客观存在的一种属性, 反映了图像中灰度的空间相关特性^[10]。灰度共生矩阵法是通过研究灰度的空间相关特性来描述图像纹理变化的方法, 反映图像灰度关于方向、相邻间隔、变化幅度的综合信息, 是分析图像的局部模式和排列规则的基础。

灰度共生矩阵定义为:灰度值为i的像素沿θ 方向, 到达与其相隔距离为d的像素点(灰度值为j)的概率。设g(x, y)为图像中位于(x, y)像素点的灰度值, 则其数学表达式为:

P(i, j, d, θ )=#{(x₁, y₁), (x₂, y₂)|f(x₁, y₁)=i, f(x₂, y₂)=j, |(x₁, y₁)-(x₂, y₂)|=d, ∠(x₁, y₁), (x₂, y₂)=θ }。(1)

式(1)中, #表示根据式(1)的条件, 图像中位于(x₁, y₁)灰度值是i, 位于(x₂, y₂)灰度值是j的像素对个数。通常, d={1, 2, 3, 4}, θ ={0° , 45° , 90° , 135° }。

为了便于二次统计特征的提取, 需对灰度共生矩阵进行归一化处理:

P(i, j, d, θ )= P(i, j, d, θ)R,

R= N(N-1) θ=0°或θ=90°(N-1)2θ=45°或θ=135°。(2)

Haralick等^[11]根据纹理信息的特点, 提出了用于分析灰度共生矩阵的多个特征, 用于挖掘图像深层信息。常用的特征参数简述如下。

(1)角二阶矩

式(3)中, k为灰度级[式(4)~(6)同]。角二阶矩也称为能量, 是灰度共生矩阵各个元素值的平方和, 反映了图像灰度分布均匀程度和纹理粗细度。

(2)对比度

对比度, 也称主对角线惯性矩, 代表检测像素灰度值与其领域像素灰度值的亮度反差, 反映了图像的清晰度和纹理沟纹深浅的程度。

(3)逆差矩

逆差矩用于度量图像纹理局部变化的多少, 反映图像纹理的同质性, 其值大说明图像纹理的不同区域间缺少变化, 局部非常均匀。

(4)熵

熵是图像所具有的信息量的随机性度量, 反映了图像中纹理的非均匀程度或复杂程度。

1.3.3 主成分分析法

主成分分析是通过线性变换将多个指标简化为少数指标的统计分析方法, 用少数互不相关的指标尽可能多地表征图像信息, 以保证原始信息损失量最小^[12]。其基本原理是:将原始具有相关性的多个指标X₁, X₂, …, X_p(p个指标), 重新组合成一组互不相关的综合指标F_m来代替原来指标。

设X=(X₁, X₂, …, X_p)^T是p维随机向量, 其线性变化为:

F₁=a₁₁X₁+a₁₂X₂+… +a_1pX_p; (7)

采用新变量F₁来表征原始p个向量, F₁所包含的信息应尽可能多地反映所有变量信息, 则其为第一主成分。若第一主成分不足以代表足够多的原始信息, 则引入第二主成分, 以此类推。其中, 新变量F₁, F₂, …, F_m要求满足以下条件:

(1)变量所提取的信息量由方差度量, 即要求方差Var(F_m)尽量大;

(2)F_i与F_j互不相关, 即COV(F_i, F_j)=0。

主成分分析的主要目的是简化数据量, 因此, 所提取的主成分个数m应小于p。m的值依据各主成分的累计方差贡献率P(m)来确定。其数学表达式为:

P(m)=∑i=1mλi/∑k=1pλk。(8)

式(8)中:λ 为各个主成分对应的特征值; k为全部的主成分个数; i为选定主成分个数。

一般地, 当累积贡献率大于85%时, 就认为足够反映原来变量的信息^[13], 也可根据实际应用的需求自主选择累积贡献率的范围。

1.4.1 验证方法确定

由于土壤CT扫描图像数据库中单个样本单次冻融循环的图像为7幅, 数据量相对较少, 需选用适合于小样本的验证方法。留一验证法的思想是选择一幅作为测试图像, 剩余的作为训练图像, 以此循环, 直至所有图像都参与测试与训练。该方法对图像的利用率高, 能提高实验结果的精确度, 故本研究选用该方法用于对土壤CT图像的判别与分类^[14]。

1.4.2 欧氏距离计算

欧氏距离指的是在m维空间中2个点之间的真实距离, 通过欧氏距离的数值可准确定义2个特征向量间的相似度。

1.4.3 识别正确率计算

若测试图像对于本组内欧氏距离较其他组小, 则为判别正确, 反之则为错误, 最终正确判断个数与总判断数的比值即为识别正确率。

2.1.1 灰度共生矩阵特征

灰度共生矩阵法的主要变量是方向和特征个数。选定0° 、45° 、90° 、135° , 提取能量、熵、对比度、逆差矩共4个特征对土壤CT图像纹理信息进行描述。由于篇幅有限, 随机选取经历不同冻融次数的同一土壤CT图像进行部分特征数据的展示, 结果如图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 灰度共生矩阵特征值Fig.2 The eigenvalue of gray symbiotic matrix

能量和逆差距主要反映土壤纹理信息的同质性。由图2可知, 不同方向的能量和逆差距先随着冻融次数的增加而增加, 当冻融次数达到6次后趋于减小, 说明土壤内部结构在前6次冻融过程中变化比较强烈, 之后变化减弱, 内部物质结构渐趋均匀。两个方向的熵和对比度数值较小, 且呈现相反的变化趋势, 这说明当土壤图像纹理信息趋于随机时, 土壤内部物质之间具有一定相似性。4个特征综合起来表征了土壤内部物质随着冻融次数的变化, 在第9次冻融循环后内部物质趋于均匀化。

2.1.2 主成分分析特征

采用主成分分析方法, 对不同冻融循环次数的土壤CT图像进行主要特征分量的提取。选用4种冻融循环次数下的图像, 分别记为A₁、A₂、A₃、A₄, 设定各主成分的累计方差贡献率大于0.85, 最多选取16维主成分特征。由于篇幅有限, 选取4幅图像的10维特征的贡献率进行展示, 如表1所示, a_i表示提取的第i维特征。

表1

Table 1

表1(Table 1)

表1 主成分特征贡献率 Table 1 Contribution rate of principal component 图像Image a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 A1 0.1915 0.1787 0.0923 0.0519 0.0463 0.0439 0.0419 0.0329 0.0324 0.0276 A3 0.2647 0.1980 0.0631 0.0483 0.0406 0.0386 0.0312 0.0278 0.0272 0.0242 A6 0.2732 0.1629 0.0889 0.0757 0.0663 0.0454 0.0336 0.0316 0.0258 0.0221 A9 0.2922 0.1683 0.0929 0.0651 0.0516 0.0364 0.0316 0.0272 0.0255 0.0248

表1 主成分特征贡献率 Table 1 Contribution rate of principal component

由于图像所包含的信息不同, 图像提取的每维主成分贡献率亦不相同, 累计方差贡献率大于0.85的主成分特征维数也不相同, 最大有16维, 最小12维。为保证图像识别的正确率, 所有图像均选取16维主成分特征参与后续的计算分析。

由于篇幅有限, 仅举一例用于说明判别过程。首先, 将原状饱和水样本经历1、3、6、9次冻融循环的扫描图像分别编组, 标记为Y₁、Y₃、Y₆、Y₉, 下标中的数字代表其所经历的冻融次数; 然后, 选择Y₃组中的7幅图像(即经历3次冻融的原状土图像), 依次命名为I₁~I₇, 作为测试样本, 分别测算其每一幅图片的灰度共生矩阵特征向量与Y₁、Y₃、Y₆、Y₉组均值特征向量之间的欧式距离, 与哪一组的欧式距离数值最小, 则自动判别其属于该组。如表2所示, 属于Y₃组的7幅图像中, 除I₄被自动判别为Y9组(即经历9次冻融), 判别错误外, 其他判断均正确, 则该测试样本的判别正确率为(6÷ 7)× 100%=85.7%。

表2

Table 2

表2(Table 2)

表2 样本欧氏距离 Table 2 The Euclidean distance of samples 组别Group I1 I2 I3 I4 I5 I6 I7 Y1 9.9794 6.6998 8.5730 9.2161 11.4059 9.1988 9.5257 Y3 6.0742* 5.0340* 5.6900* 6.6188 5.7590* 6.9309* 6.0481* Y6 10.7891 6.7512 9.3236 10.4523 11.1434 9.9628 11.0759 Y9 6.0932 5.2856 5.7721 6.0670* 6.7939 7.0306 6.0547 * indicated the smallest Euclidean distance within the image tested in the column and the corresponding groups. 每列中加* 的表示所在列图像与该组的欧氏距离最小。

表2 样本欧氏距离 Table 2 The Euclidean distance of samples

基于本研究所用土壤CT图像数据库, 对全部样本进行自动判别, 最后测定得出采用灰度共生矩阵法自动判别的正确率是81.4%, 主成分分析法自动判别的正确率是61.0%。