生物扰动作用(bioturbation)是指底栖生物的摄食、爬行、排泄和驻穴等生命活动对沉积物初级结构发生的改变, 通常伴随着沉积物的转移、混合以及上覆水体涌入穴道的现象(Meysman et al, 2006; 覃雪波, 2010; Kristensen et al, 2012)。生物扰动构造是表征海底生态环境健康的重要指标, 生物扰动作用促进沉积物与介质之间的物质交换, 并可引起沉积物的扩散混合从而影响沉积物中记录的古环境信号(Lyttle, 2013)。当前沉积记录的研究分辨率可达百年尺度和年代际尺度, 由于生物扰动作用造成的测年偏移会使岩心的年代框架产生错乱, 最终影响到整个研究结果(杨群慧等, 2008)。所以, 在进行高分辨率沉积记录研究时, 应尽可能避免生物扰动作用对同位素测年的影响或者寻求其他手段消除这个影响(范德江等, 2008)。正因如此, 海洋研究者越来越重视生物扰动作用机制的研究(Mermillod-Blondin, 2011)。

受到观测手段和研究方法的限制, 现代海洋沉积中生物扰动构造的研究一直处在定性-半定量化的阶段, 这是由于生物扰动作用存在不确定性(龚一鸣等, 2009; 张翔宇, 2018), 以及扰动构造后期改造严重所致。21世纪之后数字图像处理技术得到了极大的发展, 在地质学研究领域取得了一系列优秀成果。计算机断层成像技术是X射线技术的延伸, 通过螺旋扫射的方式, 可以获得样品三维立体的密度信息, 同时还继承了X光负片无损性的特点(Hamblin, 1962)。故本文尝试利用计算机断层成像技术实现生物扰动构造的三维重构, 完成生物扰动强度的定量表征, 并对长江口现代海洋沉积物中的生物扰动构造进行定性和定量研究, 初步分析生物扰动构造的发育特征和影响因素, 为生物扰动作用机制以及环境效应研究提供基础。

1 样品采集及分析方法 1.1 样品采集

沉积物样品于2020年10月由向阳红18号科学考察船执行东海共享航次时获得。通过箱式取样器获得底质沉积物, 然后进行插管取样, 共计获得4根沉积物柱状样。其中, SF-1站位处在长江水下三角洲前缘斜坡处, 受长江冲淡水作用明显。S01-1和S01-2站位在前三角洲边界地带, S02-3站位进一步向海延伸, 处在过渡带之上。柱状样品采用长60 cm、横截面内径直径7.5 cm的聚氯乙烯(polyvinyl chloride, PVC)管采集。站位信息见表 1, 站位分布如图 1所示。

1.2 分析方法 1.2.1 沉积物岩心CT图像获取

插管岩心样品CT (computed tomography)数据是使用美国通用公司64排128层螺旋GE OPTIMA 660CT仪器通过平扫的方式获取的。经多次尝试使用120 kVp、50 mA的工作环境能够使图像质量达到最高水准。本次扫描分辨率是0.625 mm, 即每隔0.625 mm获取一次切片的CT数据; 空间分辨率是512×512, 即最后成图的水平像素点和垂直像素点个数均为512个(图 2)。最终将图片以Dicom格式保存, 供后期进行数字图像处理。4根岩心共获得2 692张切片。

1.2.2 沉积物粒度和含水率测定

沉积物粒度测定: 选取SF-1、S02-3、S01-2沉积物岩心进行分样, 取样间隔0.25 cm, 共计获得442个沉积物样品。使用英国Malvern公司生产的Mastersizer3000型激光粒度仪进行粒度分析, 具体步骤如下: 取约1 g的样品放入50 ml离心管当中, 加入5 ml浓度为30%的双氧水用以去除沉积物中的有机质, 摇晃均匀, 此时会发现气泡反应, 静置24 h直至气泡消失, 如若气泡不消失则继续静置。然后加入3 ml的0.5 mol/L的六偏磷酸钠, 使沉积物颗粒分散, 最后超声30 min后进行上机测试。粒度分级采用Ф值标准, 粒度数据以0.25Ф的采集间隔导出, 粒度参数采用McManus矩法公式计算。

含水率测定: 选取各个岩心中分布沉积物样品进行含水量测定, 将样品放在电子恒温干燥箱中于110 ℃干燥至恒质量, 用公式(1)计算含水率。

其中, ω为含水率(%), W₁为干燥前湿样质量, W₂为干燥后样品质量。

2 计算机断层扫描的三维重构技术 2.1 基于CT值的数字图像处理技术

CT值(记为V_CT)的标准定义是具等厚的物质对X射线衰减的影响变化, 它是计算机断层扫描中普遍使用的无量纲单位。计算公式如下(汪家旺等, 2004; 彭文献等, 2010):

式中, μ为衰减系数, a为标度因素, 当a=1 000时, V_CT为Hounsfield单位。

前人研究指出物质衰减系数与物质密度成线性正相关(章程辉等, 2006), 因此推断CT值受到物质密度的控制, 可以作为密度大小的反映。在沉积学的研究当中, 沉积物密度是一个非常重要的参数, 受到含水率、孔隙度、组分含量等因素影响。同时CT成像技术相较于传统X光技术可以获得高分辨率的信息。对基于CT值的数字图像处理技术包括图像预处理、原始CT值提取等工作流程。

2.1.1 图像预处理

利用ImageJ软件打开沉积物岩心CT图像的Tag信息列表, Tag通过数据集的形式保存了拍摄时间、工作环境、扫描分辨率、空间分辨率、窗宽数值、三维空间位置等信息, 这对后续处理工作非常重要。获取Tag信息之后就需要将原始CT切片进行预处理工作, 选用Python编程语言作为工具开展预处理工作, 包含图片排序和图片裁剪两个步骤。

(1) 图片排序: 医用级CT仪器对X光切片的排列顺序非常混乱, 与我们常用的排列顺序逻辑不一致, 不同的设备之间的排序方法也不一致, 有时会出现切片缺失的情况, 对于生物扰动构造的识别造成不便, 因此需要根据实际层厚(0018, 0050 Thinkness)和切片实际位置(0020, 1041 Slice Location)进行重新排序。

(2) 图片裁剪: 因沉积物柱状样CT图像拍摄时为多根一组同时扫描, 故需要将CT切片进行批量裁剪, 并且在裁剪过程中需要保持坐标系统不变, 这里运用到Python的os、mumpy和matplotlib三个工具包, 在可视化页面下将图片裁剪为矩形。最后将矩阵数据保存成mat格式, 以便后期在Matlab环境下进行操作。

2.1.2 原始CT值的获取

CT值会被计算机通过直线映射的方式转换为灰度值, 故此我们需要通过逆映射的方式获取CT值, 逆映射公式如下所示:

式中, p为像素pixel的取值; s为斜率slope的取值; i为截距intercept的取值。

每个点的像素可以在Python中通过get函数获得, 斜率和截距均能在Tag信息中获取。至此原灰度图像已经转换为CT值的矩阵, 将所获得的矩阵保存为mat格式。

2.2 CT值与沉积物平均粒径和含水率的关系

在CT切片图像上取40×40像素点大小的矩阵, 获取每一个像素点的CT值, 利用Matlab中的mean2函数求取二维矩阵的平均值, 这样就获得每一个切片的平均CT值, 利用沉积物平均粒径和含水率的高分辨率记录, 分析CT值与平均粒径和含水量的关系。选用S02-3站位沉积物数据, 绘制CT平均值、含水率、平均粒径变化趋势图(图 3), 提取相同深度的CT平均值与平均粒径和含水率分别进行统计学分析(图 4)。

如图 3和图 4所示, CT值与含水率有较为显著的线性负相关关系, 符合线性回归方程Y=1.035 6- 0.000 5X, 系数在99.9%的置信区间上显著, 拟合优度R²=0.551, 2.5 cm处含水率较低而该处CT值显著升高, 5 cm处含水率较高CT值则非常低, 整个沉积物岩心20 cm以上区域的含水率约在50%以上, 20 cm以下的区域含水率较低, 约为40%, CT值在20 cm也出现了变化, 底部区域CT值显著高于顶部区域。

CT值与沉积物平均粒径的波动趋势也有着一定的对应关系。经过统计学分析, 发现平均粒径与CT平均值符合二次项回归方程Y=303.9–0.060 87X+ 0.000 3X², 所有系数均在99.5%的置信区间上显著, 拟合优度R²=0.15, 二次项回归方程的对称轴在CT平均值为1 015处, 即当CT平均值小于1 015时, CT平均值上升而平均粒径变小; 当CT平均值大于1 015时, CT平均值升高则平均粒径变大。选用数据缩尾法对回归结果进行稳健性检验, 分别去除5%的极值数据后重新进行回归分析, 发现平均粒径与CT平均值仍然符合二次项关系, 所有系数在99.5%的置信区间上显著, 拟合优度R²=0.17, 拟合优度的上升说明回归分析通过稳健性检验。从地质学角度来讲, 浅表层区域的沉积物较为松散, 固结程度差, 平均粒径大引起沉积物孔隙度变大, 导致密度降低。20 cm以下的深部区域, 沉积物固结程度相对较高, 孔隙度对于沉积物密度的影响减小, 此时CT值与平均粒径呈现微弱的正相关关系, 这也揭示了该站位样品由松散到紧密的压实作用。

S02-3站位沉积物CT值具有显著的分段性, 20 cm以上部分CT值较低且变化幅度较大, 20 cm以下深度CT值较高且变化幅度不大。结合平均粒径和含水率的变化状况不难看出, 沉积物CT值对于含水率和平均粒径具有一定的指示意义, 可以很好地展现沉积物密度变化。

2.3 CT图像三维重构的实现

三维重构技术是指将多组二维图片按照一定顺序输入计算机, 得到几何物体的各项参数信息后重构其几何状态(王宗彦等, 2002)。针对连续切片的CT图像常采用体重构法(volume rendering), 即将CT获得的连续切片按照比例放大, 再按照原层厚组合起来, 实现三维空间的重构。该技术已经广泛应用于医学、材料学等领域, 并且有多种软件可以实现三维重构, 例如3D Slicer、3D Doctor等。本文利用了3D Slicer软件用以沉积物CT图像三维重建。

3D Slicer是一款开源的医学图像处理软件, 基于VTK、Python、Qt的基础上研发, 支持Windows系统, 提供了多模式数据处理的可视化操作功能(曾文晔, 2013)。基于3D Slicer的沉积物CT图像三维重构工作流程如图 5所示。利用该方法对沉积物样品的CT图像进行三维重构, 获得4个重构模型。

2.4 生物扰动强度的定量表征

生物扰动强度通常以扰动比例来表示(王英国等, 1999), 即生物扰动构造体积与沉积物体积的比值, 数学表达式如下:

式中, V_扰动表示生物扰动构造体积; V_沉积物表示沉积物体积。本文为实现生物扰动强度的定量化表达, 在高分辨率的条件下将扰动构造近似看做规则圆柱体, 生物扰动比例的计算公式可以得到如下变换:

其中, i代表切片数量; k_i表示第i个切片的生物扰动强度; S_i^b代表第i个切片上生物扰动的面积; S_i^s表示第i个切片上沉积物的面积; h表示圆柱体高度即切片层厚。因为沉积物的面积和切片层厚是固定值, 故全样沉积物生物扰动强度可以进行式(6)的变化:

即全样沉积物生物扰动强度为各切片生物扰动强度的平均数。通过式(6)就可以将生物扰动强度的体积比转换为面积比。下面就面积圈定法的工作步骤进行阐述。

前文已将沉积物CT图像转换为矩阵数据进行保存, 利用Matlab进行下一步分析处理, 处理之前结合生物扰动构造的三维重构图像, 本文统计了生物扰动构造特征要素的CT值, 如表 2所示。

图像增强是数字图像处理的关键步骤, 其含义是根据矩阵数据进行数字变换使想要研究的区域特点放大, 改善图像的视觉效果(解睿, 2019)。常采用阈值调节的手段, 本文曾尝试多种算法调节阈值凸显生物扰动构造特征, 同时减小沉积物、贝壳、含水率的影响, 但是结果并不如人意, 这是由生物扰动构造的随机性和环境的复杂性所决定的。故根据特征要素的CT值统计结果, 划分CT值梯度用以强化生物扰动构造的展示效果。梯度选择如下: (-1 000, -500, 0, 100, 150, 200, 250, 300, 350, 400, 450, 500, 600, 650, 700, 900, 1 300, 1 500)。在Matlab中加载预处理完成后的数据矩阵, 根据CT值梯度画横剖面和纵剖面的等值线填充图像, 如图 6和图 7所示。共计获取横剖面实验用图2692张, 纵剖面实验用图678张。

通过图 6和图 7可以很清晰地观察到生物扰动的虫孔构造。结合虫孔构造的生长性特征, 并用三维重构图像和纵剖面图像可以将生物扰动构造与沉积物区分开, 之后使用ImageJ软件通过标定感兴趣区域获取生物扰动面积。根据计算公式即能实现生物扰动强度的定量化表达。

3 长江水下三角洲的生物扰动构造 3.1 长江水下三角洲概况

长江水下三角洲是由三角洲平原、三角洲前缘、前三角洲组成, 呈从河口开始舌形向东南突出(陈吉余等, 1959, 1987; 陈中原等, 1986; 董永发, 1989)。该处近海环流发育, 水动力环境不稳定, 以长江冲淡水的影响最为剧烈, 它具有显著的季节变化特征(王亮, 2014)。长江入海陆源碎屑物质大约有50%在此沉积, 以细砂粉砂为主, 底质沉积物类型复杂。长江口及其邻近海域存在大型底栖生物约400余种, 主要分为四大类: 多毛类、软体动物、甲壳类、棘皮动物(刘录三, 2002)。多毛类和甲壳类整体分布广, 遍布整个长江水下三角洲区域, 软体动物主要分布在三角洲前缘斜坡和泥质沉积区当中, 棘皮动物种类较少, 多分布于浙闽沿岸和陆架斜坡处。多毛类是在生物量上占据绝对优势, 自近岸至外海逐渐降低, 是该区最主要的造迹生物。

3.2 生物扰动构造的三维重构

利用三维重构技术对各站位生物扰动构造进行重构, 可以清晰地显示各个站位沉积物中生物扰动构造发育和分布特征(图 8)。

(1) SF-1: 生物扰动构造数量较多, 各个层深均有分布, 整体扰动强烈。表层存在单侧通水的逃逸虫孔构造, 该构造长约5 cm, 垂切于沉积物表面, 其下部发育有一巨大的U型虫孔构造, 斜插在PVC管中, 该构造属于觅食虫孔构造, 视野内直线长约18 cm, 虫孔内径最大处约为1.2 cm, 虫孔周围的沉积物扰动明显, 影响到虫孔壁的观察, U型虫孔构造近沉积物界面一侧似与顶部的逃逸虫孔构造连接在一起, 形成共生组合。中下部多发育各种形态的觅食虫孔构造, 亦可见部分残存的挖掘构造, 最长一处构造痕迹长约12 cm, 内径不足0.1 cm。

(2) S01-1: 该站位生物扰动构造集中在浅表层5—10 cm处, 顶部生物扰动发育, 中下部生物扰动几乎不发育。单个扰动构造呈现U型、Y型等多种形态, 多数为逃逸虫孔构造, 垂向长度小于10 cm, 内径细小, 纤维管状的虫孔构造相互交缠, 呈树根状错综复杂的组合形态。

(3) S01-2: 该岩心生物扰动状态与S01-1相似, 岩心上部(0—10 cm)逃逸虫孔构造发育, 垂向长度约为7—8 cm, 相互交缠呈网状; 中部区域可见较多地觅食虫孔构造, 孔径较顶部的逃逸虫孔构造明显增大, 约为0.5 cm, 构造周围有不规则沉积物聚集, 推测为虫孔通水后致使周围沉积物密度发生改变所致, 顶部的逃逸虫孔构造与中部觅食虫孔应为共生组合; 下部生物扰动不发育。

(4) S02-3: 整体上来看该岩心生物扰动作用较为微弱, 且扰动构造主要集中发育在顶部, 以微小纤维管形式出现, 向下生物扰动构造零星出现; 底部出现一生物遗迹构造, 呈现月牙形, 由许多球形颗粒堆积而成, 推测为生物排泄物的堆积体。多数虫孔构造都有重新充填的痕迹, 与周围沉积物发生混合作用, 密度特征效应降低, 导致形态展示模糊。

3.3 生物扰动强度及其垂向分布

按照5张切片厚度(3.125 mm)的分辨率, 利用面积圈定法对四根沉积物岩心进行生物扰动强度的计算, 根据深度绘制垂向分布图像。

岩心SF-1、S01-1、S01-2、S02-3的平均生物扰动强度分别为1.325%、0.889%、1.680%、0.471%, 扰动强度的极大值点出现在S01-2岩心的表层0.3 cm处, 约为10%, 极小值为0, 表示无扰动状态, 多出现在岩心中下部, 四个不同岩心的生物扰动强度差异显著。同时扰动强度的垂向分布特征也不尽相同(图 9), S01-1和S01-2站位生物扰动强度垂向上呈现向下衰减的形式, 但是衰减速度不一致, S01-1站位扰动强度的衰减速度较慢, 呈现垂向递减形式。S01-2站位扰动强度的衰减速度较快, 呈现指数衰减的形式。SF-1站位沉积物中部生物扰动强度较高而两侧较低, 呈现形式为中间隆起的橄榄形式。S02-3站位的生物扰动强度微弱, 垂向上呈现波动减少的特点。

3.4 生物扰动构造发育影响因素的初步分析

刘卫东等(2009)在进行生物扰动构造研究的过程中曾经指出温度、光照、食物是影响底栖生物活动的重要因素, 而王慧中(1985)指出沉积速率较快的河口地区不适合生物扰动构造的发育。本文以底质沉积物的类型和上覆水的温盐数据为支撑, 对生物扰动构造发育的影响因素进行初步分析。

3.4.1 底质类型

底质类型会影响大型底栖生物的栖息环境, 并且朝利于生物活动的环境迁移。生物扰动构造在浅表层发育显著, 故取顶部(10 cm以上)沉积物生物扰动强度与沉积物组分含量进行分析, 将生物扰动强度和沉积物组分含量汇总见表 3。

结合图 9和表 3可知, 生物扰动强度与底质中黏土含量密切相关, 黏土含量高的站位(如SF-1)生物扰动强度较低, 而黏土含量低的站位(如S01-2)生物扰动强度较高。黏土组分在生物扰动作用中代表的主要是颗粒细、孔隙度小、固结程度强的沉积物, 黏土组分含量低则砂质+粉砂质含量高, 含水率和含氧量均会随砂质+粉砂质的含量上升而升高, 可能形成利于生物栖息、觅食的环境, 促进生物扰动作用。

从生物扰动构造类型的角度来将, 砂+粉砂质含量较高的沉积物中易发育网状、根系状虫孔构造, 虫孔交织缠绕, 对沉积物扰动强烈。而黏土组分含量高的沉积物中生物扰动构造发育程度较弱, 主要发育挖掘构造, 底栖生物向深部掘穴获得生命所需营养物质、空气和水。沉积物粒级组成与生物扰动强度具有显著的相关性(覃雪波, 2010)。在加拿大芬迪湾的生物扰动作用研究中Dashtgard等(2008)建立粒度参数与总有机碳的关系以及野外实测调查解释这一现象, 同时也指出相同环境下生物扰动易在砂质含量高处发生。本文实测结果对该结论加以印证。

3.4.2 水团和沉积速率的影响

长江水下三角洲上覆水团受河流和海洋的混合作用, 温度、盐度变化较大, 在河口汊道区以淡水为主, 水流湍急, 而前三角洲和陆架过渡区则以海水为主, 受环流、冲淡水的共同作用, 不同的水团环境所发育的生物扰动构造和扰动强度也不尽相同。相对而言, 在受长江冲淡水影响的海域, 因为接收了丰富的长江来源的营养盐(董书航, 2015), 有利于生物的生存, 可以促进底栖生物的发育, 导致该区S01-1和S01-2的生物扰动强度值较高, 生物扰动类型以网状虫孔构造和根系状虫孔构造为主, 集中发育在顶部沉积物区域; 而S02-3站位于正常陆架, 受长江来源影响小, 生物扰动作用较前两处站位微弱, 接近于不发育状态。同样受长江冲淡水强烈影响的SF-1站, 总体上生物扰动强度反而较低, 或是因为高速的沉积速率不利用大型底栖生物生活。

统计前人在该研究区域²¹⁰Pb、¹³⁷Cs测年的研究结果, 通过插值法获得四个站位的沉积速率。SF-1站的沉积速率约为1.14 cm/a, S01-1站约为0.69 cm/a, S01-2站约为0.35 cm/a, S02-3站约为0.623 cm/a (王永红, 2003; 庄克琳等, 2005; Liu et al, 2006; 王昕等, 2013; Qiao et al, 2017)。关于沉积速率对生物扰动作用的影响, 学者持有不同观点, Rhoads (1982)认为生物扰动强度与沉积速率呈现负相关关系, Moore等(1957)在密西西比三角洲研究时也指出当沉积速率超过4 cm/a时生物扰动几乎不发育, 覃雪波(2010)则认为高沉积速率可以带来更加丰沛的营养物质从而促进生物扰动作用进行。结合本文生物扰动强度数据以及前人实测的沉积速率来看, 沉积速率制约生物扰动作用, 当沉积速率过快时生物扰动构造发育微弱, 在长江水下三角洲及其邻近海域最适宜生物扰动构造发育的沉积速率约为0.3—0.6 cm/a。

4 结论

利用计算机断层扫描技术实现生物扰动构造的三维重构, 并以此为基础对长江水下三角洲生物扰动构造进行三维重构, 初步探讨了生物扰动构造的影响因素, 取得以下认识。

(1) 形成基于CT数字图像的重构沉积物岩心中三维生物扰动构造的方法和处理流程。

(2) 揭示了沉积物岩心CT值与含水率、沉积物粒级构成的关系。经过统计学分析, 发现CT值与沉积物含水率呈现显著的负相关关系; 与沉积物平均粒径的关系较为复杂, 与CT值存在二次项关系, 在浅表层与沉积物平均粒径呈现负相关, 在深部区域与沉积物平均粒径呈现正相关。

(3) 重构了长江水下三角洲4个站位的三维生物扰动构造, 揭示该区出现挖掘构造、觅食虫孔构造、逃逸虫孔构造、生物遗迹构造等显性扰动构造类型和垂向发育特征, 定量计算生物扰动强度。

(4) 该区生物扰动构造的类型和强度与底质类型密切相关, 黏土含量高的泥质沉积物不适宜生物扰动构造的发育, 而在砂质和粉砂质含量高的沉积物中易发育网状和根系状虫孔构造。沉积速率以及上覆水团性质也对生物扰动构造有着重要的影响, 沉积速率过快不利于生物扰动构造的发育, 正常温盐的海相环境中生物扰动构造较发育, 而在温盐梯度变化大的混合水团环境中生物扰动作用微弱。