CN116822189B - 一种考虑地表耦合的可控震源振动器输出力的计算方法 - Google Patents

Abstract

本发明公开了一种考虑地表耦合的可控震源振动器输出力的计算方法，包括如下：构建振动器平板‑大地接触振动模型，采用表面分形模型构建大地表面形貌；基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组的求解问题；利用迭代法对线性方程组进行求解；基于最小二乘理论求解平板与大地之间接触力的关系式；建立考虑地表形貌的可控震源振动器‑大地耦合振动动力学模型；基于伽辽金法对可控震源振动器‑大地耦合振动动力学模型进行离散并利用四阶Runge‑Kutta法对其进行求解，获得振动输出力。本发明解决了现有可控震源振动器‑大地系统模型无法真实反映振动器与大地之间的相互作用关系的问题，提高了模型的适用性和精确性。

Description

一种考虑地表耦合的可控震源振动器输出力的计算方法

技术领域

本发明涉及油气勘探技术领域，具体涉及一种考虑地表耦合的可控震源振动器输出力的计算方法。

背景技术

随着油气勘探的主战场逐步向深层、非常规领域转移，勘探地表条件的日趋复杂化，对可控震源的勘探性能提出了更高的要求。由于可控震源的激发原理是通过振动器与大地之间的相互作用产生振动输出力，振动器-大地相互作用性质直接决定了振动器的动态响应，影响了地震波的激发效果，关系着可控震源的勘探性能。振动器-大地相互作用与地表性质息息相关，但现有可控震源振动器输出力计算模型对地表性质影响的考虑还远远不够，难以反映振动器与大地之间的真实相互作用关系，导致可控震源的实际输出力计算精度较低，难以满足可控震源在深层、非常规油气领域的应用。

针对上述问题，研究者们虽然意识到了地表性质对可控震源振动器输出力计算模型准确性的影响，但在实际建模过程中对耦合地表的考虑仍然比较缺乏，未能根据可控震源振动器的工作特点找到其关键影响因素。基于现有可控震源振动器输出力计算模型难以满足可控震源在深层、非常规油气领域的应用，亟需提出一种考虑地表耦合的振动器输出力计算方法，增强模型的适用性和计算精度，以符合工程实际，进而促进可控震源在深层、非常规油气领域的应用。

发明内容

为解决现有技术中存在的问题，本发明提供了一种考虑地表耦合的可控震源振动器输出力的计算方法，解决了现有可控震源振动器输出力计算模型未考虑地表形貌、地表耦合效应等的影响，难以满足工程实际的问题，提高了模型的计算精度和适用性。

为实现上述目的，本发明提供如下技术方案：一种考虑地表耦合的可控震源振动器输出力的计算方法，包括如下步骤：

S101、构建振动器平板-大地接触振动模型，将振动器平板与大地的接触等效成一个刚性平板和一个弹性粗糙表面的接触；

S102、采用表面分形模型构建大地表面形貌；

S103、基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组的求解问题；

S104、利用Gauss-Seidel迭代法对线性方程组进行求解，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移；

S105、基于最小二乘理论求解平板与大地之间接触力的关系式；

S106、利用薄板振动理论引入大地表面的耦合效应，采用Pasternak地基模型模拟深层大地的反应，建立考虑地表形貌的可控震源振动器-大地耦合振动动力学模型；

S107、基于伽辽金法对可控震源振动器-大地耦合振动动力学模型进行离散；

S108、利用四阶Runge-Kutta法对可控震源振动器-大地耦合振动动力学模型进行求解，获得振动输出力。

优选的，在步骤S101中，振动器平板-大地接触振动模型为单自由度的质量-弹簧-阻尼模型。

优选的，在步骤S102中，表面分形模型具有自相似和尺度独立的特性，能够从本质上描述粗糙表面的性质；弹性粗糙表面的大地表面形貌采用修正的两参数Weierstrass-Mandelbrot函数生成；所述修正的两参数Weierstrass-Mandelbrot函数的表达式为：

其中，L为样本长度，D为分形维数，G为分形粗糙度，M为生成分形表面的脊线数量，n为频率因子，n_max＝int[log(L/L_S)/logγ]_max，L_s为截断长度，Φ_m,n∈[0,2π]为随机相位，γ为缩放参数，γ>1；大地表面的粗糙度由D和G决定；表面分形维数D其物理意义为表示粗糙表面所占据的空间程度大小，D值越大对应于越密集的表面形态(越光滑的表面形貌)；分形粗糙度G是高度尺度参数，G值越大则表面形貌越粗糙。

优选的，在步骤S103中，基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组求解的问题，具体包括如下：

S301、设弹性材料的应变余能和总应变余能V^*的表达式分别为：

其中，为两个弹性体的应变余能，p为接触压力，/>和/>为两个接触体的预定位移，/>为基于几何干涉假定的接触区域内两接触体的总预定位移，/>为待计算接触区域内的复合表面位移，等于两个接触体总的接触变形；

根据弹性材料的应变余能在数值上等于弹性应变能U_E，将弹性材料的应变余能代入两个粗糙表面接触的总应变余能表达式里面，获得总应变余能V^*的表达式为：

S302、利用布森内斯解将表面位移和接触压力p联系起来，得到法向载荷作用下的接触压力p与表面位移/>之间的关系式为：

其中，E^*为复合杨氏模量，E₁,E₂分别为土壤和平板的杨氏模量，ν₁，ν₂分别为土壤和平板的泊松比；

S303、将表面位移代入接触压力p与表面位移/>之间的关系式，将求总应变余能最小值的接触问题转化为求解相应积分式关于接触压力p(x′,y′)的最小值问题，相应积分式为：

其中，为2a×2b矩形单元中心的表面位移，E^*为复合杨氏模量；

S304、将积分式在整个接触区域内进行离散，同时采用分段插值的方法获得每个离散单元的接触压力，从而将积分式转化为求和式；将积分式离散处理后获得的求和式为：

其中，为2a×2b矩形单元中心的表面位移，Ω_k为单元接触面积子域，c_num为根据平板与大地之间的刚体位移求得的几何干涉区域内的初始接触点，索引k和索引l分别代表接触压力产生的位置和表面位移发生的位置，而C_kl为影响矩阵，其表达式为：

S305、将接触压力向量p与表面位移之间的关系式离散化处理后，求得总应变余能V^*为：

将总应变余能V^*化为矩阵式：

其中,p^T＝[p₁,p₂,...,p_k,...p_M]为接触面的法向接触压力向量，C为影响矩阵，为法向位移向量；

S306、将矩阵式看作为二次规划问题，进而将振动器平板-大地接触问题演变为在满足相应边界条件情况下，对线性方程组的求解问题，线性方程组为：

p^*＝C^-1u；

满足的边界条件为：

接触区域内：

接触区域外：

其中，h(x,y)＝|z₁(x,y)-z₂(x,y)|为平板与大地接触前的初始间距，δ为平板在载荷作用下靠近大地的刚体位移，g(x,y)为变形发生后接触表面间的间隙。

优选的，在步骤S104中，利用Gauss-Seidel迭代法对线性方程组进行求解，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移，具体包括如下：

S401、导入离散化的平板与大地接触表面形貌数据，接触面材料属性以及法向载荷W，计算影响矩阵C，下三角矩阵L_matrix，对角阵D_matrix，上三角矩阵并满足关系式：

S402、移动平板，在平板与大地产生初始接触后，根据几何干涉计算干涉区域内的预定位移矩阵，利用Gauss-Seidel迭代法求解接触压力，为满足边界条件，当某一接触点的接触压力为负值时，强制将其设为零，此过程表达为：

S403、计算接触区域内产生接触压力的合力∑p，当接触合力与法向载荷满足关系式W＝∑p时，则迭代终止，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移。

优选的，在步骤S105中，基于最小二乘理论求解平板与大地之间接触力的关系式，具体包括如下：

S501、计算振动器平板与大地接触力-位移变化曲线；

S502、利用最小二乘法对曲线进行拟合，获得接触力的表达式：

接触力的表达式为：

式中，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，z为位移。

优选的，在步骤S106中，利用薄板振动理论引入大地表面的耦合效应，采用Pasternak地基模型模拟深层大地的反应，建立考虑地表形貌的可控震源振动器-大地耦合振动动力学模型，具体包括如下：

S601、在已建立的振动器平板-大地接触振动模型的基础上，引入重锤与平板之间的关系，建立完整的振动器振动模型，并将步骤S105中求得的接触力f_n的表达式代入振动器方程得到：

式中，m_r为重锤质量，m_b为平板质量，z₁为重锤位移，z₂为平板位移，z₂₀为静载作用下平板的静态位移，k₁为液压刚度，c₁为液压阻尼，c₂为接触阻尼，F_t为静态载荷，F(t)为扫描信号，该信号的频率在一个频段内随时间连续变化，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，为大地表面的位移，/>为静态载荷作用下大地表面的位移；

S602、振动过程中，振动器会捕获一部分大地与其发生耦合振动，利用薄板振动理论引入大地表面的耦合效应，其振动方程为：

式中，为大地表面的弯曲刚度，h为捕获大地深度，q(x,y,t)为振动器与大地之间的相互作用力，q(x,y,t)为地基反力；

S603、利用Pasternak模型描述深层大地的反应，其表达式为：

式中，k_1s为地基刚度，c_s为地基阻尼，G_s为切变系数；

S604、将振动器与大地之间的相互作用力看作为均布载荷q(x,y,t)，其表达式为：

式中，ξ表示Heaviside函数，a₁为二分之一平板宽度，b₁为二分之一平板长度，q₀表达式为：

S605、振动器-大地耦合振动动力学模型为：

本发明的有益效果是：本发明在建立可控震源振动器-大地耦合振动动力学模型时同时考虑了大地表面形貌引发的平板-大地接触不均以及大地表面的耦合效应等非线性问题，解决了现有可控震源振动器-大地系统模型难以反映振动器与大地之间的真实相互作用关系的问题，提高了模型的适用性和精确性，对准确预测计算振动器的实际输出力具有积极的现实意义。

附图说明

图1为本发明方法步骤流程示意图；

图2为振动器平板-大地相互作用模型示意图；

图3为振动器平板-大地接触示意图；

图4为接触算法迭代流程图；

图5为不同大地表面形貌下平板与大地之间接触刚度随载荷的变化曲线，图左为不同D下接触刚度曲线，图右为不同G下接触刚度曲线；

图6为不同大地表面形貌下平板-大地接触力随位移的变化曲线，图左为不同D下接触力曲线，图右为不同G下接触力曲线；

图7为振动器-大地耦合振动动力学模型图，图左为振动器结构示意图，图右为振动器-大地耦合振动动力学模型图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明的具体实施例之前，为使本发明的方案更加清楚完整，首先对本发明中出现的缩略语和关键术语定义进行说明：

表面分形模型：粗糙界面表面形貌的一种描述方法具有自相似和尺度独立的特性，能够从本质上描述粗糙表面的性质。

变分原理：变分原理是自然界静止(相对稳定状态)事物中的一个普遍适应的数学定律，也称最小作用定理，也是物理学的一条基本原理，以变分法来表达；针对发明中的变分原理，把一个力学问题用变分法化为求泛函极值(或驻值)的问题，就称为力学问题的变分原理。

请参阅图1-图7，本发明提供一种技术方案：一种考虑地表耦合的可控震源振动器输出力的计算方法，如图1所示，包括如下步骤：

步骤一：构建振动器平板-大地接触振动模型，在保证模型计算精度的前提下，忽略振动器平板的变形，将振动器平板与大地之间的接触等效成一个刚性平板和一个弹性粗糙表面的接触；

步骤二：采用表面分形模型构建大地表面形貌；

步骤三：基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组p^*＝C^-1u的求解问题；

步骤四：利用Gauss-Seidel迭代法对线性方程组p^*＝C^-1u进行求解，根据分析需要，输出平板与大地之间的真实接触面积，接触压力分布，接触面位移等；

步骤五：基于最小二乘理论求解平板与大地之间接触力的关系式；

步骤六：建立同时考虑地表形貌和大地表面耦合效应的可控震源振动器-大地耦合振动动力学模型；

步骤七：基于伽辽金法对可控震源振动器-大地耦合振动动力学模型进行离散，采用四阶的Runge-Kutta法对可控震源振动器-大地耦合振动动力学模型进行求解。

在本发明实施例中，步骤一中，振动器平板-大地接触振动模型如图2所示，可控震源振动器平板与大地之间的相互作用问题归根结底属于接触问题，正确描述平板与大地之间的接触作用关系，找到影响其作用的本质因素，是建立可控震源振动器平板-大地接触模型的基础。前期的研究表明，平板变形对振动器输出信号的影响较为有限，故忽略平板在振动过程中的变形，振动器平板-大地接触振动模型可以利用单自由度的质量-弹簧-阻尼模型进行描述；

在本发明实施例中，步骤二中，振动器平板与大地之间的接触面实际为具有一定表面形貌的三维粗糙表面，当忽略平板的变形时，平板与大地的接触可以等效成一个刚性平板和一个弹性粗糙表面的接触，其实际接触的区域是粗糙表面上的微凸体，等效后的振动器平板-大地接触示意图如图3所示。

在本发明实施例中，步骤二中，常用的粗糙表面的表面形貌描述方法包括统计学参数模型和表面分形模型，统计学参数模型依赖于测量仪器的精度和样本的尺寸等，具有尺度依赖性。相反，表面分形模型具有自相似和尺度独立的特性，能够从本质上描述粗糙表面的性质，研究表明，大地表面具有自相似和尺度独立性，因此可以采用表面分形模型进行描述。

表面分形模型具体的描述方法为采用修正的两参数Weierstrass-Mandelbrot函数生成弹性粗糙表面，修正的两参数Weierstrass-Mandelbrot函数的表达式为：

其中，L为样本长度，D为分形维数，G为分形粗糙度，M为生成分形表面的脊线数量，n为频率因子，n_max＝int[log(L/L_S)/logγ]_max，L_s为截断长度，Φ_m,n∈[0,2π]为随机相位，γ(γ>1)为缩放参数；大地表面的粗糙度由D和G决定；表面分形维数D表示粗糙表面所占据的空间程度大小，D值越大对应于越密集的表面形态；分形粗糙度G是高度尺度参数，G值越大则表面形貌越粗糙。

本发明实施例中，为考虑振动器平板的实际形状与尺寸大小，获得振动器平板-大地接触过程中接触力-变形的关系，从而推导出与位移相关的接触刚度公式，考虑采用数值求解法对振动器平板-大地接触模型进行求解。数值求解法通常采用先进的数字化技术获得表面形貌的数字化数据，或者基于相关数值方法采用软件编程的形式来生成随机粗糙表面。随后，在构造好表面形貌的基础上，利用计算机技术模拟粗糙表面的接触过程。

随着计算机性能的提升和粗糙表面测量技术的进步，数值接触模拟技术被广泛应用到粗糙表面接触过程的研究中。例如，随着有限元分析方法的发展，研究人员通过建立粗糙表面接触的有限元仿真模型来分析其接触特性。但是，对于平板-大地接触问题而言，由于接触面的名义接触面积较大，为充分还原大地表面形貌需要产生大量的网格，从而大大的增加工作量和计算时间，增加了有限元求解的难度。为克服上述困难，本发明采用一种基于变分原理结合数学规划的方法，该方法有效地提高了计算效率，同时保证了计算结果的唯一性，适用于振动器平板-大地接触作用问题。步骤三包括以下子步骤：

S31、根据弹性材料的应变余能在数值上等于弹性应变能U_E，将弹性材料的应变余能/>代入两个粗糙表面接触的总应变余能表达式里面获得总应变余能V^*。

S32、利用布森内斯解将表面位移和接触压力p联系起来，从而获得法向载荷作用下的接触压力p与表面位移/>之间的关系式。

S33、将表面位移代入接触压力p与表面位移/>之间的关系式，从而将求总应变余能最小值的接触问题归结为求解相应积分式关于接触压力p(x′,y′)的最小值问题。

S34、将积分式在整个接触区域内进行离散，同时采用分段插值的方法获得每个离散单元的接触压力，从而将积分式转化为求和式。

S35、将接触压力向量p与表面位移之间的关系式离散化处理后，求得总应变余能V^*并将其化为矩阵形式。

S36、将矩阵式看作为二次规划问题，进而将振动器平板-大地接触问题演变为在满足相应边界条件式的情况下，对线性方程组的求解问题。

进一步的，步骤S31中，弹性材料的应变余能和总的应变余能V^*的表达式分别为：

其中，为两个弹性体的应变余能，p为接触压力，/>和/>为两个接触体的预定位移，/>为基于几何干涉假定的接触区域内两接触体的总预定位移，/>为待计算接触区域内的复合表面位移，等于两个接触体总的接触变形；

步骤S31中将弹性材料的应变余能代入两个粗糙表面接触的总应变余能表达式中得：

进一步的，步骤S32中接触压力p与表面位移之间的关系式为：

其中，E^*为复合杨氏模量，E₁,E₂分别为土壤和平板的杨氏模量，ν₁，ν₂分别为土壤和平板的泊松比。

进一步的，步骤S33中相应积分式为：

其中，为2a×2b矩形单元中心的表面位移。

进一步的，步骤S34中将积分式离散处理后获得的求和式为：

其中，为2a×2b矩形单元中心的表面位移，Ω_k为单元接触面积子域，c_num为根据平板与大地之间的刚体位移求得的几何干涉区域内的初始接触点，索引k和索引l分别代表接触压力产生的位置和表面位移发生的位置，而C_kl为影响矩阵，其表达式为：

进一步的，步骤S35中，将接触压力向量p与表面位移之间的关系式离散化处理后，求得总应变余能V^*为：

步骤S35中总应变余能的矩阵形式为：

其中,p^T＝[p₁,p₂,...,p_k,...p_M]为接触面的法向接触压力向量，C为影响矩阵，为法向位移向量。

进一步的，步骤S36中线性方程组为：

p^*＝C^-1u；

步骤S36中接触压力需满足的边界条件为：

接触区域内：

接触区域外：

其中，h(x,y)＝|z₁(x,y)-z₂(x,y)|为平板与大地接触前的初始间距，δ为平板在载荷作用下靠近大地的刚体位移，g(x,y)为变形发生后接触表面间的间隙。

在本发明实施例中，步骤四中，关于线性方程组的数值解法一般分为两类：直接法和迭代法。其中，直接法是经过有限步算术运算，从而求得线性方程组精确解的方法，但实际运算过程中由于舍入误差的存在和影响，该方法也只能求得线性方程组的近似解。直接法一般适用于求解低阶稠密矩阵方程组。但是，对于振动器平板与大地之间的接触问题，计算过程中容易产生零元素较多的大型稀疏矩阵方程组，从而引起计算消耗巨大和矩阵存储困难等问题。相反，迭代法是采用某种极限方式去逐步逼近线性方程组的精确解，具有对计算机存储要求较低，程序设计简便，原始稀疏矩阵在计算过程中保持不变等优点，更适用于求解本发明的问题。目前，已有多种不同的线性方程组的迭代求解格式，包括Gauss-Seidel迭代法、Jacobi迭代法、Successive Over Relaxation迭代法、Conjugate Gradient法等，这些方法各有其优缺点，但是求解精度完全能满足一般的精度要求。本发明采用Gauss-Seidel迭代法，因为其具有比简单迭代法收敛快的优点，同时由于接触问题中的影响矩阵为对称正定，Gauss-Seidel迭代法必收敛。

步骤四中，Gauss-Seidel迭代法的基本求解原理如下：

设有线性方程组：

简记为Ax＝b。

选取分裂矩阵M为系数矩阵A的下三角部分，即选取M＝D-L(下三角矩阵)，A＝M-N，于是可得到解Ax＝b的Gauss-Seidel迭代法：

其中B＝I-(D-L)^-1A＝(D-L)^-1U≡G，f＝(D-L)^-1b。称G＝(D-L)^-1U为解Ax＝b的Gauss-Seidel迭代法的迭代矩阵。

记Gauss-Seidel迭代法的分量计算公式为：

由上式可得：

(D-L)x^(k+1)＝Ux^(k)+b

化为：

Dx^(k+1)＝Lx^(k+1)+Ux^(k)+b

因此求解线性方程组Ax＝b的Gauss-Seidel迭代法计算公式为：

在本发明实施例中，步骤四中，Gauss-Seidel迭代法的实施过程，具体包含以下内容：

S41、导入离散化的平板与大地接触表面形貌数据，接触面材料属性，以及法向载荷W，计算影响矩阵C，下三角矩阵L_matrix，对角阵D_matrix，上三角矩阵并满足关系式：

S42、移动平板，在平板与大地产生初始接触后，根据几何干涉计算干涉区域内的预定位移矩阵，利用Gauss-Seidel迭代法求解接触压力，为满足边界条件，当某一接触点的接触压力为负值，强制将其设为零，此过程可表达为：

S43、计算接触区域内产生接触压力的合力Σp，当接触合力与法向载荷满足关系式W＝Σp时，则迭代终止，根据分析需要，输出平板与大地之间的真实接触面积，接触压力分布，接触面位移等。

步骤四中，Gauss-Seidel迭代法的迭代算法如图4所示。

在本发明实施例中，步骤五中，基于最小二乘理论求解平板与大地之间接触力和接触阻尼的关系式，具体包含以下内容：

S51、利用步骤四计算振动器平板与大地的接触力-位移变化曲线。

S52、利用最小二乘法对上述曲线进行拟合，获得接触力的表达式。

所述步骤S52中，接触力的表达式为：

式中，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，z为位移。

在本发明实施例中，步骤S51中，根据计算结果绘制振动器平板与大地接触力-位移变化曲线。

为了确定不同大地表面形貌下平板与大地之间的接触刚度特性，建立具有不同大地表面形貌的平板-大地接触模型，求解载荷作用下接触表面的位移，绘制接触力-位移曲线，曲线的斜率即为平板与大地之间的接触刚度。图5为不同大地表面形貌下平板与大地之间接触刚度随载荷的变化曲线。

在本发明实施例中，如图6所示，为不同大地表面形貌下平板-大地接触力随位移的变化曲线，平板与大地之间的接触力-位移曲线存在非线性特性，曲线的斜率随接触力的增加而减小。

在本发明实施例中，步骤S52中，接触力的表达式中待定系数数值，如表1所示。

表1

在本发明实施例中，步骤六中，由可控震源振动器工作原理可知，液压系统输出的液压力作用在重锤上，使其产生上下振动，同时其产生的反作用力作用在活塞杆上，通过活塞杆、上板以及立柱传递给平板，并最终传递给大地。因此，振动器可以采用二自由度的质量-弹簧-阻尼模型进行描述，其中：重锤为一个自由度，而活塞杆、上板、立柱以及平板共同组成第二个自由度。重锤与平板之间的作用关系与进入重锤腔体的液压油有关。由于液压油刚度和液压油阻尼的计算过程繁琐复杂，以往的研究大都忽略液压系统的非线性，同时，本发明的主要关注点为振动器与大地之间的相互作用，因此在模型建立中采用线性弹簧和线性阻尼来描述重锤与平板之间的关系为考虑大地表面的耦合效应，将大地表面看作为一个放置在Pasternak地基上的薄板，通过平板与大地之间的相互作用力，将二自由度的振动器模型与大地模型进行耦合，最终建立的振动器-大地耦合振动动力学模型如图7所示。

将步骤五中求得的接触力f_n的表达式代入振动器方程得：

式中，m_r为重锤质量，m_b为平板质量，z₁为重锤位移，z₂为平板位移，z₂₀为静载作用下平板的静态位移，k₁为液压刚度，c₁为液压阻尼，c₂为接触阻尼，F_t为静态载荷，F(t)为扫描信号，该信号的频率在一个频段内随时间连续变化，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，为大地表面的位移，/>为静态载荷作用下大地表面的位移。

振动过程中，振动器会捕获一部分大地与其发生耦合振动，利用薄板振动理论引入大地表面的耦合效应，其振动方程为：

/>

式中，为大地表面的弯曲刚度，h为捕获大地深度，q(x,y,t)为振动器与大地之间的相互作用力，q(x,y,t)为地基反力；

利用Pasternak模型描述深层大地的反应，其表达式为：

式中，k_1s为地基刚度，c_s为地基阻尼，G_s为切变系数；

将振动器与大地之间的相互作用力看作为均布载荷q(x,y,t)，其表达式为：

式中，ξ表示Heaviside函数，a₁为二分之一平板宽度，b₁为二分之一平板长度，q₀表达式为：

引入大地的耦合效应，振动器-大地耦合振动动力学模型为：

/>

在本发明的实施例中，步骤七中，针对振动器-大地耦合振动动力学模型，采用伽辽金法对其进行离散，得到M×N+2个常微分方程组为：

其中：i＝1,2,3,…,M；j＝1,2,3,,N。

在本发明的实施例中，步骤七中，考虑可控震源振动器-大地耦合振动动力学模型属于非线性振动模型，其振动方程为非线性微分方程，除了极少数的特殊情况，难以采用解析法求得其精确解。因此，本发明采用数值方法对系统的振动响应进行求解。非线性振动的数值求解基础为常微分方程组的数值解法，其包括：Runge-Kutta法，线性加速度法，Euler法，Wilson-θ法等，考虑四阶Runge-Kutta法可达四阶精度O(h5)，数值稳定性较好，且只需知道一阶导数，无需明确定义或计算其他高阶导数就可进行迭代，因此本发明选用四阶Runge-Kutta法。

四阶Runge-Kutta法的具体实施方法为：

考虑一阶常微分方程的初值问题为：

y′＝f(x,y)

y(x₀)＝y₀

用数值方法求函数y，其基本思想是从y(x₀)点出发计算当

x＝x₀+ΔH,x₀+2ΔH,...x₀+nΔH

其中ΔH＝x_n+1-x_n为步长。

四阶Runge-Kutta法的计算公式如下：

针对振动器-大地耦合振动动力学模型，其振动方程由两个二阶微分方程组成，为采用上述的四阶Runge-Kutta法，需先化为2×(M×N+2)个一阶微分方程，其满足的初始条件为：

由于振动方程中存在非线性项，为了平衡计算时间并保证结果精度，步长ΔH取为1×10^-4s，对振动方程进行迭代求解，求得振动器-大地耦合系统的振动响应，并最终获得振动输出力。

本发明在建立可控震源振动器-大地耦合振动动力学模型时考虑了大地表面形貌引发的平板-大地接触不均以及大地表面的耦合效应等非线性问题，解决了现有可控震源振动器-大地系统模型无法真实反映振动器与大地之间的相互作用关系的问题，提高了模型的适用性和精确性，对准确预测计算振动器的实际输出力具有积极的现实意义。

尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于，包括如下步骤：

S101、构建振动器平板-大地接触振动模型，将振动器平板与大地的接触等效成一个刚性平板和一个弹性粗糙表面的接触；

S102、采用表面分形模型构建大地表面形貌；

S103、基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组的求解问题；

S104、利用Gauss-Seidel迭代法对线性方程组进行求解，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移；

S105、基于最小二乘理论求解平板与大地之间接触力的关系式；

S106、利用薄板振动理论引入大地表面的耦合效应，采用Pasternak地基模型模拟深层大地的反应，建立考虑地表形貌的可控震源振动器-大地耦合振动动力学模型；

S107、基于伽辽金法对可控震源振动器-大地耦合振动动力学模型进行离散；

S108、利用四阶Runge-Kutta法对可控震源振动器-大地耦合振动动力学模型进行求解，获得振动输出力。

2.根据权利要求1所述的考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于：在步骤S102中，表面分形模型具有自相似和尺度独立的特性，能够从本质上描述粗糙表面的性质；弹性粗糙表面的大地表面形貌采用修正的两参数Weierstrass-Mandelbrot函数生成；所述修正的两参数Weierstrass-Mandelbrot函数的表达式为：

其中，L为样本长度，D为分形维数，G为分形粗糙度，M为生成分形表面的脊线数量，n为频率因子，n_max＝int[log(L/L_S)/logγ]_max，L_s为截断长度，Φ_m,n∈[0,2π]为随机相位，γ为缩放参数，γ>1；大地表面的粗糙度由D和G决定；表面分形维数D表示粗糙表面所占据的空间程度大小，D值越大对应于越密集的表面形态；分形粗糙度G是高度尺度参数，G值越大则表面形貌越粗糙。

3.根据权利要求1所述的考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于：在步骤S103中，基于变分原理结合数学规划的方法将振动器平板与大地之间的接触问题转化为对线性方程组的求解问题，具体包括如下：

S301、设弹性材料的应变余能和总应变余能V^*的表达式分别为：

其中，为两个弹性体的应变余能，p为接触压力，/>和/>为两个接触体的预定位移，/>为基于几何干涉假定的接触区域内两接触体的总预定位移，/>为待计算接触区域内的复合表面位移，等于两个接触体总的接触变形；

根据弹性材料的应变余能在数值上等于弹性应变能U_E，将弹性材料的应变余能/>代入两个粗糙表面接触的总应变余能表达式里面，获得总应变余能V^*的表达式为：

S302、法向载荷作用下的接触压力p与表面位移之间的关系式为：

其中，E^*为复合杨氏模量，E₁,E₂分别为土壤和平板的杨氏模量，ν₁，ν₂分别为土壤和平板的泊松比；

S303、将表面位移代入接触压力p与表面位移/>之间的关系式，将求总应变余能最小值的接触问题转化为求解相应积分式关于接触压力p(x′,y′)的最小值问题，相应积分式为：

其中，为2a×2b矩形单元中心的表面位移，E^*为复合杨氏模量；

S304、将积分式在整个接触区域内进行离散，同时采用分段插值的方法获得每个离散单元的接触压力，从而将积分式转化为求和式；将积分式离散处理后获得的求和式为：

其中，为2a×2b矩形单元中心的表面位移，Ω_k为单元接触面积子域，c_num为根据平板与大地之间的刚体位移求得的几何干涉区域内的初始接触点，索引k和索引l分别代表接触压力产生的位置和表面位移发生的位置，而C_kl为影响矩阵，其表达式为：

S305、将接触压力向量p与表面位移之间的关系式离散化处理后，求得总应变余能V^*为：

将总应变余能V^*化为矩阵式：

其中,p^T＝[p₁,p₂,...,p_k,...p_M]为接触面的法向接触压力向量，C为影响矩阵，为法向位移向量；

S306、将矩阵式看作为二次规划问题，进而将振动器平板-大地接触问题演变为在满足相应边界条件情况下，对线性方程组的求解问题，线性方程组为：

p^*＝C^-1u；

满足的边界条件为：

接触区域内：

接触区域外：

其中，h(x,y)＝|z₁(x,y)-z₂(x,y)|为平板与大地接触前的初始间距，δ为平板在载荷作用下靠近大地的刚体位移，g(x,y)为变形发生后接触表面间的间隙。

4.根据权利要求1所述的考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于：在步骤S104中，利用Gauss-Seidel迭代法对线性方程组进行求解，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移，具体包括如下：

S401、导入离散化的平板与大地接触表面形貌数据，接触面材料属性以及法向载荷W，计算影响矩阵C，下三角矩阵L_matrix，对角阵D_matrix，上三角矩阵并满足关系式：

S402、移动平板，在平板与大地产生初始接触后，根据几何干涉计算干涉区域内的预定位移矩阵，利用Gauss-Seidel迭代法求解接触压力，为满足边界条件，当某一接触点的接触压力为负值时，强制将其设为零，此过程表达为：

S403、计算接触区域内产生接触压力的合力∑p，当接触合力与法向载荷满足关系式W＝∑p时，则迭代终止，输出平板与大地之间的真实接触面积、接触压力分布以及接触面位移。

5.根据权利要求1所述的考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于：在步骤S105中，基于最小二乘理论求解平板与大地之间接触力的关系式，具体包括如下：

S501、计算振动器平板与大地接触力-位移变化曲线；

S502、利用最小二乘法对曲线进行拟合，获得接触力的表达式：

接触力的表达式为：f_n＝k₂zⁿ¹；

式中，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，z为位移。

6.根据权利要求1所述的考虑地表耦合的可控震源振动器输出力的计算方法，其特征在于：在步骤S106中，利用薄板振动理论引入大地表面的耦合效应，采用Pasternak地基模型模拟深层大地的反应，建立考虑地表形貌的可控震源振动器-大地耦合振动动力学模型，具体包括如下：

S601、在已建立的振动器平板-大地接触振动模型的基础上，引入重锤与平板之间的关系，建立完整的振动器振动模型，并将步骤S105中求得的接触力f_n的表达式代入振动器方程得到：

式中，m_r为重锤质量，m_b为平板质量，z₁为重锤位移，z₂为平板位移，z₂₀为静载作用下平板的静态位移，k₁为液压刚度，c₁为液压阻尼，c₂为接触阻尼，F_t为静态载荷，F(t)为扫描信号，该信号的频率在一个频段内随时间连续变化，待定系数k₂，n₁与接触面的形貌以及接触面的材料属性有关，为大地表面的位移，/>为静态载荷作用下大地表面的位移；

S602、振动过程中，振动器会捕获一部分大地与其发生耦合振动，利用薄板振动理论引入大地表面的耦合效应，其振动方程为：

式中，为大地表面的弯曲刚度，h为捕获大地深度，q(x,y,t)为振动器与大地之间的相互作用力，q(x,y,t)为地基反力；

S603、利用Pasternak模型描述深层大地的反应，其表达式为：

式中，k_1s为地基刚度，c_s为地基阻尼，G_s为切变系数；

S604、将振动器与大地之间的相互作用力看作为均布载荷q(x,y,t)，其表达式为：

式中，ξ表示Heaviside函数，a₁为二分之一平板宽度，b₁为二分之一平板长度，q₀表达式为：

S605、振动器-大地耦合振动动力学模型为：