多尺度形状可控沙粒生成方法：傅里叶合成与遗传算法优化相结合

《Computers and Geotechnics》：Multi-scale shape-controlled sand particle generation method: Fourier synthesis coupled with genetic algorithm optimization

　　摘要：沙粒的多尺度形态特征是决定其宏观力学行为的关键因素。本文提出了一种基于程序几何的颗粒生成模型，可独立控制多尺度形状参数，即伸长指数EI、角形指数AI和粗糙度指数RO。该方法首先通过分形参数控制的傅里叶谱合成结合针对性后处理，生成具有自然多尺度纹理的初始颗粒表面，再利用基于高斯核变异的遗传算法迭代优化颗粒形态。该算法通过采用二十面体网格避免球坐标系带来的畸变，并通过加权适应度函数实现多参数优化，最终生成形状特性可精确控制的颗粒模型。将所生成颗粒与真实沙子的X射线计算机断层扫描数据对比后发现，其形态参数的统计分布与真实沙子高度一致。最后，通过离散元法模拟验证了所生成颗粒在研究形态-力学关系方面的有效性，同时还利用3D打印技术成功制备出了适用于透明土壤物理模型试验的样品。本研究为构建具有可控多尺度形态的沙粒模型提供了一种可靠有效的方法。

1. 引言
天然颗粒材料在岩土工程中得到广泛应用。作为控制颗粒互锁与接触力学的关键微观参数，颗粒形状已被证实会显著影响材料的刚度、强度和抗变形能力等宏观力学性能（Cho等人，2006；Rousé等人，2008；Yang和Luo，2015；Nie等人，2022、2025；Gong等人，2024；Ma等人，2024；Shi等人，2024）。

为系统描述颗粒形状，Barrett（1980）提出了广泛采用的多尺度框架，将形状分为三个独立尺度：大尺度上的整体形态、中尺度上的圆度与角形度，以及小尺度上的表面纹理（Zheng和Hryciw，2015；Su和Yan，2018；Guida等人，2020；Liu等人，2023；Khan和Latha，2024）。大量研究表明，颗粒形状特征对宏观力学性能的影响机制在不同尺度上存在显著差异（Li，2013；Suh等人，2017；Lashkari等人，2020；Rorato等人，2020；Wai等人，2022；Mondal等人，2025；Mousavi等人，2025）。然而，由于技术手段和实验条件的限制，大多数传统实验研究仅聚焦于一两个尺度。这主要是因为天然颗粒的不规则性使得在实验室条件下难以直接量化和控制不同尺度上的颗粒形状（Altun等人，2011；Sarkar等人，2019）。幸运的是，玻璃珠和3D打印透明土壤等材料为解决这一问题提供了有效方案。高度标准化的玻璃珠能够有效消除因天然沙子形状不规则而带来的实验不确定性（Yang和Luo，2015；Alshibli和Cil，2018；Xiao等人，2019a、b、c；Miao等人，2025）。更为先进的是，3D打印透明土壤不仅能够根据数字模型精确制造出几何结构极为复杂的颗粒，还能直接无损地观察其内部三维结构及相互作用（Gupta等人，2019；Ahmed和Martinez，2021；Liang等人，2022；Anusree等人，2024；Liang等人，2024；Tong等人，2024a、b；Chen和Tong，2025；Liang等人，2025）。

与此同时，离散元法模拟也已成为研究颗粒形状效应的重要工具，因为它能够精细控制颗粒尺度特征并直接观测微观力学信息（Cundall和Strack，1979；Yimsiri和Soga，2010；Gong和Liu，2017；Feng等人，2020；Nie等人，2020a、b；Huang等人，2023；Wang等人，2025）。Xiong等人（2022）通过结合X射线CT扫描与离散元法模拟，系统揭示了不同长度尺度下颗粒形态对颗粒土微宏观力学行为的影响，并进一步建立了考虑多尺度颗粒形态特征的福建河沙CT基离散元建模框架（Xiong等人，2025）。不过，这类基于CT的离散元研究高度依赖于真实颗粒几何模型的获取，而这也是3D打印物理实验和离散元数值模拟的共同前提条件。早期研究人员开发了激光扫描（Latham等人，2008；Guo等人，2018）、计算机断层扫描（Katagiri等人，2010；Zhou等人，2015；Kootahi等人，2024）以及图像处理（Guida等人，2020；Ruidong等人，2025）等技术，用于提取和表征真实颗粒材料的形状。尽管这些方法在技术上较为有效，但仍然受到高昂成本和时间限制。

为获得大量形状可控制的颗粒模型，研究人员开发了多种随机颗粒生成方法。其中，基于球面谐波（SH）（Grigoriu等人，2006；Liu等人，2011；Zhou和Wang，2017；Su和Yan，2018；Wei等人，2018a、b；Zhou等人，2018；Sun和Zheng，2021）和傅里叶级数（Thomas等人，1995；Bowman等人，2001；Mollon和Zhao，2013、2014）等数学描述符的方法被广泛应用。这类方法的核心思想是通过与所选描述符相关的系数组来描述颗粒形状。其处理流程通常是从重建CT扫描得到的颗粒形状开始，然后通过特征分析算法调整这些系数以生成新系数，从而得到具有可控且相似形态特征的新颗粒。Wei等人（2018a）提出了一种基于旋转不变球面谐波分析的颗粒克隆技术，利用两个旋转不变球面谐波因子来控制颗粒形状。Mollon和Zhao（2013）基于对真实沙粒二维傅里叶谱的统计分析，建立了3D颗粒生成的参数化框架，能够根据轮廓参数合成3D颗粒。Xiong和Wang（2021）进一步将主成分分析引入球面谐波表示体系，并提出了颗粒形状的形态基因突变方法，为多尺度形态特征的调控提供了分析范式。然而，这些方法通常仅限于“星形”颗粒（即从表面到中心的任意直线都完全位于颗粒内部），限制了其几何通用性。Vlahini?等人（2017）尝试将水平集函数作为描述符，虽然这种方法获得了高保真度，但接触检测的计算成本远远高于基于简单几何结构的传统离散元方法。为了解决效率问题，人们提出了基于元球（Zhang等人，2021、2023；Zhao等人，2023）和多超椭球（Zhao等人，2020；Zhao和Zhao，2019）的基于描述符的重建方法，但这些方法仍难以准确重建复杂形状的颗粒。此外，计算这些方法的描述符及其系数需要大量资源，而且将这些系数与直观的形状参数（如角形度、粗糙度）关联起来也相当困难。

近年来，基于深度学习的变分自编码器为颗粒生成提供了新途径（Bourilkov，2019）。变分自编码器建立在概率框架之上，能够自动学习并表征复杂形态的潜在分布（Kingma和Welling，2013）。其编码器将输入颗粒映射到低维潜在空间，解码器则从该空间重建并生成新颗粒（Higgins等人，2017）。这种数据驱动的方法能够有效捕捉并再现真实颗粒形态的统计特征，生成高度逼真的样本（Tan等人，2018；Shi等人，2021）。不过，该方法依赖于大规模、高质量的数据集，而这类数据集的构建需要大量资源，且生成颗粒的多样性本质上受训练数据分布的限制。因此，亟需开发一种便捷高效的方法，用于生成具有可控多尺度形状参数且与真实沙粒高度相似的3D颗粒模型。

本文提出了一种两阶段方法，用于生成具有可控多尺度形状特征的颗粒形态。该方法首先通过分形参数控制的傅里叶谱合成结合针对性后处理，生成具有自然多尺度纹理的初始颗粒表面，然后采用集成高斯核变异的遗传算法迭代优化颗粒形态。该算法通过使用二十面体网格避免球坐标系畸变，并通过加权适应度函数实现稳定的多参数优化，从而生成几何特性可精确控制的颗粒模型。通过比较所生成颗粒与真实CT扫描沙粒的形状参数统计一致性，验证了所提方法的有效性。最后，通过在离散元模拟和3D打印物理模型实验中应用该方法，尤其是用于透明土壤建模，进一步证明了其实际应用价值。

2. 颗粒形状量化方法
2.1 颗粒形态获取
为建立可靠的基准数据集，以便后续进行形态特征表征和方法验证，本研究采用了Zhou等人（2018）提出的成熟颗粒成像与重建方法。测试用沙料、μCT扫描系统、样本尺寸以及3D图像重建参数均与原始文献设置保持一致，以确保数据的可比性和实验的可重复性。本研究选择了两种具有不同形态特征的天然沙子作为研究对象：来自英国的莱顿巴扎德沙，由于长期地质搬运作用，其形态圆润光滑；来自香港的高分解花岗岩沙，其特征为角形突出、表面粗糙度较高且表面凹坑发育明显。共对这两种沙子的120个颗粒样本进行了扫描，所用设备为Carl Zeiss Metrotom 1500 X射线μCT系统，3D图像的体素分辨率约为6.5微米。

后续的图像处理遵循文献中记载的标准流程，包括相位分离、降噪以及分离相互接触的颗粒。随后利用球面谐波函数重建每个颗粒的表面，从而得到连续且可微分的表面模型。这一重建结果为后续计算用于量化颗粒形状、角形度和表面粗糙度的3D形态描述符提供了基础。图1展示了CT扫描得到的颗粒的体积可视化图像以及通过球面谐波函数重建的模型对比图。

下载：下载高分辨率图像（291KB）
下载：下载全尺寸图像
图1. μCT扫描得到的沙粒颗粒的体积可视化图像与基于球面谐波函数重建的颗粒模型对比图

2.2 球面谐波重建
通过X射线CT扫描获得的颗粒3D表面信息表现为笛卡尔坐标系下的点云P(x, y, z)。这一信息也可用球坐标系表示，即以颗粒质心O为原点的极径r(θ, φ)，其中θ∈[0, π]为极角，φ∈[0, 2π]为方位角，r表示径向距离。在球面谐波分析中，每个顶点的径向函数r(θ, φ)可展开为一系列球面谐波基函数的组合：

(1)r(θ,φ)=∑n=0N∑m=-nnanmYnm(θi,φi)

其中N为球面谐波基函数的最大阶数；n和m分别表示阶数和阶次；anm为复数形式的球面谐波系数；Ynm为归一化的球面谐波函数，其表达式为公式(2)：

(2)Ynm(θ,φ)=(2n+1)(n-|m|)!4π(n+|m|)!Pn|m|(cosθ)sin(|m|φ)

其中Pnm为勒让德多项式，其表达式为公式(3)：

(3)Pnm(x)=(-1)m2nn!(1-x2)m/2dn+mdxn+m(x2-1)n

其中n和m分别为勒让德多项式的阶数和阶次。每个沙粒的球面谐波系数是通过最小二乘法确定的。利用这些系数，即可根据公式(1)计算出任意极角θ和方位角φ对应的极径r，进而重建颗粒表面。

在这一过程中，最大阶数N的选择至关重要，因为它决定了形态细节与计算效率之间的平衡：当N较小时，只能捕捉到大致轮廓，而增大N则有助于逐步恢复更细微的特征。不过，超过一定阈值后，继续增大N对提升精度的作用有限，但会导致计算成本急剧上升。许多研究（Wang等人，2021；Fan等人，2024a、2024b）表明，当球面谐波基函数的最大阶数≥15时，已足以表征颗粒表面的局部形态特征。因此，本研究中选择N=15。

2.3 颗粒形态量化
对颗粒形态进行系统量化是岩土材料多尺度表征的基础。依据Barrett提出的三尺度框架（图2），沙粒形态在三个不同尺度上进行表征：整体形态（宏观尺度）、角形度（中观尺度）以及表面纹理（微观尺度）。图3详细展示了这种多尺度定量分析的示意图。在本研究中，选择伸长指数（EI）、角度指数（AI）和粗糙度指数（RO）作为粒子生成的核心输入参数。这三个参数直接对应于Barrett经典框架中的三个层次结构，能够全面覆盖从整体轮廓到细微表面纹理的粒子形态特征。每个参数分别反映不同空间频率下的形态信息，且彼此之间具有较好的独立性。这一三参数体系已被广泛应用于岩土工程中的粒子形态研究（Su和Yan，2018；Fan等人，2024a）。本研究采用的三个形状描述符定义如下。下载：高分辨率图像（82KB）下载：完整尺寸图像图2. 沙粒的多尺度形态。下载：高分辨率图像（625KB）下载：完整尺寸图像图3. 多尺度形态的量化 (a) AI，(b) RO。粒子形态通常通过伸长指数来量化，该指数是根据粒子最小边界长方体的最长轴（Lc）和中间轴（Lb）的比值计算得出的。伸长指数（EI）由公式(4)定义：(4) EI = Lb/Lc。粒子角度描述了其尖锐程度或圆润程度，属于中尺度特征。在各种角度指数中，Fan等人（2024a）提出的基于半径的角度指数（AI）因具有较高的计算效率和理论严谨性而被选用：(5) AI = 1/k∑i=1^k r(θi,φi) / [rEE(θi,φi)]，其中k为节点总数；θi和φi分别为节点处的方位角和极角；r(θi, φi)表示节点i处重构粒子的极径，rEE(θi, φi)则表示相同球面坐标下等效椭球的极径。该等效椭球是通过使其主轴与粒子主轴重合得到的。图3(a)展示了AI的定义。为避免表面粗糙度对角度计算的影响，粒子表面是使用最大阶数为N=8的SH函数进行重构的。定义为微观尺度纹理的粒子粗糙度，指的是粒子表面的细微不规则性。与上述角度指数的计算方式类似，粗糙度指数也是通过在不同最大阶数下进行SH函数重构来计算的：平滑基准面对应的N=8，高分辨率重构面对应的N=15。图3(b)展示了RO的示意图。粗糙度是通过平滑基准面与粒子重构面之间的偏差来量化的，该偏差可通过公式(6)计算：(6) RO = 4π/3 * V^(3/1) * S / ∑i=1^k d_i * S_i，其中V为基准面的体积；S表示基准面的面积；d_i为平滑基准面与实际表面之间的局部偏差；S_i表示第i个三角形网格的面积。3. 方法论本章提出了一种用于生成具有多尺度形态特征的3D粒子模型的综合方法框架，如图4所示。该框架以目标EI、AI和RO作为输入。在傅里叶谱合成阶段，通过分形控制的傅里叶谱生成四个表面分量矩阵，然后对这些矩阵进行线性叠加并经过后处理，从而得到具有自然多尺度纹理的初始粒子。接着，通过双线性插值将初始粒子的球坐标网格映射到四阶精化的二十面体网格上，以消除极向畸变。进入遗传算法阶段后，会生成初始粒子群体。对于每一个个体，都会根据其形态优化需求进行相应处理：对过于平滑的粒子施加椭圆高斯粗糙化处理，而对过于粗糙的粒子则施加高斯平滑处理。经过变异操作后，会计算形状参数的适应度，再通过精英选择保留表现最佳的粒子。该算法不断重复这种变异-选择循环，直到适应度收敛到预定的阈值，最终输出标准STL格式的3D粒子模型。下载：高分辨率图像（557KB）下载：完整尺寸图像图4. 整体方法流程图。3.1 初始表面高度矩阵的生成本节介绍了通过合成多尺度地形特征来构建初始粒子表面高度矩阵的方法。该过程主要包括两个阶段：(a)通过基于傅里叶的合成方法，利用受控的分形参数生成不同尺度的多个表面分量矩阵；(b)对这些分量进行线性叠加和处理，从而得到初始高度矩阵Zinit。(i) 通过分形傅里叶合成得到单组分表面在球坐标系下，通过基于傅里叶的谱面合成方法，可以创建出初始的表面高度矩阵Z(θ, φ)，其中θ∈[0, π]表示极角（纬度），φ∈[0, 2π]表示方位角（经度）。矩阵中的每个元素Zi,j(θi, φj)表示角度网格点(θi, φj)处的径向距离，其中i和j分别对应极向和方位向的网格索引。需要明确的是，在本文中，i和j被用作默认的计数索引，始终代表连续的正整数。共创建了四个独立的组分矩阵：Zi（i=1, 2, 3, 4），每个矩阵都由特定的分形参数组合决定，这些参数包括表面粗糙度方差（σ2）、控制分形维数的赫斯特指数（H）、截止波矢量（qr）、x方向的长度（Lx）以及角度分辨率（mmesh, nmesh）。我们首先定义空间离散化方式，这是谱面合成的基础。计算域被划分为规则的网格，x方向上的空间分辨率（像素大小）为Δ=Lx/mmesh，y方向上的总长度为Ly=Δ×nmesh。为了能够在频域中对高度矩阵H?进行谱分析，我们通过公式(7)定义波矢量分量：(7a) qx[i] = 2πΔ/mmesh*i - Δ2/mmesh2，(7b) qy[j] = 2πΔ/nmesh*j - Δ2/nmesh2，其中i=0, 1…, mmesh?1，j=0, 1…, nmesh?1。随后，通过离散傅里叶变换可以得到二维谱表示H?(qx,qy)：(8) H?(qx,qy) = ∑i=0^mmesh?1∑j=0^nmesh?1 Z[j,i]exp(?i(qx*xi + qy*yy))。对于此处所有的二维矩阵，其元素用[行, 列]表示，行索引对应y轴/极角（θ），列索引对应x轴/方位角（φ）。因此，Z[j, i]表示网格点(i, j)处的极径。为了合成具有分形特性的表面，需要根据赫斯特指数H和截止波矢量qr生成功率谱密度矩阵Cq，其表达式为公式(9)：(9a) Cq(ρ) = {qr?2(H+1)}, ρ P90则取|Zinit[j,i]|·|P90|；其他情况下则保持Zinit[j,i]不变，其中P35和P90分别代表Zinit的第35百分位数和第90百分位数。其次，为消除高频噪声和较小的随机波动，同时保留主要的形态特征，采用了二维中值滤波器。该滤波器通过用局部核内的强度中值替换每个像素值，从而有效去除孤立的尖峰，而不模糊清晰的边缘。核的大小会根据所选的网格分辨率进行调整，使其在球面上对应的角度范围约为15 - 20°。这样的尺度既能去除数值上产生的、在地质意义上无关紧要的伪影，又能保留有意义的形态特征。第三，为确保初始粒子群体具有与目标形态极为相似的几何特性，同时减少后续优化所需的迭代次数，需要对球坐标下的粒子表面高度函数r(θ, φ)进行调制，使其符合目标EI的要求，调制后的高度r′(θi, φj)由公式(18)给出：(18) r′(θi,φj) = EI·[cos2(θi) + (EI·sin(θi))2]/[r(θi,φj)rmax]，其中θi和φj分别是离散的极角和方位角；rmax表示复合高度矩阵的最大值。最后，为实现无缝的球面映射，通过余弦加权混合的方式来保证边界连续性。在两极区域，数据会平滑过渡到方位平均高度；而在经向边界上，则通过与对侧区域进行混合来确保周期性。这一处理是在一个规定的过渡区域内进行的，如图7所示，其公式如下：(19a) w(si) = 1 ? cos(πsi)2，(19b) Zf = 1 ? w(si)·Zori + w(si)·Ztar，其中si是过渡区域内的归一化位置；w(si)是边界区域中位置si处的权重系数；Zf为最终的高度矩阵；Zori是要连接并平滑的区域的高度；Ztar是对侧区域用于连接和平滑的目标高度。下载：高分辨率图像（490KB）下载：完整尺寸图像图7.南向与极向连续性处理。该处理的结果为具有几何适配性及平滑边界的高度矩阵Zf(θ, φ)。图8(a)展示了以极角θ和方位角φ为函数的表面径向距离r的2D分布情况，而图8(b)则呈现了对应的3D球坐标模型，该模型可作为种群中高质量初始个体，用于后续的遗传算法优化。下载：下载高分辨率图像（177KB）下载：下载全尺寸图像图8。(a)所生成的初始表面高度矩阵的2D示意图；(b)对应初始粒子的3D示意图。3.2基于遗传算法的高斯核变异形状优化：第3.1节中生成的初始粒子表面具有天然的多尺度纹理特征，但其形态参数（EI、AI、RO）呈随机分布，无法与真实沙粒的目标值相匹配。为此，采用基于高斯核变异的遗传算法优化方法对形状进行迭代细化。同时，由于初始表面的球坐标网格存在极点奇异性畸变，会导致形态参数出现系统性偏差，因此引入均匀的二十面体网格来消除参数畸变。以下为遗传算法优化的步骤。(i)二十面体网格映射首先，将球坐标网格转换为均匀分布的二十面体网格，如图9所示。球坐标中的极点奇异性会导致顶点密度出现显著变化，进而引发面积畸变和计算误差。在接近极点的区域，相同的角度增量（Δθ、Δφ）对应的表面积会大幅减小，这会导致变异操作产生畸变的表面，并使形状、尖锐度和粗糙度等描述指标出现系统性偏差。二十面体网格通过均匀的顶点分布消除了此类参数畸变，从而实现空间上一致的变异和无偏的几何评估。这种方法为基于遗传算法的粒子重建建立了稳定且精确的框架。下载：下载高分辨率图像（410KB）下载：下载全尺寸图像图9.相同粒径下的球坐标网格与二十面体网格。重建过程从一个由20个等边三角形面和12个顶点构成的单位二十面体开始。随后，该网格通过L级递归细分进行均匀细化。在每一级中，每个三角形都会通过在边中点插入新顶点的方式被分割成四个，然后再投影到单位球面上。这样就能得到一个高度均匀的网格，其顶点数Nv和三角面数F分别满足公式(20a)Nv=10×22L+2以及(20b)F=20×4L。由于更高的L值能生成更密集的网格，从而在更高的计算成本下捕捉更多细节，因此选择L值时需要权衡利弊。为了在形态细节与计算效率之间取得平衡，本研究始终采用L=4的均匀二十面体网格。这种分辨率已被证明能够为基于SH的粒子重建提供最佳平衡，既能有效捕捉沙粒、道碴粒和卵石粒的几何特征，又能将形状描述指标的平均误差控制在0.5%以下，粗糙度指标的误差控制在3%以内（Fan等人，2025年）。为了在该网格上表示粒子，需利用双线性插值将规则球坐标网格上的高度数据映射到网格的顶点上。具体而言，首先将球面高度矩阵Zf(θ, φ)沿经度方向周期性扩展，覆盖三个完整周期（方位角从-540°到540°），以避免边缘伪影。然后，在这个扩展后的区域内对每个网格顶点(θi, φi)的径向距离ri进行插值。所有仍以球坐标系描述的顶点Zi(ri θi, φi)都会以Z′i(xi, yi, zi)的形式系统地存储在生成的二十面体高度矩阵中。(ii)粒子形状变异操作优化过程首先生成由Np个候选粒子组成的初始种群，每个粒子用Zi(θ, φ)表示。为了逐步将形状演化为目标形态，会对均匀的二十面体网格应用基于高斯核的变异操作。高斯核提供了一种模拟自然表面演变的数学工具，能够生成平滑且随距离呈衰减的扰动。这样一来，修改效果会随着与变异中心的距离增加而减弱，从而保持几何一致性，避免出现人为的不连续性，如图10所示。具体的操作细节如下。下载：下载高分辨率图像（216KB）下载：下载全尺寸图像图10.球面高度场的高斯变异：(a)变异前；(b)变异后。首先，随机选择一个顶点作为高斯核变异的中心。对于这个坐标为rc(θc, φc)的中心顶点vc，变异的空间范围和幅度由参数ks和kc控制。这些参数可根据优化目标进行设定。例如，若要细化微观粗糙度，则可降低kc的值。受影响区域是通过在网格邻接关系上执行广度优先搜索（BFS）来确定的。所有满足3D欧几里得距离dg < 2σg的顶点vv rv (θv, φv)都包含在变异区域内：(21a)x=rsinθcosφ y=rsinθsinφ z=rcosθ(21b)dcv=[(xc-xv)2+(yc-yv)2+(zc-zv)2](21c)σg=U1ksZˉ其中，σg表示高斯核的半径，U1是介于0.3到1之间的随机数，Zˉ则是矩阵Z的均值。为了使不同粒子的优化能够更快收敛，根据粒子表面粗糙度的特性，设计了两种变异方案：(1)椭圆高斯粗化操作当需要增强表面粗糙度时，可采用各向异性的高斯扰动。图11(a)展示了粗化前后粒子的对比情况。一般的变异公式为：(22)Zv′=Zv+ΔZ(vv,vc)其中，ΔZ (vv, vc)表示相对于中心vc，顶点v处的变异扰动。标准的（各向同性的）高斯核定义为：(23)f(x)=e(-x2/(2σg2))其中x是偏离中心的距离。这种扰动在各个方向上都是均匀施加的。这种各向同性特性适用于模拟在均匀条件下形成的特征。然而，天然沙粒的形态往往受到各向异性的物理过程影响，比如优先风化、沉积物输送方向以及分层的矿物组成，这些因素会导致沙粒表面出现拉长或定向的特征。为了更准确地再现这些自然的、不规则的纹理模式，我们采用了椭圆高斯核函数在球面上实现局部变形。对于操作范围内的所有顶点Vv，椭圆高斯核值ΔZ(vv, vc)的计算方式如下。下载：下载高分辨率图像（842KB）下载：下载全尺寸图像图11.变异操作前后的粒子：(a)高斯粗化操作，(b)高斯平滑操作。首先，如图12所示，在Vc的切平面上构建正交基向量u和v，其中法向量nc则沿径向指向：(24)nc=[xc,yc,zc]/rc接着，计算Vv相对于点Vc的位移向量在切平面上的投影：(25a)dcv=xc-xv yc-yv zc-zv(25b)pu=dcv·u(25c)pv=dcv·v下载：下载高分辨率图像（472KB）下载：下载全尺寸图像图12.椭圆高斯核。随后，为模拟各向异性的变异模式，我们引入了两个随机参数：控制椭圆偏心率的纵横比αg ∈ [0.4, 1.0]，以及定义椭圆方向的旋转角θg ∈ [-180°, 180°]。坐标(u, v)会通过θg进行旋转，使其与椭圆的主轴对齐：(26)u′v′=cosθgsinθg - sinθgcosθguv接着，计算标准化的椭圆距离Dp：(27)Dp2=u′σg2+v′αgσg2最后，高斯核值可通过以下公式得到：(28)ΔZ(vv,vc)=ZˉMg*exp[-Dp2/2]其中，Mg在[-kc/2, kc/2]范围内，是随机幅度系数。(2)高斯平滑操作作为进一步粗化的补充手段，高斯平滑采用标准核函数，旨在在保留重要形态特征的同时，减轻表面的过度尖锐度和高频粗糙度。平滑的效果如图11(b)所示。具体操作步骤如下。对于操作范围内的所有顶点Vv，其高斯核权重为：(29)wcv=exp(-dcv2/2σg2)对于权重wv > 0.01的任意顶点，需确定其两层扩展域N(v)。对该域内的每个点vj，按照公式(30)计算局部高斯权重wvj和综合权重vvj：(30a)wvj=exp(-dcv2/2(σg2)2)(30b)vvj=wvj·[1+0.3wcj]综合权重不仅能确保平滑中心的全局影响，还能保证相邻点的连续性。顶点v的加权平均高度为：(31)Zˉv=∑j∈N(v)vvj·Zj/∑j∈N(v)vvj此外，顶点v的平滑后的高度为：(32a)βv=4Zˉkc·wcv(32b)Zv′=(1-βv)·Zv+βv·Zvˉ(iii)适应度评估与迭代遗传算法优化每次变异操作之后，都会使用公式(4)-(6)计算不规则粒子的3D形状参数EI、AI和RO。在遗传算法框架中，单个粒子fi的适应度被定义为计算得到的参数与目标值之间的加权绝对偏差。由于这些参数（EI、AI和RO）的数值量级通常差异较大（例如，对于大多数粒子，EI ∈ [0.2, 1.0]，AI ∈ [0, 0.2]，RO ∈ [0, 0.04]）（Fan等人，2024a），如果直接对它们的绝对偏差求和，数值较大的参数就会主导适应度评估，从而使搜索过程出现偏差。为让每个形状参数对总适应度的贡献相等，消除尺度偏差，我们采用了基于范围归一化的权重策略，其核心理念是：每个参数可能出现的最大绝对偏差（即其完整的数值范围）应对应相等的总适应度贡献；具体而言，我们首先根据天然沙粒的统计特性，计算出每个形态描述指标的完整数值范围，得到伸长率指标的范围ΔEI = 1.0 - 0.2 = 0.8，尖锐度指标的范围ΔAI = 0.2 - 0 = 0.2，粗糙度指标的范围ΔRO = 0.04 - 0 = 0.04，然后根据参数数值范围的倒数确定权重比例，即ω1:ω2:ω3 = 1/ΔEI:1/ΔAI:1/ΔRO = 1.25:5:25，最后通过将所有数值除以1.25，将最小的权重标准化为1，从而得到最终的默认权重系数ω1 = 1，ω2 = 4，ω3 = 20。基于上述归一化权重方法，可构建出加权适应度函数：(33)fi=ω1|EI-EIc|+ω2|AI-AIc|+ω3|RO-ROc|其中，EIc、AIc和ROc分别是伸长率指标、尖锐度指标和粗糙度指标的目标值；ω1、ω2和ω3则是相应的加权系数，可根据特定研究场景中的目标形状参数进行灵活调整。适应度fi越低，说明粒子的形态与目标形态的契合度越高。基于这一指标，会采用精英选择策略保留表现最佳的粒子，这些粒子会被带入下一代，继续进行遗传操作。通过不断的变异和选择，种群逐渐向具有目标形状参数的粒子靠拢。如图13所示，迭代会持续进行，直到最佳粒子的适应度降至预定的阈值以下（例如0.02），此时该粒子就会被导出为STL格式的最终目标3D模型。图14证明了所提出的算法能够生成具有多种形态特征的粒子，并且在可定制的形状控制方面表现出极高的精度。需要指出的是，从算法层面来看，这三个形状参数是相互独立的。在天然沙粒的合理统计范围内（EI ∈ [0.2, 1.0]，AI ∈ [0, 0.2]，RO ∈ [0, 0.04]），通常可以稳定地生成合理的星形粒子。但是，若参数组合随意，就无法得到有效的粒子结果。那些参数间差异过大的极端跨尺度组合，会给几何协调带来一定困难，也无法保证稳定的收敛。例如，在近乎光滑的球形基底上叠加过多的高频纹理，很容易打破星形的几何约束，导致表面出现局部自相交现象，同时遗传算法的优化效率也会显著下降。因此，在优化过程开始之前，可以进行初步的可行性检测，确保输入的目标参数处于合理且可实现的范围内。下载：下载高分辨率图像（190KB）下载：下载全尺寸图像图13.粒子通过遗传算法经过迭代处理。下载：下载高分辨率图像（268KB）下载：下载全尺寸图像图14.具有不同形状参数的生成粒子的形态：(a)具有不同EI值的粒子，(b)具有不同AI值的等轴粒子，(c)具有不同RO值的等轴粒子。4.所提方法的验证为了验证所提出的基于遗传算法的随机粒子生成方法的有效性和准确性，我们使用了通过X射线微计算机断层扫描（X-μCT）数据重建的真实沙粒进行了对比研究。为了重建真实沙粒的形态，首先使用SH展开法处理从X-μCT扫描中获取的3D表面数据（见2.1节）。随后，为所有120个重建得到的粒子计算了2.3节中定义的三个形状描述指标（El、AI、RO）。这些通过实验得到的形状参数分布被用作目标输入，利用本研究开发的程序化方法生成了一组对应的虚拟粒子。最后，计算出这些生成粒子的形状参数，以便与它们的自然对应物进行直接比较。图15直观展示了通过SH重建的真实沙粒与其生成对应物之间的对比。所提出的方法能够真实再现自然沙粒的典型轮廓特征——即由伸长率指标决定的特征，以及包括中观尺度角度特征和微观尺度表面纹理在内的多尺度形态本质。该图中的成对颗粒存在视觉上的差异，这是因为本研究开发的算法并非旨在精确复制单个真实沙粒的几何形状，而是依据天然沙子的形态统计分布生成全新的随机颗粒样本。每个生成的单一颗粒都能准确匹配目标形状参数，且颗粒集合的整体参数分布与真实沙粒高度一致。下载：下载高分辨率图像（597KB）下载：下载全尺寸图像图15. SH重建的天然沙粒（上）与所提形态学算法生成的颗粒（下）的对比。为定量验证统计一致性，我们首先检测生成的颗粒是否在集中趋势上引入了系统性偏差。针对每个形状描述符（EI、AI、RO），使用独立样本t检验公式计算原始颗粒集与生成颗粒集均值之间的差异（公式(34)）。(34)t=μgen-μorigvgen2ngen+vorig2norig其中μ表示样本均值，v2表示样本方差，n表示样本量。α表示显著性水平，计算出的t统计量绝对值越小，说明两组之间的相似度越高。由于每组的样本量均为120（属于大样本），t分布可近似为标准正态分布。在α=0.05的显著性水平下，双尾检验的临界值约为1.96。因此，若计算出的t值绝对值小于1.96，则无法拒绝零假设，即认为两组均值之间不存在统计学上的显著差异（Liu和Shao，2014）。此外，还通过比较虚拟颗粒与天然沙子的全形态分布累积分布函数来评估其再现精度。累积分布函数能够完整描述形状参数的概率特征。对于某一形状描述符的给定数值，F(x)表示随机选取的颗粒其EI值小于或等于x的概率。与基本的统计矩（如均值和方差）不同，累积分布函数能够全面体现分布的集中趋势、离散程度及形状特征。因此，对整个数值范围进行比较，可以有效地评估典型特征及整体变异性的再现精度。不过，传统的分布比较方法（如Kolmogorov–Smirnov检验）（Nie等人，2020a）在量化累积分布函数之间的精确点差异方面存在局限。为解决这一问题，我们采用了Fan等人（2023）提出的决定系数R2。该指标用于衡量原始颗粒与生成颗粒的参数值有序序列之间的吻合程度。对于某一形状描述符，R2的计算公式如下：(35)R2=1-∑i=1N(Xgen(i)-Xorig(i))2∑i=1N(Xorig(i)-Xˉorig)2其中Xorig(i)和Xgen(i)分别表示原始累积分布函数与生成累积分布函数中相同累积概率处的参数值，Xorig(i)为原始值的均值，N为数据点数量。R2值为1时表示完全吻合，数值越接近1，说明对分布的再现精度越高。原始颗粒与生成颗粒的形状参数分布的定量比较结果总结于图16和表2。如表2所示，所有计算出的t值在95%置信水平下均低于1.96的临界值，这说明对于三种形状描述符而言，原始颗粒集与生成颗粒集的均值之间不存在统计学上的显著差异。下载：下载高分辨率图像（353KB）下载：下载全尺寸图像图16. 原始颗粒与随机生成的沙粒的形状参数累积分布：(a) EI，(b) AI，(c) RO。表2. 形状描述符的决定系数R2及t检验结果。形状描述符R2tEI0.9820.614AI0.970?0.771RO0.993?0.283图16展示了120个真实母颗粒及其生成的模拟颗粒的形态参数累积分布函数。生成颗粒的受控形态指标的累积分布函数与真实母颗粒的累积分布函数高度吻合，这从定量角度证明了所提方法能够有效生成与真实沙粒形态一致的多样化沙粒。在本研究中，采用了一种稳健且具有普遍适用性的方法，即运用平衡的适应度函数与保守的变异策略。所提方法所达到的高一致性，充分证明了其在再现真实沙粒形态特征方面的有效性。5. 应用5.1. 在DEM模拟中的应用为验证所提颗粒生成方法在DEM模拟中的应用潜力，进行了多组颗粒柱坍塌试验。每种形状参数选取了三个数值进行模拟测试，共计七种情况，具体见表3。表3. DEM颗粒柱坍塌试验中各颗粒组的形状参数。颗粒组EIAIROⅠ0.50.050.01Ⅱ0.70.050.01Ⅲ0.90.050.01Ⅳ0.70.10.01Ⅴ0.70.150.01Ⅵ0.70.10.02Ⅶ0.70.10.03该颗粒模型是在3D颗粒流代码（PFC3D）中通过团块技术实现的，其中构成球体的数量由角分辨率φ和最小半径与最大半径之比ρ决定。为在保持计算效率的同时提升形状精度，角分辨率设定为φ=150°，半径比设定为ρ=0.2。各颗粒组的团块模型如图17(a)所示。研究发现，颗粒的棱角程度和表面粗糙度越高，团块模型所需的填充球体数量就越多。DEM模拟中使用的关键物理参数总结于表4，这些数值取自Wei等人（2018）的研究。图17(b)展示了颗粒在自然重力作用下的初始状态。每组包含1000个颗粒，这些颗粒生成在一个尺寸为15毫米×15毫米×120毫米的矩形柱体内，颗粒直径在1毫米到3毫米之间均匀分布。系统达到平衡后，立即移除容器右侧的壁面，以此模拟颗粒柱的坍塌过程。下载：下载高分辨率图像（406KB）下载：下载全尺寸图像图17. (a) 各组的团块结构；(b) 颗粒柱在自然重力作用下的初始状态。表4. DEM模拟中使用的物理参数。参数数值密度（千克/立方米）2600重力加速度（米/秒2）9.8颗粒直径（毫米）1.0–3.0摩擦系数0.5阻尼比0.7颗粒法向与剪切刚度（牛顿/米）2×10?墙体法向与剪切刚度（牛顿/米）2×10?图18展示了各组颗粒坍塌后的位移云图及最终坡度。显然，颗粒的棱角程度/粗糙度越高，伸长指数越低，坍塌后的坡度角度就越陡。粗糙度和伸长程度的影响尤为显著，案例I的坡度角度高达54.78°，案例VII的坡度角度为56.27°。这种现象的原因是，圆形光滑的颗粒之间容易滑动，因此形成的坡度较缓，而棱角分明的颗粒之间相互咬合紧密，从而形成更陡且更稳定的结构。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图18. 不同形状参数下颗粒坍塌后的位移云图：(a) EI变化情况，(b) AI变化情况，(c) RO变化情况。图19展示了各颗粒组摩擦活化指数Im的概率分布，该指数定义为接触点处切向力ft与最大滑动阻力μfn的比值。结果表明，所有样本的Im值都集中在0.9至1.0之间。不过，随着AI和RO值的增加以及EI值的减小，Im值落在该区间内的概率逐渐上升。这一结果表明，当颗粒的伸长程度更高、棱角更明显且表面更粗糙时，更多接触点会处于高滑动阻力状态。这进一步体现了颗粒形状对接触摩擦行为及咬合作用的直接影响。下载：下载高分辨率图像（323KB）下载：下载全尺寸图像图19. 不同形状参数下颗粒组摩擦活化指数Im的概率分布：(a) EI变化情况，(b) AI变化情况，(c) RO变化情况。5.2. 用于3D打印透明土壤的制备本研究提出的算法可有效应用于3D打印透明土壤试样的制备。为验证其可行性，首先将算法生成的具有多样化形态的颗粒模型导入Formlabs PreForm软件中进行预处理，包括模型方向优化、支撑结构生成以及切片参数校准。随后，使用Form 4打印机（图20(a)）以及Clear Resin V5材料打印这些颗粒，层厚设定为0.025毫米，以确保最佳的几何精度。打印完成的颗粒经过标准的后处理工序，包括残留树脂清除、支撑结构去除以及紫外线固化，最终得到具有稳定物理性能的完整试样（图20(b)）。在25摄氏度条件下，打印颗粒材料的折射率被测定为1.5355，这一数值是匹配孔隙流体折射率的关键参考值——而匹配孔隙流体折射率是制备透明土壤的基本前提。根据Arago-Biot体积加权折射率混合规则（Pretorius等人，2021），通过按4.68:1的体积比例混合苯甲基硅油（20摄氏度时的折射率为1.552）和15号白油（20摄氏度时的折射率为1.460），制备出了折射率匹配的混合孔隙流体。下载：下载高分辨率图像（314KB）下载：下载全尺寸图像图20. (a) Formlabs Form4 3D打印机，(b) 具有不同形态的3D打印沙粒，(c) 经过激光照射后的3D打印透明土壤试样的散斑场。透明土壤试样的制备步骤如下：首先，将两种油按照规定的体积比例充分混合，然后使用阿贝折射仪精确调整混合液的折射率，使其与打印颗粒的折射率高度匹配。接着，将孔隙流体倒入一个10×10×10厘米的亚克力容器中，再用筛子将3D打印的颗粒缓慢放入容器中，同时用玻璃棒轻轻搅拌。之后，将装有颗粒与流体混合物的整个容器放入真空室中，抽气2小时，以彻底消除残留的微气泡。最后，让试样静置，这样就得到了具有良好光学透明度的透明土壤样本。通过激光片照射来评估所制试样的光学性能。实验中使用的是输出功率为4瓦、厚度约为1毫米的532纳米连续波聚焦绿光激光片，用该激光片扫描亚克力容器，即可清晰地观察到容器内部的截面结构（图20(c)）。结果显示，容器内部截面呈现出均匀分布的散斑场，这进一步证明了使用3D打印合成颗粒制备的透明土壤试样具有令人满意的光学性能，同时也证实了所生成颗粒在物理建模方面的实际应用价值。6. 结论本研究提出了一种两阶段流程模型，可用于生成具有精确且可独立控制的多尺度形态特征（EI、AI、RO）的沙粒。该模型通过分形控制的傅里叶谱合成方法构建多尺度纹理表面，再通过高斯核变异遗传算法进一步优化颗粒几何形状。通过累积分布函数比较及统计检验，从定量角度验证了该方法与天然沙子的形态一致性，同时通过DEM模拟和3D打印透明土壤试验验证了该方法的实用性。主要结论如下：(1) 构建了一个用于生成随机沙粒的集成多尺度生成框架。分形控制的傅里叶谱合成方法与遗传算法的结合，能够生成逼真的初始颗粒表面。在二十面体网格上的优化避免了几何变形，而高斯粗糙化与平滑化策略则实现了对形态特征的独立调控。设计的权重函数能够平衡多尺度参数的影响，从而在确保几何连贯性和计算效率的同时实现高精度的优化。(2) 所提模型能够实现与天然沙子高度一致的形态特性及统计一致性。目标形态参数是根据真实粗砂和细砂颗粒的X射线μCT扫描结果确定的。生成的颗粒的EI、AI和RO的累积分布函数与真实沙粒的高度一致，所有的决定系数均大于0.97。独立的t检验进一步证实，在95%的置信水平下，两种颗粒的形态均值之间不存在显著的统计学差异。(3) 所提方法已在DEM模拟和透明土壤实验中得到验证。在颗粒柱坍塌试验中，颗粒的棱角程度和粗糙度越高，伸长程度越低，最终形成的坡度就越陡，摩擦活化指数也越高，这些结果与理论预期相符。此外，高精度的3D打印颗粒还被用于制备具有良好透光性且内部结构清晰可见的透明土壤，这进一步证明了所生成颗粒在物理建模中的可行性。不过，当前方法也存在一些局限性。首先，由于形态参数和几何约束是基于球面谐波表示方法建立的，所以所提方法仅适用于星形颗粒，而不适用于具有复杂凹面结构的非星形颗粒。其次，遗传算法在每一代都需要过滤掉大量中间几何形态，这就导致了频繁的无效计算，进而增加了计算成本。第三，目前的随机突变操作缺乏精确的维度控制，难以独立调整某一形态属性。为今后改进，可开发定向目标突变策略，以便在特定尺度上独立调节某一形态参数，例如仅修改表面粗糙度而不影响弹性模量和抗压强度。该方案有望减少冗余迭代次数，加快优化算法的收敛速度。在工程应用方面，该方法可为颗粒材料的离散元模拟提供形态可控的数字粒子样本，支持对沙土的介观力学研究及典型岩土工程的数值分析，还能用于制造定制化的3D打印实物颗粒以用于实验室测试，同时可为各类工程颗粒材料的形态设计与性能优化提供参考。

CRediT作者贡献说明：
周鹏：写作——审阅与编辑、指导、项目管理、资金获取、概念构思。万俊杰：写作——初稿撰写、验证、方法论、实验研究、正式分析、数据整理。徐长杰：指导、项目管理、资金获取。丁宣明：指导、概念构思。童龙勇：写作——审阅与编辑、指导、方法论、实验研究、概念构思。

周鹏|万俊杰|徐长杰|丁宣明|童龙勇
中国南昌市330013，华东交通大学土木建筑学院，山区土木工程安全与韧性国家重点实验室