本发明公开了一种基于奇异边界法的波动类型动态数据重构方法，直接应用波动方程时间基本解作为核函数，用源点强度因子代替源点奇异项，无数值积分、奇异积分，无网格划分，直接应用波动方程时间基本解作为核函数，无复杂的数学变换，其简单高效的特点切合波动类型动态数据重构需要数据重构技术快速、稳定、精确的特征要求；实验对比表明，应用本发明所提出的技术，处理波动类型动态数据重构问题，在取得相似精度条件下，一般耗时只需传统线性边界元算法的10%左右，具有精度高，计算快，数学简单，程序简便的特点，可应用于声波减噪，海浪消能，地震预测等波动型数据重构，图像处理等工程技术领域。
1.一种基于奇异边界法的波动类型动态数据重构方法，其特征在于：包括以下步骤：(1)在所考察物质的内部和边界配置若干测试点，获得这些测试点的波动数据值；
(2)直接采用三维波方程的时间基本解作为核函数，建立波传播问题相应的插值矩阵；
(3)利用经验公式计算处理波动数据重构的源点强度因子；
(4)将源点强度因子代入插值矩阵，计算插值矩阵的未知系数；
其中，Ω表示所考察物质的区域，(x，y，z)为空间坐标，t为(x，y，z)所对应的时刻，c代表波速，u代表势函数，代表边界条件，u0和υ1表示初值条件，Δ表示Laplace算子，在声波传播中表示声压值；
其中δ表示狄拉克函数，c表示波速，t和τ表示配点和源点的时刻，r表示距离；由于波的传播需要时间，基本解G仅仅用延迟时刻的相应声压值的相关组合便可表示当前时刻的声压值，延迟时刻取决于源点与待测点之间的距离和波的传播速度；因此，对于测试点，只要求得关于源点{sj}的延迟时刻的未知系数{αj}，则测试点所求时刻的声压值即可求出；
原方程u被看作一个初边值问题，应用叠加原理，将其拆分为u1和u2，如式(3)和式(4)所示，即 其中u1被看作一个边值问题，u2被看作一个初值问题：式(4) 用三维泊松方程直接求出：
其中 表示半径为ct以M为圆心的球面；在式(4)中，仅需要求出测试点M影响区域中的声压值，所以在影响域 中布置初始数据采集点；
其后，考虑方程(3) 由惠更斯原理“波传播中所到达的每一点都可以作为一个新的次波源，所有这些次波所形成的包络面构成下一时刻的新波面”，假定在边界存在一系列随时间变化的点波源 这些点波源所发出的次波的总和形成计算域Ω中其后任意时刻的声压值 其中 表示声压强度，xm表示边界点，则计算域Ω中的声压被表示为
在边界Γ上布置数据采集点，以Δt为时间间隔，v表示空间积分变量，tn∈((n-1)Δt，nΔt)，t表示计算时刻，t的下标n表示时间层数，声压 被基本解G的一系列线性组合来近似：其中，αj(tR)表示对于源点sj的在延迟时刻tR的未知系数，Gij表示在延迟时刻tR的三维波方程基本解，Qi表示在延迟时刻tR的源点强度因子，xi表示第i个配点，xj表示第j个源点，N为数据采集点总数；
其中，0≤tmΔt-tR＜Δt，αj(tR)表示延迟时刻未知系数， 表示分离时间变量t后的延迟时刻未知系数，m和n表示时间层数，如果m＜n，则未知系数αj(tR)已经被前步求出，如果tR＜0，则表示波未传至配点处， ω表示波的频率，k＝ω/c表示波数；
式(5)、式(7)和式(8)共同构成了奇异边界法处理波传播问题的差值矩阵。
2.根据权利要求1所述的基于奇异边界法的波动类型动态数据重构方法，其特征在于：步骤(1)中，初始时刻在所考察物质的内部配置N1个测试点，获得初始时刻这些测试点的波动数据值ui，i＝1，...，N1；在所考察物质的表面配置N2个测试点，随着波动的进行，获得不同时刻这些点的波动数据值ui，i＝N1+1，...，N1+N2。
3.根据权利要求1所述的基于奇异边界法的波动类型动态数据重构方法，其特征在于：步骤(3)中，计算源点强度因子的经验公式为：
其中Aj表示源点sj的影响区域，S表示整个计算区域的表面积。
4.根据权利要求3所述的基于奇异边界法的波动类型动态数据重构方法，其特征在于：根据步骤(3)所推得源点强度因子，未知系数由式(7)和(8)求解得出：通过式(7)和(8)，求出在第n时间层，nΔt秒的未知系数 进而由式(8)求得延迟时刻的未知系数αj(tR)。
5.根据权利要求4所述的基于奇异边界法的波动类型动态数据重构方法，其特征在于：步骤(5)在时刻t任意内点x的声压u通过无积分计算公式求得：
一种基于奇异边界法的波动类型动态数据重构方法
技术领域
[0001] 本发明涉及一种波动类型动态数据重构方法，具体涉及一种奇异边界法的波动类型动态数据重构方法。
背景技术
[0002] 波动现象作为自然界最常见的自然现象之一，在工程应用领域具有广泛应用与影响，如声波减噪、海浪消能、地震波预测等科研工程领域，均需处理大量波动类型动态数据。由于波动现象基本控制方程的独特性，如波传播的滞后性、三维无后效性等，使波动方程基本解具有与扩散方程、拉普拉斯方程基本解有明显的不同形式，这使得波动类型动态数据的重构极为困难。
[0003] 传统上处理波动类型动态数据重构，一般有有限元方法、边界元方法、有限差分法等。但作为传统区域型网格方法，有限元在处理如声波散射，声波辐射等一系列无限域问题的数据重构时，往往会出现网格划分困难、计算速度慢、计算荷载大等难以克服的问题。而边界元方法作为边界型网格方法，虽然在一定程度上克服了有限元方法需要划分区域网格的缺点，提高了计算效率，但由于在计算过程中需要处理大量奇异与近奇异积分，严重影响了数据的重构速度。基本解方法作为一种无网格方法一定程度上克服了传统网格类算法需要划分网格的缺陷，但由于在基本解方法中为克服奇异性引入了虚假边界，使得该方法在处理一些复杂工况条件下的数据重构问题时，常出现数值不稳定的现象(见文献1.Young D L,Gu M H,Fan C M.The time-marching method of fundamental solutions for wave equations[J].Engineering Analysis with Boundary Elements,2009,33(12):1411-1425.)。
发明内容
[0004] 发明目的：本发明的目的在于针对现有技术的波动类型动态数据重构技术过程复杂、计算荷载大计算时间长的缺陷，提供一种简单、高效、无积分、无网格的基于奇异边界法的波动类型动态数据重构方法，从而达到节约重构时间，提高重构效率的目的。
[0005] 技术方案：本发明提供了一种基于奇异边界法的波动类型动态数据重构方法，包括以下步骤：
[0006] (1)在所考察物质的内部和边界配置若干测试点，获得这些测试点的初始数据和边界数据；
[0007] (2)直接采用三维波方程的时间基本解作为核函数，建立波传播问题相应的插值矩阵；
[0008] (3)利用经验公式计算处理波动数据重构的源点强度因子；
[0009] (4)将源点强度因子代入插值矩阵，计算插值矩阵的未知系数；
[0010] (5)计算任意时刻任意内点的波动数据值。
[0011] 进一步，步骤(1)中，初始时刻在所考察物质的内部配置N1个测试点，获得初始时刻这些测试点的波动数据值ui,i＝1,...,N1；在所考察物质的表面配置N2个测试点，随着波动的进行，获得不同时刻这些点的波动数据值ui,i＝N1+1,...,N1+N2。
[0012] 进一步，步骤(2)中采用的波动问题控制方程为：
[0013]
[0014] 其中，Ω表示所考察物质的区域，(x,y,z)为空间坐标，t为(x,y,z)所对应的时刻，c代表波速，u代表势函数，代表边界条件，u0和υ1表示初值条件，Δ表示Laplace算子，在声波传播中表示声压值；
[0015] 三维波传播问题对应基本解为：
[0016]
[0017] 其中δ表示狄拉克函数，c表示波速，t和τ表示配点和源点的时刻，r表示距离；由于波的传播需要时间，基本解G仅仅用延迟时刻的相应声压值的相关组合便可表示当前时刻的声压值，延迟时刻取决于源点与待测点之间的距离和波的传播速度；因此，对于特定测试点，只要求得关于源点{sj}的延迟时刻的未知系数{αj}，则测试点所求时刻的声压值即可求出；
[0018] 原方程u可以被看作一个初边值问题，应用叠加原理，将其拆分为u1和u2，如式(3)和式(4)所示，即 其中u1可以被看作一个边值问题，u2可以被看作一个初值问题：
[0019]
[0020]
[0021] 式(4) 用三维泊松方程直接求出：
[0022]
[0023] 其中 表示半径为ct以M为圆心的球面；在式(4)中，仅需要求出测试点M影响区域 中的声压值，所以在影响域 中布置初始数据采集点；
[0024] 其后，考虑方程(3) 由惠更斯原理“波传播中所到达的每一点都可以作为一个新的次波源，所有这些次波所形成的包络面构成下一时刻的新波面”，假定在边界存在一系列随时间变化的点波源 这些点波源所发出的次波的总和形成计算域Ω中其后任意时刻的声压值 其中 表示声压强度，xm表示边界点，则计算域Ω中的声压 可以被表示为
[0025]
[0026] 在边界Γ上布置数据采集点，以Δt为时间间隔，v表示空间积分变量，tn∈((n-1)Δt,nΔt)，t表示计算时刻，t的下标n表示时间层数，声压 可以被基本解G的一系列线性组合来近似：
[0027]
[0028] 其中，αj(tR)表示对于源点sj的在延迟时刻tR的未知系数，Gij表示在延迟时刻tR的三维波方程基本解，Qi表示在延迟时刻tR的源点强度因子，xi表示第i个配点，xj表示第j个源点，N为数据采集点总数；
[0029] 在简谐波动中，假定
[0030]
[0031] 其中，0≤tmΔt-tR＜Δt，αj(tR)表示延迟时刻未知系数， 表示分离时间变量t后的延迟时刻未知系数，m和n表示时间层数，如果m＜n，则未知系数αj(tR)已经被前步求出，如果tR＜0，则表示波未传至配点处， ω表示波的频率，k＝ω/c表示波数；
[0032] 式(5)、式(7)和式(8)共同构成了奇异边界法处理波传播问题的差值矩阵。
[0033] 进一步，步骤(3)中，计算源点强度因子的经验公式为：
[0034]
[0035] 其中Aj表示源点sj的影响区域，S表示整个计算区域的表面积。
[0036] 进一步，根据步骤(3)所推得源点强度因子，未知系数由式(7)和(8)求解得出：通过式(7)和(8)，求出在第n时间层，nΔt秒的未知系数 进而由式(8)求得延迟时刻的未知系数αj(tR)。
[0037] 进一步，步骤(5)在时刻t任意内点x的声压u通过无积分计算公式求得：
[0038]
[0039]
[0040]
[0041] 有益效果：本发明仅需边界配置边界数据采集点，部分区域配置初始数据采集点，无数值积分、奇异积分，无网格划分，直接应用波动方程时间基本解作为核函数，无复杂的数学变换，其简单高效的特点切合波动类型动态数据重构需要数据重构技术快速、稳定、精确的特征要求；实验对比表明，应用本发明所提出的技术，处理波动类型动态数据重构问题，在取得相似精度条件下，一般耗时只需传统线性边界元算法的10％左右，具有精度高，计算快，数学简单，程序简便的特点，可应用于声波减噪，海浪消能，地震预测等波动型数据重构，图像处理等工程技术领域。
附图说明
[0042] 图1为本发明波动类型动态数据重构方法流程图；
[0043] 图2(a)为三维波传播问题边界数据采集点布置示意图，图2(b)为区域数据采集点布置示意图；
[0044] 图3为边界数据采集点sj影响区域示意图；
[0045] 图4为轮胎形区域数据采集点布置示意图；
[0046] 图5为奇异边界法数据重构结果随边界数据采集点点数增加变化收敛图；
[0047] 图6为奇异边界法数据重构结果随时间延续变化趋势图；
[0048] 图7为声辐射无量纲实部声压重构数值结果随波数变化图；
[0049] 图8为声辐射无量纲虚部声压重构数值结果随波数变化图。
具体实施方式
[0050] 下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
[0051] 实施例：如图1所示，一种基于无网格奇异边界法重构波动类型动态数据，具体步骤如下：
[0052] (1)初始时刻在所考察物质的内部配置N1个测试点，获得初始时刻这些测试点的波动数据值ui,i＝1,...,N1；
[0053] 在所考察物质的表面配置N2个测试点，随着波动的进行，获得不同时刻这些点的波动数据值ui,i＝N1+1,...,N1+N2。
[0054] (2)根据奇异边界法的基本思想，直接采用三维波方程的时间基本解作为核函数，建立波传播问题相应的奇异边界法插值矩阵：
[0055] 应用叠加原理，原数据u可分为初始数据u2和边界数据u1，即
[0056] 初始数据 可以用三维泊松方程直接求出
[0057]
[0058] 其中 表示半径为ct以M为圆心的球面。在式(S.1)中，奇异边界法仅需要求出测试点M影响区域 中的声压值，所以奇异边界法仅需要在影响域 中布置初始数据采集点，如图2所示。
[0059] 其后，考虑边界数据 奇异边界法在边界Γ上布置边界数据采集点，如图2所示，以Δt为时间间隔，t的下标表示时间层数。声压 tn∈((n-1)Δt,nΔt)，可以被基本解G的一系列线性组合来近似，如式(S.2)所示：
[0060]
[0061] 其中αj(tR)表示对于源点sj的在延迟时刻tR的未知系数，Gij表示在延迟时刻tR的三维波方程基本解，Qi表示在延迟时刻tR的源点强度因子。
[0062] 在简谐波动中，我们假定
[0063]
[0064] 其中0≤tmΔt-tR＜Δt，如果m＜n，则未知系数αj(tR)已经被前步求出，如果tR＜0，则表示波未传至配点处， ω表示波的频率，k＝ω/c表示波数。
[0065] 式(S.1)，式(S.2)和式(S.3)共同构成了奇异边界法处理波传播问题的插值矩阵。
[0066] (3)利用经验公式计算奇异边界法的源点强度因子：三维波动方程的源点强度因子计算公式如式(S.4)所示。
[0067]
[0068] 其中Aj表示源点sj的影响区域，如图3所示，S表示整个计算区域的表面积。
[0069] (4)将源点强度因子代入奇异边界插值矩阵，计算奇异边界法求解方程的未知系数。应用步骤(3)所推得源点强度因子，步骤(4)中奇异边界法的未知系数可由式(S.2)和(S.3)求解得出。通过式(S.2)和(S.3)，求出在第n时间层，nΔt秒的未知系数 进而由式(S.3)求得延迟时刻tR的未知系数αj(tR)。
[0070] (5)根据奇异边界法公式，计算任意时刻任意内点的波动数据值。步骤(5)中，通过步骤(4)求得的未知系数，在时刻t任意内点x的声压u可通过下面奇异边界法公式求得：
[0071]
[0072]
[0073]
[0074] 实施例1：考虑如图4所示的轮胎形区域波动数据重构，区域方程为[0075]
[0076] 其中R＝0.8，r＝0.2，边界数据采集点如图3所示。
[0077] 波动数据控制方程及精确解为：
[0078]
[0079] 本例中，我们应用时间步长Δt＝2×10-1，波速c＝10，对每个测试点布置区域数据采集点Nf＝1255，测试点被布置在半径为0.8，t＝1s的时间层上，图5给出了奇异边界法随边界采集数据点点数增加的数值精度收敛图，其中波数分别为(k1＝0.5，k2＝5，k3＝10)。可以发现奇异边界法以2阶收敛速度快速收敛，当k＝5时，收敛速度甚至达到了C＝5.0。
[0080] 其后，应用时间步长Δt＝2×10-1，波速c＝10，布置边界数据采集点Ns＝972，对每个测试点布置区域数据采集点Nf＝1255，图6给出了在不同波数情况下(k1＝0.5，k2＝5，k3＝10)奇异边界法数据重构最大绝对误差随时间变化情况，可以看到随时间延续，奇异边界法最大相对误差仍与精确解保持一致，并未随时间变化出现离散情况。
[0081] 实施例2：考虑在单位球声压辐射数据重构，其中半径a＝1，波速c＝v0，精确解可以被表示为：
[0082]
[0083] 其中z0＝ρ0c，表示介质阻抗，ρ0为介质密度，c为波速，ω＝kc为波频。图7和图8给出了无量纲声压实部解 和虚部解 在边界数据采集点Ns＝400，时间步长Δt＝0.5s条件下的数值结重构结果。可以发现奇异边界法重构结果和精确解精确拟合，相比之下当虚边界d＝0.5时，基本解方法重构结果和精确解拟合良好，但是当虚边界取为d＝0.9时，数据重构结果发生了明显离散。同时，需要指明的是，对于声辐射，声散射等外域问题的波动数据重构，本发明方法无需采集区域波动数据，大大提高了波动类型数据的重构效率。
[0084] 综上，本发明波动类型动态数据重构方法以奇异边界法为基础，直接应用波动方程时间基本解作为核函数，用源点强度因子代替源点奇异项，避免了数值积分和网格，仅需边界布点，无网格、无数值积分和数学变换，提高了数据重构效率。相较于现有技术，无需做复杂的快速傅里叶变换或繁琐的多级波形迭代求解，直接在时域上运用时间依赖基本解作为插值基函数，通过应用源点强度因子重构波动类型数据，其简单高效的特点切合波动类型动态数据重构需要计算工具快速、稳定、精确的特征要求，适用于声波，水波，地震波等标量波传播动态数据的重构问题。
[0085] 实验对比表明，应用本发明方法处理波动类型动态数据重构问题，在取得相似精度条件下，一般耗时只需传统线性边界元算法的10％左右，具有精度高，计算快，数学简单，程序简便的特点。该技术为波动型动态数据的重构提供了新的，简单高效的技术路线。
