东北黑土区侵蚀沟空间分布与形态的地质约束：以黑河地区为例

《CATENA》：Geological constraints on the spatial distribution and morphology of Erosion gullies in the black soil region of Northeast China: a case study of the Heihe area

　　

土壤侵蚀是指地表土壤及母质在水力、风力、重力、冻融及相关外力共同作用下发生的剥离、搬运和沉积过程（Shi等，2020；Chen等，2022；Sun等，2026）。这一现象被广泛视为环境问题（Lal，2004；Pimentel，2006；Van Oost等，2007；Zhang等，2023），会降低土壤肥力、加速土地退化，并威胁农业生产力和可持续土地利用（Peng和Wang，2012；Chi等，2021；Wei等，2022；Duan等，2026）。在全球尺度上，耕作区的侵蚀速率往往超过自然成土速率，导致土壤资源的持续损失和长期土地生产力的下降（Verheijen等，2009；Zhang等，2010；Lal，2015；Smith等，2016；IPCC，2019）。 东北黑土区是中国最重要的商品粮基地之一（Guo等，2025）。长期耕作、季节性降雨集中、冻融扰动以及坡地耕种导致该区域水土流失严重（Liu等，2020；Jiao等，2022）。沟蚀尤为具有破坏性，因为它将耕地切割成破碎斑块，移除肥沃表土，加速径流汇集，并减少可用农田面积（Zhang等，2016；Zhang和Li，2020；Liu等，2022；Li等，2024；Ma等，2024）。活跃的溯源侵蚀和耕型沟谷的岸坡崩塌还增加了生态恢复与土地整治的难度和成本（Duan等，2011；Qin等，2021；Wen等，2021；Wang等，2023）。 现有关于黑土区沟蚀的研究主要考察了地形因子（Jiao等，2025；Zhang等，2025；Zhang等，2026）、降雨侵蚀力（Li和Fang，2016；Xin等，2024；Dong等，2019）、植被覆盖（Hu等，2026）、土地利用及综合环境驱动因素（Zhu等，2025；Xu等，2026）。这些研究阐明了沟壑启动和扩展的重要外部驱动力。然而，沟壑是在地质框架内发育的：基岩风化成风化层，风化层供给成土母质，由此产生的近地表介质在质地、厚度、渗透性和黏聚力上存在差异，进而影响入渗、地下渗流以及沟道边界抵抗下切和坍塌的能力（Dardis和Beckedahl，1991；Li等，2002；Jin等，2023）。母质还塑造了土壤质地、结构、黏粒含量和抗蚀性（He等，2024；Peng等，2024）。在本研究中，研究人员使用术语“地质约束”（geological constraint）特指一个有限的意义：即地质基底对沟壑发育施加的物质和结构边界条件。该术语并不主张地质学在所有其他因素之外具有独立的因果效应；特别是，不同地质单元在高程和地势上也存在系统性差异，因此任何表观的地质关联都可能部分或完全由地形所中介。 这种混杂性是核心的方法论难题。简单地将沟壑清单叠加在地质图上无法区分：（i）在可比地形条件下仍存在的基底信号与（ii）地质仅仅通过差异性地貌演化继承而来的地形信号。另一个复杂因素是，岩性图与成土母质图并非独立的数据源：母质类别很大程度上源自下伏岩性，因此将两者作为独立证据处理存在重复计算的风险。 本研究以黑土区典型沟蚀影响的黑河地区为例，利用普查尺度的沟壑清单（5407条沟壑）和区域地质数据来探讨三个问题：（1）经面积归一化后，沟壑发生强度（密度）在不同岩性和母质单元间的差异有多大？（2）在考虑了地形、水文和土地利用协变量之后，地质变量是否能增加对沟壑发生强度和沟壑形态的解释力？（3）岩性与母质相对于彼此携带了多少独立信息，以及地形变异本身有多少与地质单元相关联？本研究的贡献不在于证明沟壑属性在不同地质单元间存在差异——叠置分析已反复证明了这一点——而在于量化了在地形框架之外，地质变量所增加的额外解释力，这通过面积归一化密度、效应量以及沟壑尺度和网格尺度的地形调整模型来实现。

研究背景与科学问题

土壤侵蚀是全球性的环境问题，会导致土壤肥力下降、土地退化和农业生产力受损。东北黑土区是中国重要的商品粮基地，但由于长期耕作、降雨集中、冻融作用和坡地开垦，该区域水土流失严重，其中沟蚀（gully erosion）尤为突出，它会切割耕地、移除表土、加剧径流并减少有效农田面积。已有研究多聚焦于地形、降雨、植被和土地利用等外部驱动因素，但沟蚀实际上发生在特定的地质框架内——基岩的风化产物决定了近地表介质的质地、厚度、渗透性和黏聚力，从而影响入渗、潜流和沟壁的抗侵蚀能力。然而，地质单元与地形之间存在天然耦合（如高地势区常对应抗蚀性强的基岩），使得简单叠加分析难以区分地质本身的独立作用还是地形中介效应。此外，岩性图（lithological map）与成土母质图（parent material map）高度同源，若同时使用可能造成信息重复计数。为此，研究人员以黑龙江省黑河地区的北安、五大连池和嫩江三市为研究区，基于5407条长度≥100 m的普查沟壑数据和1:25万地质图，提出了三个核心问题：（1）面积归一化后，沟蚀强度在不同岩性和母质单元间的差异有多大？（2）在地形、水文和土地利用协变量基础上，地质变量能否增加对沟蚀发生强度和形态的解释力？（3）岩性与母质的信息冗余程度如何，以及地形变异本身在多大程度上受控于地质单元？该研究发表于《CATENA》。

主要技术方法

研究人员使用了以下关键技术方法：①利用ArcGIS和Python（geopandas、rasterio、WhiteboxTools）进行空间数据处理与叠置分析；②采用12.5 m分辨率的ALOS PALSAR RTC数字高程模型（DEM）提取高程、坡度、局部起伏和汇流累积量，并划分子流域（subwatershed，共864个，阈值5 km2）；③使用第三次全国国土调查（2019年）和ESA WorldCover 2021 v200（10 m分辨率）确定土地利用和耕地掩膜；④基于Kruskal-Wallis秩和检验（rank-sum test）比较各形态指标（长度、纵向梯度、面积、宽度、弯曲度）在地质组间的差异，并用Dunn事后检验和Cliff's delta计算配对效应量；⑤建立广义可加混合模型（GAMM，通过R语言的mgcv包中的bam函数实现），在沟壑尺度上控制地形协变量（高程、缓冲坡度、局部起伏、汇流面积）、土地利用类别、二维空间平滑项和子流域随机效应后，分别加入岩性或母质作为分类固定效应，通过ΔAIC和偏差增量评估地质变量的附加解释力；⑥在1 km和2 km网格尺度上，采用Tweedie广义可加模型（GAM，log链接）拟合每格网内的总沟长，模型中纳入格网平均高程、坡度、起伏、耕地比例和二维空间平滑项，再添加岩性分类变量，同样通过ΔAIC和偏差增量衡量地质信号的贡献。

研究结果

4.1 沟蚀发生强度在地质单元间的差异

经面积归一化后，沟蚀密度在不同地质单元间差异显著。在主要岩性单元中，长度密度（D_L）从更新世玄武岩的0.174 km·km?2（富集比R=2.31）到晚三叠世—早侏罗世花岗岩的不足0.001 km·km?2（R=0.011），跨度超过两个数量级。当分母限制为耕地面积后，更新世玄武岩的耕地限制富集比（R_crop）仍高达2.52，而晚三叠世—早侏罗世花岗岩即使在其少量耕地上也几乎无沟（R_crop=0.065），表明低耕地暴露不能单独解释花岗岩单元的缺沟现象。

4.2 沟壑形态在岩性和母质单元间的差异

未调整的Kruskal-Wallis检验显示所有五个形态变量在地质组间均有显著差异，但效应量差别很大：纵向梯度的效应最大（ε2=0.101），其次是宽度（ε2=0.057）和面积（ε2=0.024），长度（ε2=0.006）和弯曲度（ε2=0.003）的效应极小。具体而言，早白垩世基性火山岩上的沟壑中位梯度最高（8.77°），嫩江组最低（3.08°）；宽度中位数从林西组的25.1 m到孙吴组的17.9 m不等。母质分组得到相似排序且效应略小。

4.3 地形调整后的关联

4.3.1 沟壑形态：地形基线模型（M0）解释了纵向梯度77.3%、宽度56.8%、面积42.7%和长度31.1%的偏差。加入岩性（M1）或母质（M2）后，增加的解释偏差最多仅为0.24%（长度），宽度甚至无改善（ΔAIC≈0）。由于200 m缓冲坡度与纵向梯度响应本身高度相关（r=0.78），研究人员认为宽度和面积的调整后结果更能反映形态主要由地形中介的结论。

4.3.2 沟蚀发生强度：在1 km网格尺度上，加入岩性使地形+耕地模型的解释偏差从39.63%升至40.90%（增加1.27个百分点，ΔAIC=191.7）；2 km尺度增加0.98个百分点（ΔAIC=54.3）。调整后的岩性系数显示：孙吴组（+0.605）、嫩江组（+0.418）和更新世玄武岩（+0.462）在同等地形和耕地暴露下为正关联，而晚石炭世—早二叠世酸性火山岩为负关联（?0.715）。增大空间平滑基维数后ΔAIC降至111.3，表明部分地质信号存在空间混淆。

4.4 岩性—母质冗余与地形的地质继承性

岩性与母质分类的Cramér's V=0.993，两者近乎完全冗余。地质单元隶属关系自身可解释网格高程方差的46.3%、局部起伏的34.0%和坡度的12.9%，表明地形协变量本身已继承了大量的地质信息。

4.5 敏感性分析

改变跨边界沟壑分配的主导阈值（60%/70%/80%）和最小分组样本量（30/50/100条沟），主要统计结果保持稳定。

讨论与结论

讨论部分指出，本研究的核心实证结果是沟蚀发生与形态在地质关联上的不对称性：发生强度在调整地形和耕地暴露后仍存在微弱但可检测的地质信号（1 km网格ΔAIC=192），而形态差异则基本可由地形中介解释。孙吴组在原始密度上低于区域均值（R=0.61），但在调整模型中表现为正关联，这是因为其占据平缓地形——地形模型本就预测低密度——一旦控制这种有利地形，其弱固结砂岩-砾岩反而成为高于期望的产沟单元。花岗岩单元的近乎无沟状态在耕地限制密度下依然成立（R_crop=0.065），加强了基底抗蚀性解释，但耕地选址的内生性问题（陡峭薄层花岗岩本身不适宜耕作）意味着这不是完全的因果分离检验。岩性与母质的近完全冗余（V=0.993）表明二者构成同一地质信号的两个分类层级而非独立证据。研究人员提出了一个田野约束的概念模型，将各环节按证据等级分为直接观测（A）、地形调整统计关联（B）、剖面约束的过程解释（C）和待验证假说（D），其中花岗岩缺沟、密度对比属于A级，孙吴组和嫩江组的调整后系数属于B级，孙吴组渗透层覆铁质胶结层的剖面结构与潜蚀机制的关联属于C级，而冻融与不同基质的交互作用等属于D级。对于黑土保护，建议将地质图纳入侵蚀风险区划以补充地形信息，但工程设计参数仍应以地形属性为主。

研究结论翻译如下：

(1) 结合5407条普查沟壑、1:25万岩性和母质图、12.5 m DEM及2019年土地利用数据，本研究通过面积归一化密度、效应量和地形调整模型，量化了黑河地区沟蚀发生与形态的地质关联。

(2) 经面积归一化后，沟蚀发生强度在地质单元间差异显著，从更新世玄武岩的0.174 km·km?2（R=2.31）到晚三叠世—早侏罗世花岗岩的低于0.001 km·km?2。孙吴组在原始计数中的表观优势是露头面积假象（R=0.61）。

(3) 网格模型检测到超出所包含地形和耕地协变量的地质信号：岩性在1 km（ΔAIC=192）和2 km（ΔAIC=54）尺度上改善了Tweedie模型，酸性火山岩呈负调整关联，更新世玄武岩及弱固结的孙吴组和嫩江组呈正调整关联；花岗岩单元即使在自身耕地上也几乎无沟（耕地限制富集比低至0.065）。

(4) 沟壑形态在未调整检验中单元间差异显著——纵向梯度（ε2=0.10）和宽度（ε2=0.06）最强，长度仅可忽略（ε2=0.006）——但一旦控制地形、土地利用和子流域隶属关系，岩性或母质最多仅增加0.24%的偏差。鉴于地质单元可解释46%的高程方差和34%的起伏方差，观察到的形态格局与地质继承性地形的实质性中介作用一致。

(5) 岩性与母质分类近乎完全冗余（Cramér's V=0.993），构成同一地质信号的两个分类层级，而非两条独立证据链。

(6) 对于黑土保护，底质图为基于地形的沟蚀风险区划增添了识别沟壑集中区域的信息，而地形属性仍是预测沟壑形态的主要依据；这些推论需通过监测的沟壑发育进行验证。