瑞雷波是地震勘探中常用的地震波类型之一，利用其在层状介质中的频散特性能够获得地下介质的信息，如地层厚度、岩性、介质的变化等^[1-4]。瑞雷波勘探的流程通常包括数据的采集、处理和解释，瑞雷波频散曲线的提取是数据处理中的关键步骤，其对后续的地层结构反演起到至关重要的作用^[5]。通过频散曲线可以将瑞雷波的相速度转化为横波速度，进而实现对地下结构的速度分层，这在工程地质勘查^[6-8]、环境探测^[9]、城市地下空间探测^[10] 等众多领域得到了广泛应用。

瑞雷波频散曲线反演为典型的多参数、多极值、高度非线性的地球物理优化问题，频散曲线的高精度、高效率反演方法为重点研究内容。常用的频散曲线反演方法主要分为局部线性化方法和全局非线性方法。局部线性化方法中最常用的是最小二乘法^[11]，其在初始模型建立不够精准的情况下存在难以通过反演获得准确地层模型结构的不足。以模拟退火算法（Simulated Annealing，简称SA）^[12]、遗传算法（Genetic Algorithm，简称 GA）^[13]、蒙特卡罗法（Monte Carlo，简称 MC）^[14]、群体智能优化算法^[15-17] 等为代表的全局非线性方法对初始模型依赖性较小，可进行全局寻优，是反演方法的主要发展方向。具体到瑞雷波频散曲线反演方法，其受样本规模的影响，收敛速度不确定，尤其在面对复杂地质模型时，还需合理设置初始参数，避免反演搜索过早陷入局部极小值^[18]。

针对上述问题，LU 等对模拟退火算法进行改进，提出了热浴模拟退火算法（Heat-bath SA），并对理论以及实际地震记录进行反演^[19]。SOCCO 等提出了一种改进的蒙特卡罗（Improved Monte Carlo，简称 IMC）算法，利用瑞雷波正演无量纲化的特性对模型空间的采样过程进行优化，通过统计以得到最终反演结果^[20]。LÜ 等提出了改进的粒子群（Im‐ proved Eliminate Particle Swarm Optimization，简称 IEPSO）算法，通过加入局部-全局信息共享项以提高算法的全局搜索能力^[21]。吴聪等提出了改进自适应遗传算法（New Improved Adaptive Genetic Algo‐ rithm，简称 NIAGA），将原算法中的交叉概率和变异概率进行自适应化，使其能够根据个体适应度、迭代次数等参数的变化进行自我调整，达到抑制算法早熟收敛、提高全局探索的目的^[22]。

蜣螂优化（Dung Beetle Optimizer，简称 DBO） 算法是 2022 年提出的一种全新的群体智能优化算法^[23]。该算法受蜣螂生活习性的启发，通过模拟蜣螂种群的滚球、跳舞、产卵、觅食、偷窃 5 种行为模式，实现对空间的高效探索与开发。DBO算法含参数较少，对模型参数初始设置没有特殊要求，该算法通用性较强，在处理复杂高维优化问题时，具备高收敛速度、高精确度、高稳定性的特点。然而，DBO算法在种群初始化的处理中，可能存在因种群聚集而降低算法收敛性能的不足。为了提高初始化种群的质量，保证频散曲线反演的效率和准确性，笔者提出基于改进DBO算法的瑞雷波频散曲线反演方法，该方法采用 Halton序列初始化种群的位置，使种群的初始分布更加均匀；然后将种群分为4 种不同的子种群，每个子种群对应不同的搜索策略，避免了反演搜索陷入局部最优解，实现了算法的快速收敛。首先简要介绍瑞雷波频散曲线反演目标函数以及 DBO 算法的基本原理；然后针对 DBO算法初始化种群分布不均匀的缺点，提出基于改进 DBO 算法的瑞雷波频散曲线反演方法的处理流程；最后，使用3种典型的理论地质模型和实际资料对方法的有效性进行了验证。

1 方法原理

1.1 瑞雷波频散曲线反演目标函数

瑞雷波频散曲线反演本质上是通过某种策略不断更新地层模型参数，使目标函数最小化的过程。传统瑞雷波频散曲线反演目标函数一般设定为均方差函数，其表达式为：

此外，反演目标函数也可以通过实测与理论计算的频散曲线值的相对误差计算得到，其表达式为：

该目标函数受所有采样频率点信息的共同影响，经过多次试验测试，采用 25～35 个频率点即可获得较为准确的反演结果，且可以保证计算效率。本文采用由（2）式确定的瑞雷波频散曲线反演目标函数F₂。

1.2 DBO算法

蜣螂是自然界常见的昆虫，在自然界中扮演着分解者的角色。蜣螂依赖粪球生存，根据蜣螂个体对粪球利用方式的不同，可以将蜣螂种群划分为滚球蜣螂、产卵蜣螂、小蜣螂、窃贼蜣螂4个子种群，不同的蜣螂子种群对应不同的行为模式。

1.2.1 滚球蜣螂行为模式

滚球蜣螂负责将粪便滚成球，运送到安全的地方。滚球蜣螂在滚动粪球的过程中，可以实现对全空间的探索。滚球蜣螂的行为模式有滚球行为和跳舞行为2种。

当前方无障碍时，实施滚球行为，其移动路径受到光源强度的影响，位置更新的计算公式为：

当前方出现障碍时，实施跳舞行为，获取新的滚球方向，位置更新的计算公式为：

1.2.2 产卵蜣螂行为模式

产卵蜣螂负责将粪球滚到适合产卵的安全区域，产下卵球以繁殖后代。产卵蜣螂的移动范围被严格限制在安全区域内，实现对区域的局部开发。当产卵蜣螂处于安全区域时，实施产卵行为，位置更新的计算公式为：

需要注意，产卵蜣螂的位置被严格限制在安全区域内。若X^t_spawning（ j）>U^∗ _b（ j ），则 X^t _spawning（ j）= U^∗ _b（ j）；若X^t _spawning（ j ）<L^∗ _b（ j ），则X^t _spawning（ j）= L^∗ _b（ j）。

1.2.3 小蜣螂行为模式

小蜣螂需要寻找食物，其活动范围被严格限制在最佳觅食区。当小蜣螂处于该区域时，实施觅食行为，位置更新的计算公式为：

其中：

1.2.4 窃贼蜣螂行为模式

窃贼蜣螂会与其他个体竞争食物，偷取粪球。窃贼蜣螂位于最佳食物源附近时，实施偷窃行为，位置更新的计算公式为：

综上所述，蜣螂优化算法正是通过对蜣螂种群的滚球、跳舞、产卵、觅食和偷窃5种行为进行建模，然后基于 4 个子种群执行不同的搜索方法，从而进行全局探索和局部开发。

1.3 DBO算法改进

DBO算法的种群初始化采用随机法生成，无法保证种群中的个体均匀分布，可能导致种群搜索速度慢，算法难以收敛。针对上述问题，引入 Halton 序列初始化种群个体的位置，可以更好地控制个体的分布情况，避免在搜索空间的某些区域过于密集或过于稀疏的情况，使得算法更加完善和稳定；同时，由于 Halton 序列具有较好的遍历性，能够更好地覆盖整个搜索空间，因此可以更快地发现全局最优解，提高算法的收敛速度和精度。

二维Halton序列的实现过程为选定2个不同的质数 p1和 p2，对于任意的正整数 n，分别以 p₁和 p₂作为基础量不断对其切分和组合，以得到二维 Halton 序列H （ n），切分过程可以描述为：

为验证 Halton 序列初始化的优越性，取 Halton 序列产生的初始种群与随机方法产生的初始种群进行对比。种群个数取 100，取值范围为[ 0，1]。由随机方法和 Halton序列初始化种群分布对比（图1） 可以看出，Halton序列产生的种群分布更加均匀，初始解质量更高。

2 方法流程

基于改进 DBO 算法提出一种新的频散曲线反演方法，具体流程（图2）为：①设定待反演模型的基本结构，包括采样频率 $f_M={\left[f_1，f_2，\cdots ，f_M\right]}^{\mathrm{T}}$，地层数量 d，地层横波速度${\mathit{V}}_{\mathrm{s}}={\left[V_{{\mathrm{s}}_1}，V_{{\mathrm{s}}_2}，\cdots ，V_{{\mathrm{s}}_4}\right]}^{\mathrm{T}}$，地层纵波速度$V_{\mathrm{P}}={\left[V_{{\mathrm{P}}_1}，V_{{\mathrm{P}}_2}，\cdots ，V_{{\mathrm{P}}_d}\right]}^{\mathrm{T}}$，地层密度$\rho ={\left[{\rho }_1，{\rho }_2，\cdots ，{\rho }_d\right]}^{\mathrm{T}}$，地层厚度$\mathit{H}={\left[H_1，H_2，\cdots ，H_d\right]}^{\mathrm{T}}$，以及横波速度的搜索范围U_b ^VS 和L_b ^VS 与地层厚度的搜索范围 U_b ^H 和 L_b ^H。初始化 DBO 算法参数，包括蜣螂种群总数 N，空间搜索维度 D = 2d − 1，最大迭代次数T_max；根据具体问题，确定滚球蜣螂、产卵蜣螂（卵球）、小蜣螂、窃贼蜣螂的数量以及权重系数 k，b，S 等。②利用 Halton 序列初始化蜣螂种群的位置X⁰，每只蜣螂个体对速度模型进行正演模拟，采用快速矢量传递算法^[24] 计算得到瑞雷波相速度的理论值 c_cal，进而计算该值与实测值 c_obs 的目标函数 F₂。③利用 DBO 算法更新横波速度-地层厚度模型，具体步骤为：1）根据计算得出的种群中每只蜣螂个体的目标函数值，得到最小目标函数值对应的全局最优位置 X^best 和最大目标函数值对应的最差位置 X^worst。 2）更新滚球蜣螂的位置。若该蜣螂个体所属种类为滚球蜣螂，通过概率法确定该蜣螂个体接下来的行为继续是滚球还是跳舞，给定一个滚球概率阈值 λ，令 δ = rand（1），若 δ <0.9，则令 η = rand（1），α = 1（ η >λ）或 α =− 1（ η ≤ λ），并按照（3）式更新位置，否则按照（4）式更新位置。进而确定当前局部最优位置X^∗。3）更新其他3个子种群蜣螂的位置。若该蜣螂个体所属种类为产卵蜣螂，则按照（5）式更新位置；若该蜣螂个体所属种类为小蜣螂，则按照（9） 式更新位置；若该蜣螂个体所属种类为窃贼蜣螂，则按照（12）式更新位置。4）计算当前位置蜣螂个体的目标函数值，得到最小目标函数值对应的全局最优位置X^best。5）判断是否满足迭代终止条件。若种群中最优蜣螂个体目标函数值F_b <∈，或者迭代次数达到T_max时，迭代终止，输出全局最优位置X^best（横波速度和地层厚度）；否则，跳转至步骤2）进行下一次迭代寻优。④输出最小目标函数值对应的横波速度模型结构。

3 理论模型数据测试

为了说明改进 DBO算法的有效性，对 3个理论地层模型（即四层速度递增模型、四层含低速软夹层模型和四层含高速硬夹层模型）进行瑞雷波频散曲线的反演测试，并与 NIAGA 算法在搜索效率和搜索精度上进行对比。为了达到兼顾反演准确度和运算时长的目的，经过多次测试，设定所有算法的蜣螂种群总数N为210，最大迭代次数T_max为150。改进 DBO 算法的参数设置如下：滚球蜣螂、产卵蜣螂、小蜣螂以及窃贼蜣螂的数量比例为 6∶6∶7∶11；k 取值为0.1，b取值为0.3，S取值为0.5。NIAGA算法的参数设置如下：交叉概率 p_c，p_c1，p_c2的取值分别为 0.6，0.6，0.3；变异概率 p_m，p_m1，p_m2 的取值分别为 0.01，0.01，0.002。在测试中，对于模型Ⅰ和模型Ⅲ，在 5～99 Hz 频率范围内计算基阶频散曲线并进行反演；对于模型Ⅱ，考虑实际情况，仅在5～60 Hz频率范围内计算基阶频散曲线并进行反演。

3.1 模型Ⅰ（四层速度递增模型）

四层速度递增模型在勘探中最为常见，地质参数和反演搜索范围如表1所示。对提取的频散曲线利用改进 DBO 算法和 NIAGA 算法分别进行 50 次独立运算，由频散曲线反演结果（图3）可以看出，2 种算法计算得到的频散曲线与理论值均无明显偏差；但通过与频散曲线理论值进行统计分析发现，改进 DBO 算法的结果与理论值的平均相对误差为 0.01%，而 NIAGA 算法的结果与理论值的平均相对误差为 0.11%。2 个算法的目标函数随迭代次数的收敛曲线（图3c）显示，改进DBO算法的收敛过程更加连续且稳定，在第 50 次迭代时已经趋于收敛，且速度更快，同时也能获得比 NIAGA 算法更加精确的解。

通过改进 DBO 算法得到的地层横波速度结构集中地分布在真实地层周围（图4a）。对比 NIAGA 算法得到的地层横波速度结构（图4b）可以看出， NIAGA 算法的结果与真实模型存在较为明显的差异。统计分析50次反演结果（表2）发现，改进DBO 算法计算的地质参数与真实模型的平均相对误差为0.76%，显著低于NIAGA算法的8.60%，说明改进 DBO算法所得结果更为准确。

3.2 模型Ⅱ（四层含低速软夹层模型）

模型Ⅱ与人工铺就路面的地下横波速度结构相似，其地质参数和反演搜索范围如表3 所示。对提取的频散曲线利用改进DBO算法和NIAGA算法分别进行 50 次独立运算，所得频散曲线反演结果 （图5）显示，改进 DBO算法的结果可以准确刻画反映低速软夹层的“之”字形结构，而 NIAGA 算法的结果在曲线形态上虽呈现出反映低速夹层的“之” 字形结构，但仍能观察到反演结果在 17～49 Hz 频段处与理论值存在一定程度的偏差。通过与频散曲线理论值统计分析发现，改进DBO算法的结果与理论值的平均相对误差为0.03%，而NIAGA算法的结果与理论值的平均相对误差为1.06%。由目标函数随迭代次数的收敛曲线（图5c）可以看出，2 种算法都能快速收敛，改进 DBO 算法比 NIAGA 算法具有更高的收敛速度，NIAGA 算法收敛至 0.2 左右时即无法进一步收敛，而改进 DBO 算法可收敛至 0.08，具有更高的求解精度。

通过改进 DBO 算法得到的地层横波速度结构可以准确地反映出地层含有低速夹层的特征，且各独立运算结果也比较接近真实值（图6a）。对比NI‐ AGA 算法得到的地层横波速度结构（图6b）可以看出，NIAGA算法的反演结果虽能反映出地层含有低速夹层的特征，但并未准确计算出第 1层、第 2层的横波速度以及第 2 层、第 3 层的地层厚度。统计分析50次反演结果（表4）可以看出，改进DBO算法计算的地质参数与真实模型的平均相对误差为 1.30%，显著低于 NIAGA 算法的 10.27%，说明改进 DBO算法所得结果更为准确。

注：V_S1，V_S2，V_S3，V_S4的单位为m/s；H₁，H₂，H₃的单位为m。

3.3 模型Ⅲ（四层含高速硬夹层模型）

模型Ⅲ为四层含高速硬夹层模型，其地质参数和反演搜索范围如表5所示。对提取的频散曲线利用改进 DBO 算法和 NIAGA 算法分别进行 50 次独立运算，所得频散曲线反演结果（图7）显示，2 种算法计算得到的频散曲线与理论值均具有较高的拟合度，通过与频散曲线理论值统计分析发现，改进 DBO 算法的结果与理论值的平均相对误差为 0.02%，而 NIAGA 算法的结果与理论值的平均相对误差为 0.78%（表6）。由目标函数随迭代次数的收敛曲线（图7c）可以看出，NIAGA算法虽然搜索效率很高，但在第 30 次迭代后陷入早熟状态，改进 DBO算法可以收敛至全局最优解附近，且收敛过程更加连续稳定。

通过改进 DBO 算法得到的地层横波速度结构可以准确地反映地层包含高速夹层的特征，且各独立运算结果比较接近真实值（图8a）。对比 NIAGA 算法得到的地层横波速度结构（图8b）可以看出， NIAGA 算法仅对第 1 层的横波速度和地层厚度的计算结果较为准确，对第 2 层及以下地质参数的计算结果与真实值存在较为明显的偏差。统计分析 50次反演结果（表6）发现，改进DBO算法计算的地层参数与真实模型的平均相对误差为 1.51%，显著低于 NIAGA 算法的 7.07%，说明改进 DBO 算法所得结果更为准确。

注：V_S1，V_S2，V_S3，V_S4的单位为m/s；H₁，H₂，H₃的单位为m。

总之，以上 3 个模型的频散曲线反演结果分析表明，改进DBO算法可以有效反演地层的横波速度

4 实际资料测试

为了说明提出的方法在实际资料频散曲线反演中的适用性，对冰岛西南部 Arnarbæli地区（图9） 的地震数据进行了处理和分析。该地区地表以时代较新（距今地质年龄小于 0.7 Ma）的溢流玄武岩和火山玻璃角砾岩互层为主^[25]。地震数据由质量为 6.3 kg的重锤锤击地面激发作为主动震源，24个4.5 Hz的垂向检波器以1 m的间隔、最小炮检距为10 m 的线性排列构成观测系统采集得到。图9a 为采集到的时长为0.8 s的地震记录，采样间隔为1 ms。对该地震记录利用相移法得到瑞雷波频散能量图（图9b），可以观察到明显连续的基阶频散能量。

在实际面波勘探中，通常最大勘探深度为瑞雷波最大波长的三分之一至二分之一。在对该实际资料频散曲线的反演过程中，以8～40 Hz的基阶频散曲线为研究对象，并选择12 m作为本次反演的最大勘探深度，反演最大勘探深度以上地层的横波速度和地层厚度以及模型的纵波速度和密度设为合理的固定值。由于缺少Arnarbæli地区的测井数据，将 Olafsdottir等反演得到的横波速度结构作为基准值^[26]，据此将勘探空间划分为6层，各层横波速度和地层厚度的反演搜索范围、纵波速度和密度的设定如表7 所示。

使用改进DBO算法与NIAGA算法对提取的基阶频散曲线进行反演并对比结果，算法的参数设置与理论地层模型测试保持一致。改进 DBO 算法频散曲线反演得到的结果在整个频段都能很好地与基准值吻合（图10a），30 次独立运算结果的平均相对误差为 0.16%。NIAGA 算法频散曲线反演的结果在10 Hz以下与基准值的拟合度稍差（图10b），平均相对误差为0.31%。由目标函数随迭代次数的收敛曲线（图10c）可以看到，改进 DBO 算法在迭代初期，最小目标函数值即优于NIAGA算法；在第50次迭代后 2种算法都趋于稳定收敛，NIAGA算法最终收敛于 0.36，而改进 DBO 算法可以收敛至 0.06，具备更高的收敛精度，其反演结果更加可靠。

改进DBO算法反演得到的横波速度（图11a）显示，各独立运算结果的平均值与基准值接近。NIA‐ GA算法得到的横波速度结构（图11b）显示，与基准值的拟合情况并不理想，特别是第 3层、第 4层和第 5 层的地质参数计算结果偏差明显。2 种算法结果对比表明，改进 DBO 算法所得结果优于 NIAGA 算法，反演所得速度模型的可靠性更高。

注：V_S1，V_S2，V_S3，V_S4的单位为m/s；H₁，H₂，H₃的单位为m。

5 结论

提出一种基于改进 DBO 算法的频散曲线反演方法。该方法利用 Halton序列进行种群初始化，优化了种群的空间分布，对不同子种群设置不同的搜索策略及更新规则，有利于算法在局部和全局的搜索。通过3个典型的地质模型测试验证方法的有效性和稳定性。对实际资料的处理结果表明，基于改进蜣螂优化算法的瑞雷波频散曲线反演方法能够有效获得准确的地层横波速度，与 NIAGA 算法结果对比表明，改进 DBO 算法具有效率更高、结果更可靠的优点。

当然，改进 DBO 算法也存在一些不足，比如子种群的数量设置基于经验，实测数据计算耗时相对较长，对这些问题进行深入研究将大幅提高算法的求解能力。

符号解释

b——滚球系数，常数，且b ∈（0，1）

b₁，b₂——产卵系数，为1 × D维独立随机向量；

c_cal ——瑞雷波相速度的理论值，m/s；

c_cal（ m）——瑞雷波相速度的理论值的第m个分量，m/s；

c_obs（ m）——瑞雷波相速度的实测值的第m个分量，m/s；

C₁——觅食系数，服从正态分布的随机数；

C₂——觅食向量，属于（ 0，1）的1 × D维随机向量；

d——地层数量；

D——空间搜索维度；

f_M——第M个采样频率，Hz；

f_M——采样频率，Hz；

F₁——传统瑞雷波频散曲线反演目标函数；

F₂——本文采用的瑞雷波频散曲线反演目标函数；

F_b——蜣螂种群中最优蜣螂个体目标函数值；

g——偷窃向量，服从正态分布的1 × D维随机向量；

H——地层厚度，m；

H_d——第d层地层厚度，m；

H （ n ）——二维Halton序列；

j——蜣螂个体位置序号；

k——偏转系数，常数，且k ∈（ 0，0.2 ]；

L_b——搜索空间的下界，维度为D；

L_b^best ——觅食区域的下界，维度为D；

L_b ^H——地层厚度搜索范围的下界，m；

L_b ^VS ——横波速度搜索范围的下界，m/s；

L_b^∗——安全区域的下界，维度为D；

L_b ^∗（ j）——安全区域的下界的第j个分量；

m——采样频率序号；

M——频散曲线的频率点总数；

n——任意正整数；

N——蜣螂种群总数；

p，p₁，p₂——大于等于 2 的质数，表示 Halton 序列基础量；

p_c，p_c1，p_c2——交叉概率；

p_m，p_m1，p_m2——变异概率；

p_r，p^-r-1 ——Halton序列切分项；

r——由n和p决定的常量；

R——区域控制系数；

S——偷窃系数，常数；

t——当前迭代次数；

T_max——最大迭代次数；

U_b——搜索空间的上界，维度为D；

U_b^best ——觅食区域的上界，维度为D；

U_b ^H——地层厚度搜索范围的上界，m；

U_b ^VS ——横波速度搜索范围的上界，m/s；

U_b^∗（ j）——安全区域的上界的第j个分量；

V_P——地层纵波速度，m/s；

V_Pd ——第d层纵波速度，m/s；

V_S——地层横波速度，m/s；

V_Sd ——第d层横波速度，m/s；

X_best ——全局最优位置向量，维度为D；

${\mathit{X}}_{\text{molling}}^{t-1}$——第 t − 1 次迭代时滚球蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{molling}}^t$——第t次迭代时滚球蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{rolling}}^{t+1}$——第 t + 1 次迭代时滚球蜣螂个体的位置向量，维度为D；

$\left|{\mathit{X}}_{\text{rolling}}^t-{\mathit{X}}_{\text{molling}}^{t-1}\right|$——位置差异信息，其为滚球蜣螂个体在第 t次迭代与第 t − 1次迭代时的位置差异；

$\left|{\mathit{X}}_{\text{malling}}^t-{\mathit{X}}^{\text{unest}}\right|$——光源强度信息，其为模拟自然环境中光源强度的变化，数值越高代表光源强度越弱；

${\mathit{X}}_{\mathrm{s}\mathrm{m}\mathrm{a}\mathrm{l}}^t$ ——第 t次迭代时小蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{small}}^{t+1}$——第 t + 1 次迭代时小蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{spaunings}}^t$——第 t 次迭代时产卵蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{spanaing}}^t（j）$——第t次迭代时产卵蜣螂个体的第j个位置分量；

${\mathit{X}}_{\text{spauning}}^{t+1}$——第 t + 1 次迭代时产卵蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{thief}}^t$——第 t 次迭代时窃贼蜣螂个体的位置向量，维度为D；

${\mathit{X}}_{\text{thief}}^{t+1}$——第 t + 1 次迭代时窃贼蜣螂个体的位置向量，维度为D；

X^worst ——全局最差位置向量，维度为D；

X^∗ ——当前局部最优位置向量，维度为D；

X⁰ ——蜣螂种群初始位置；

y₀， y₁，y_r ——常数，Halton序列切分项系数；

α——自然环境系数，整数，取值为-1 或 1，分别代表偏离和未偏离原来方向；

δ——滚球概率系数，属于（ 0，1）的随机数；

ϵ——允许误差阈值，常数；

η——滚球偏转概率系数，属于（ 0，1）的随机数；

θ——偏转角度，且θ ∈ [ 0，π]；

ρ——地层密度，g/cm³；

ρ_d——第d层地层密度，g/cm³；

λ——滚球概率阈值，常数；

Φ_p （n）——以p为基础量的Halton序列函数；

Φ_p1（ n ）——以p₁为基础量的Halton序列函数；

Φ_p2 （ n ）——以p₂为基础量的Halton序列函数。

参考文献