WO2024119366A1 - 一种基于路径跟踪的声波直边衍射计算方法 - Google Patents

一种基于路径跟踪的声波直边衍射计算方法 Download PDF

Info

Publication number
WO2024119366A1
WO2024119366A1 PCT/CN2022/136882 CN2022136882W WO2024119366A1 WO 2024119366 A1 WO2024119366 A1 WO 2024119366A1 CN 2022136882 W CN2022136882 W CN 2022136882W WO 2024119366 A1 WO2024119366 A1 WO 2024119366A1
Authority
WO
WIPO (PCT)
Prior art keywords
path
diffraction
nodes
probability
edge
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/CN2022/136882
Other languages
English (en)
French (fr)
Inventor
曹春晓
任重
安自立
周昆
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University
Original Assignee
Zhejiang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University filed Critical Zhejiang University
Priority to PCT/CN2022/136882 priority Critical patent/WO2024119366A1/zh
Publication of WO2024119366A1 publication Critical patent/WO2024119366A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S11/00—Systems for determining distance or velocity not using reflection or reradiation
    • G01S11/14—Systems for determining distance or velocity not using reflection or reradiation using ultrasonic, sonic or infrasonic waves
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions

Definitions

  • the invention relates to the technical field of acoustic simulation, and in particular to a method for calculating straight-edge diffraction of sound waves based on path tracking.
  • Sound propagation simulation technology is widely used in fields such as architectural design and virtual reality.
  • the method based on geometric acoustics is more common in practical applications due to its high computational efficiency.
  • the geometric acoustics method simplifies the wave equation and assumes that sound waves propagate in a straight line. This simplification requires independent processing of phenomena such as diffraction and interference.
  • the lack of diffraction components means that when the listener enters the "shadow" area blocked by the object, the sound propagation simulation results will show obvious discontinuity, which is inconsistent with the daily acoustic experience and greatly reduces the realism of the simulation results.
  • Straight edge diffraction refers to diffraction near two infinite half-planes connected by a straight edge and forming a certain angle.
  • geometric diffraction theory a method to describe convex straight edge diffraction in geometric optics.
  • Keller's results are close to the exact solution in most directions, but this model fails near the "shadow boundary", that is, the discontinuity of the geometric optics solution.
  • the improved new theory is called the uniform geometric diffraction theory.
  • MCNAMARA D A PISTORIUS C W I, MALHERBE J A G.
  • Biot-Tolstoy-Medwin (BTM) model Another description model of straight edge diffraction can be obtained, called the Biot-Tolstoy-Medwin (BTM) model.
  • BTM Biot-Tolstoy-Medwin
  • SVENSSON U P FRED R I
  • VANDERKOOY J An analytic secondary source model of edge diffraction impulse responses [J]. Journal of the Acoustical Society of America, 1999, 106: 2331-2344.
  • UTD and BTM are the two most common geometric diffraction models in sound propagation simulation.
  • the UTD (Consistent Geometric Theory of Diffraction) model has been extensively studied because it appeared earlier. For plane and spherical incident waves, ideal half-planes or curved half-planes, various boundary conditions, and other situations, relevant research articles can be found. However, the Kirchhoff approximation is used in the construction of UTD, and the straight edge length is assumed to be infinite. Therefore, in practical applications, the results of UTD are bound to be different from the exact solution. In addition, the representation of the UTD solution is not in the time domain like the general impulse response, but in the frequency domain. This makes it difficult to coordinate with other parts of the geometric sound propagation simulator.
  • the purpose of the present invention is to provide a method for simulating acoustic diffraction phenomena based on path tracking in view of the deficiencies of the prior art.
  • a method for calculating straight-edge diffraction of acoustic waves based on path tracking comprises the following steps:
  • the paths are composed of a series of nodes and line segments connecting the nodes in the scene, the starting node of the path is the sound source, the ending node is the listener, and the remaining intermediate nodes are located on the convex straight edges of the geometric body in the scene; except for the nodes constituting the paths, the other parts of the paths do not intersect with the surfaces of the geometric bodies in the scene;
  • L 0 (x 0 ⁇ x) 1 is the sound source characteristic
  • M ( xi ⁇ xi+1 ) is the sound medium characteristic at xi ⁇ xi+1
  • ⁇ (xi -1 ⁇ xi ⁇ xi+1 ) is the radial function representation of the straight edge diffraction at xi ;
  • step (2) includes the following sub-steps:
  • step (2.4) the pseudo intersection Randomly select one of the three edges of the triangle where it is located; let the triangle be ⁇ abc, the selected edge be ab, and the vertex opposite to the edge be c, then use a straight line to connect c and And extend it to intersect ab at a point, denoted as x n+1 , called the intersection point;
  • step (2.5) If the edge ab in step (2.4) is not a convex straight edge, or the line segment xn ⁇ xn+1 intersects the scene geometry at a place other than the two endpoints, the intersection fails and go to step (2.8); otherwise, add the intersection point xn+1 as a node to the end of the path x0 ⁇ x1 ⁇ ... ⁇ xn to form a new path x0 ⁇ x1 ⁇ ... ⁇ xn+1 ;
  • step (2.7) If the path length does not meet the user's requirements, go to step (2.2);
  • step (2) paths are generated simultaneously from both ends of the sound source and the listener, the path starting from the sound source is called the forward subpath, and the path starting from the listener is called the backward subpath; thereafter, connections are randomly established between the nodes of the forward and backward subpaths to form a complete path from the sound source to the listener; assuming that the probability of generating the part of the forward subpath from the sound source to the connected point is p F , the probability of generating the part of the backward subpath from the listener to the connected point is p B , and the probability of establishing a connection is p L , then the probability of generating the complete path is p F ⁇ p B ⁇ p L .
  • step (5) the contribution of the path needs to be multiplied by a factor C s for suppressing extreme values:
  • Co is the product of the Co corresponding to each randomly generated segment in the path
  • the Co of each path segment is a parameter that measures the possibility of the path containing the path segment generating an extreme value.
  • step (3) the contribution calculation needs to be divided by the product of the elements in [F j,i ] corresponding to the positions of the nodes connected in the forward and reverse paths;
  • [F j,i ] is an array of M rows and N columns, where N is the maximum number of nodes in the forward and backward paths;
  • the nodes of the path are stored in an array with the same number of rows and columns, where each row corresponds to a separate path, and the elements in the first column of the array [F j,i ] are set to 1, and the first column of the node array is the starting point of the path;
  • the generation process of the node in the i-th column of the array includes the following sub-steps:
  • step (5.2) the pseudo intersection Randomly select one of the three edges of the triangle where it is located; let the triangle be ⁇ abc, the selected edge be ab, and the vertex opposite to the edge be c, then use a straight line to connect c and And extend it to intersect ab at a point, denoted as x i-1 , called the intersection point;
  • step (5.4) If the edge ab in step (5.3) is not a convex straight edge, or the line segment x i-2 ⁇ xi-1 intersects the scene geometry at a place other than the two end points, the intersection fails and go to step (5.7); otherwise, add the intersection point x i-1 as a node to the end of the path x 0 ⁇ x 1 ⁇ ... ⁇ xi-2 in the current row to form a new path x 0 ⁇ x 1 ⁇ ... ⁇ xi-1 ;
  • step (5.6) If the path length does not meet the user's requirements, go to step (5.1) and prepare to generate the next column of nodes;
  • the method in this paper has higher computational efficiency and does not require pre-processing calculations such as pre-solving wave equations and geometric simplification when dealing with large-scale scenes with complex occlusion relationships.
  • the diffraction model used in this method does not have approximate operations such as Kirchhoff approximation and infinite straight edge assumption, so the result is more accurate. For simpler geometric bodies, the result is completely consistent with the theoretical diffraction solution.
  • the diffraction algorithm based on path tracking is a natural extension of the traditional sound propagation algorithm based on path tracking, so the method of the present invention is also easier to add to the existing sound propagation simulator.
  • FIG1 is a schematic diagram of a rectangular coordinate system in the present invention.
  • FIG2 is an absolute value image of the straight edge diffraction interface term function of the present invention.
  • FIG3 is a schematic diagram of the sampling probability distribution of the straight edge diffraction interface item of the present invention.
  • FIG4 is a schematic diagram of the arrangement of a first-order diffraction test scene
  • FIG5 is a schematic diagram of calculation results of the present invention and the control method under the first-order diffraction test scene arrangement
  • FIG6 is a schematic diagram of the arrangement of a second-order diffraction test scene
  • FIG. 7 is a schematic diagram of calculation results of the present invention and the control method under the second-order diffraction test scene arrangement.
  • the present invention provides a method for calculating straight-edge diffraction of sound waves based on path tracking, comprising the following steps:
  • the paths are composed of a series of nodes and line segments connecting the nodes in the scene, the starting node of the path is the sound source, the ending node is the listener, and the remaining intermediate nodes are located on the convex straight edges of the geometric body in the scene; except for the nodes constituting the paths, the other parts of the paths do not intersect with the surfaces of the geometric bodies in the scene;
  • L 0 (x 0 ⁇ x 1 ) is the sound source characteristic
  • M ( xi ⁇ xi+1 ) is the sound medium characteristic at xi ⁇ xi+1
  • ⁇ (xi -1 ⁇ xi ⁇ xi +1 ) is the radial function representation of the straight edge diffraction at xi ;
  • the path tracking-based straight-edge diffraction calculation method for acoustic waves is specifically as follows:
  • Path tracing The technique of path tracing in sound propagation simulation is the path tracing algorithm used in sound propagation simulation and the application of diffraction representation in the algorithm.
  • t represents time.
  • L(x′ ⁇ x, t) is the energy response or impulse response excited by the sound source x′ at position x, that is, the objective function to be solved in this equation
  • L 0 (x′ ⁇ x, t) represents the direct component, that is, the response excited directly at x by the sound source without interacting with the scene geometry
  • V(x′′, x) is the visibility function, which is 1 when x′′ and x are visible to each other, and 0 otherwise
  • ⁇ (x′ ⁇ x′′ ⁇ x) is the interface term, which represents the outgoing sound intensity generated near a point x′′ on the interface in the direction x′′ ⁇ x when the incident sound direction is x′ ⁇ x′′
  • M t (x′′, x′, t) is the medium term, which represents the influence of the sound medium between x′′ and x′, including energy absorption and propagation delay.
  • the asterisk * is
  • a is the acoustic attenuation coefficient of the medium
  • is the length of the vector x′′-x′
  • ⁇ ( ⁇ ) is the unit impulse response function
  • c is the speed of sound in the medium.
  • T represents an interaction between a sound wave and the scene geometry. Since T is essentially an integral, its value can be calculated using the Monte Carlo integration method.
  • the present invention represents TnL0 as the mathematical expectation of the path contribution; a path is composed of a series of nodes and line segments connecting the nodes in the scene, with the starting node being the sound source, the end node being the listener, and the intermediate nodes being located on the surface of the scene geometry; a path containing n intermediate nodes (n+2 nodes) is called an n-order path, corresponding to the term TnL0 in the series.
  • the mathematical expectation of its contribution needs to be divided by its generation probability when calculating it. Considering the time delay, the contribution expectation can be calculated as follows:
  • t S is the total time delay of the path.
  • the expectation value is the average of the results of each path.
  • Interface term expression of straight edge diffraction The key to using path tracking to calculate diffraction is to construct the corresponding interface term.
  • the present invention will give the form of the interface term.
  • Figure 1 the coordinate system near the straight edge and the representation method and related symbols of the incident and outgoing path directions are given.
  • the present invention calls n in the figure the normal, b the binormal, and t the tangent (the t direction coincides with the straight edge direction, both positive and negative).
  • the three constitute a frame at a point on the straight edge.
  • Half of the dihedral angle at the straight edge is denoted as ⁇ . Since this article only considers convex straight edges, the value of ⁇ should be less than ⁇ /2.
  • the incident direction vector vi and the outgoing direction vector v o are represented by a series of angles:
  • the interface term is:
  • the function constructed above satisfies the non-negativity and normalization requirements of the probability density function.
  • FIG3 is a schematic diagram of the sampling probability distribution of the straight edge diffraction interface term of the present invention, showing the sampling probability density function under the corresponding direction. It can be seen that the shapes of the functions in the two figures are very similar.
  • Diffraction path generation In the path tracing algorithm, the common path generation method is to generate each node in the path in sequence. Generally, a ray with a random direction is emitted from the end node of the current path. If it intersects with the scene geometry, the intersection point is used as the successor node of the original end node. However, when simulating diffraction, the middle node of the path must be located on the edge of the geometry, and the general ray intersection cannot guarantee that the intersection point meets this requirement. Therefore, it needs to be modified.
  • the intersection point between the ray and the geometry is called a pseudo-intersection point.
  • the intersection point between the ray and the geometry is called a pseudo-intersection point.
  • the path generation probability after adding nodes is as follows:
  • the probability of generating a complete path is pF ⁇ pB ⁇ pL .
  • pF is the probability of generating the part of the forward subpath from the sound source to the connected part
  • pB is the probability of generating the part of the backward subpath from the listener to the connected part
  • pL is the probability of establishing a connection.
  • the present invention can also adopt multiple importance sampling, which is a method for improving sample quality and is usually used in conjunction with bidirectional path tracking.
  • multiple importance sampling requires giving the generation probability of this path under all generation methods.
  • a path x0 ⁇ x1 ⁇ ... ⁇ xn is generated through bidirectional path tracking. Its connecting segment may be any segment in the path, so there are a total of n-1 possible generation methods. The method for calculating the path generation probability under these generation methods.
  • corresponding forward and backward pseudo-intersections must be generated for each node that is not the starting point or end point, where the forward pseudo-intersection is visible to the predecessor node of the node, and the backward pseudo-intersection is visible to the successor node of the node; specifically as follows:
  • the generation probability of this path is p(x 0 ⁇ x 1 ⁇ ... ⁇ xi )p( xi+1 ⁇ xi +2 ⁇ ... ⁇ x n )p L ( xi ⁇ xi +1 ).
  • x i+1 ⁇ xi+2 ⁇ ... ⁇ x n means the same as x n ⁇ x n-1 ⁇ ... ⁇ xi+1 .
  • the scene positions involved in calculating the path generation probability are listed as follows:
  • each intermediate node needs to have corresponding forward and backward pseudo-intersections. Therefore, it is necessary to actively fill in the missing pseudo-intersections in the diagram.
  • the path generation probability under all generation methods can be calculated.
  • the expected contribution of the final path should be divided by the product of the two split degrees corresponding to the forward and reverse path connection nodes. In this way, the empty spaces in the array can be fully utilized, storage space can be saved, and the node connection efficiency can be improved.
  • the inventor implemented the algorithm described above on a computer equipped with an Intel i7 2.60GHz quad-core CPU and 12GB memory.
  • the specific technology selected was bidirectional path tracking combined with multiple importance sampling, and all optimization techniques in Section 5 of the implementation process description were used.
  • the inventor compared its results with the simulation results using the UTD and BTM models.
  • Figures 4 and 6 show the first-order and second-order diffraction scenarios for comparison.
  • the scenario in Figure 4 contains a straight dihedral with a straight side length of 2m, while the scenario in Figure 6 contains two parallel right dihedrals.
  • the path from the sound source must undergo two diffractions before it can reach the listener.
  • Both scenario geometries use Dirichlet boundary conditions.
  • Figures 5 and 7 show the impulse responses calculated using different methods in the scenario.
  • the methods used here include BTM, UTD, the method of this paper, and the results reconstructed by selecting 8 frequency points in UTD (located between 62.5-8000Hz, distributed in an exponential interval).
  • the error of the UTD method in calculating the second-order and above diffraction is very large, which is not shown in Figure 7.
  • the quality of the calculation results of the algorithm in this paper in the case of first to second order diffraction is not significantly different from that of the existing methods, and is significantly better than the common simplified UTD scheme (only calculating the results at 8 frequencies). Since the BTM method is accurate for the solution of the first-order diffraction, this comparison also verifies the correctness of the method in the present invention. For large-scale complex scenarios, this method shows near real-time computing efficiency. When 12,000 paths are used, for a model containing 22,974 triangular faces, the time overhead for calculating the diffraction impulse response of the present invention is about 140 ms. For a model containing 279,133 faces, the time overhead is about 205 ms. It can be seen that the method of the present invention is suitable for fast calculation of diffraction results in complex scenes.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于路径跟踪的声波直边衍射的计算方法。该方法将凸直边衍射表示为几何体直边上径向函数的积分，并使用路径跟踪方法计算该积分，从而实现了路径跟踪框架下的高效直边衍射计算。此外，本发明针对上述衍射计算方法提出了一系列配套的优化技术，以提高计算结果的质量，减少计算结果的资源消耗。实验结果表明，本发明的方法在简单场景下的结果与凸直边衍射的精确解相符，在大型复杂场景下，该方法可以接近实时的速度向用户提供衍射模拟结果。

Description

一种基于路径跟踪的声波直边衍射计算方法 技术领域
本发明涉及声学模拟技术领域，尤其涉及一种基于路径跟踪的声波直边衍射计算方法。
背景技术
声传播模拟技术在建筑设计和虚拟现实等领域中有广泛的应用。在声传播模拟算法中，基于几何声学的方法由于计算效率较高，在实际应用中较为常见。几何声学方法对波动方程进行简化近似，假设声波沿直线传播。这种简化使得衍射，干涉等现象需要独立处理。衍射成分的缺乏使得听者进入物体遮挡的“阴影”区域时，声传播模拟结果会出现明显的不连续现象，与日常声学体验不符，极大降低了模拟结果的真实感。
将衍射现象置入几何声学的框架是一件复杂的事情。实际工程一般会在使用普通几何声学理论计算出的冲激响应之上，另外叠加一个表示衍射的冲激响应。这里表示衍射的冲激响应可使用波动方程算出。但在更多的情况下，希望衍射的计算和普通几何方法相近。一种自然的想法是先场景将场景拆分成一些较简单，容易计算衍射的结构，之后将计算结果拼起来，构成最终的衍射解。一般使用三角网格表示场景中的物体，而物体遮挡产生的阴影边界和三角网格中的边有关，因此，可以尝试计算三角网格各条边附近的衍射。这就是几何声学关心直边衍射表示的动机。
直边衍射，或尖劈衍射，指由一条直边相连的，成一定夹角的两无限大半平面附近的衍射。上世纪中，Keller首先提出了一种在几何光学中表述凸直边衍射的方法，称为几何绕射理论。Keller的结果在大多数方向上接近精确解，但在“阴影边界”，即几何光学解的不连续处附近，这一模型是失效的。之后，另一些人对GTD做出改进，解决了其在阴影边界附近不连续的问题，改进后的新理论称作一致性几何绕射理论，具体可参考MCNAMARA D A,PISTORIUS C W I,MALHERBE J A G.Introduction to the Uniform Geometrical Theory of Diffraction[M].Artech House,1990。另一种对直边衍射的表达尝试始于Biot-Tolstoy解。其给出了球面波射入凸二面角时回波的解析形式。在将Biot-Tolstoy解应用于声学研究时，Medwin基于波动方程惠更斯原理的直觉，将回波分解为以直边上各点为中心的一系列径向函数的组合，并发现这一拆分有助于解释有限长直边的衍射问题。结合Biot-Tolstoy解和Medwin分解，便可得到另一种直边衍射的描述模型，称为Biot-Tolstoy-Medwin(BTM)模型，具体可参考SVENSSON U P,FRED R I,VANDERKOOY J.An analytic secondary source  model of edge diffraction impulse responses[J].Journal of the Acoustical Society of America,1999,106:2331-2344。
UTD和BTM是声传播模拟中最为常见的两种几何衍射模型。UTD(一致性几何绕射理论)模型由于出现较早，其相关研究非常丰富。对于平面和球面入射波，理想半平面或曲半平面，各种不同边界条件等情况，都可以找到相关的研究文章。然而，UTD的构建过程中使用了基尔霍夫近似，且假设直边长度无限。因此在实际应用中，UTD的结果与精确解必然有差异。此外，UTD解的表示并非如一般冲激响应一样位于时域上，而是位于频域上。这使得它和几何声传播模拟器的其它部分难以配合。通常的应用往往只会选取UTD解数个频率点上的结果，并丢掉全部相位信息。这使得UTD应用于实际声传播模拟器时的误差远大于理论误差。UTD模型所依赖的Biot-Tolstoy解是无限长直边衍射问题的解析解，Medwin分解也并非近似。因此，BTM方法理论上可达到的精确程度比UTD要高。此外，由于Medwin分解将直边分解为一系列径向函数，BTM方法还可以处理分段直边，甚至曲线边的衍射问题。BTM不及UTD的方面在于其相关研究资料较少，且在现有声传播模拟器中的实现不如UTD有效率。
UTD和BTM在声传播模拟中均有广泛的应用。本世纪初便有文章在声学虚拟现实系统中实现了基于UTD的单次衍射，后来的研究者很快将其推广至高次衍射的情况。对于较复杂场景，也存在快速寻找可行衍射路径的方法，但需要对场景进行预先简化，直接使用原始的BTM表达式进行声传播模拟，但其计算效率过低，仅适合小规模场景的离线计算。也有以BTM为基础，计算常见几何体的衍射解，并使用这些解去近似实际环境中的物体。如体积衍射及传播法。
除UTD和BTM外，还有一些并未流行开来的衍射模型。一些方法使用海森堡不确定性原理解释声波衍射，使用衍射角概率密度函数描述衍射波沿非直线传播的特性。这种方法缺少物理基础，也难以保证结果的正确性；近来提出的方向性线源模型是一个有潜力的模型。它使用少量的近似解释了大量不同情况下的直边衍射，参考MENOUNOU P,NIKOLAOU P.Analytical model for predicting edge diffraction in the time domain.[J].The Journal of the Acoustical Society of America,2017,142 6:3580。但目前DLSM在实际声传播模拟器中的应用还没有出现。
因此，我们亟需一种几何框架下，大规模场景中计算高精度直边衍射的解决方案。
发明内容
本发明的目的在于针对现有技术的不足，提供一种基于路径跟踪的声波衍射现象模拟方法。
本发明是通过以下技术方案来实现的：
一种基于路径跟踪的声波直边衍射计算方法，包括如下步骤：
(1)读取与衍射相关的声学环境状态，所述声学环境状态包括声源与听者位置、场景几何描述、以及几何体表面的介面特征；
(2)在场景中随机生成大量衍射路径，所述路径由场景中一系列节点与连接节点的线段构成，路径的起始节点为声源，末尾节点为听者，其余中间节点位于场景中几何体的凸直边上；除各构成路径的节点外，路径的其它部分与场景中的几何体表面不相交；
(3)计算每一路径的最终贡献，设有路径包含n+1个节点，各节点从声源到听者依次编号为x 0到x n，记节点之间线段为x i→x i+1，完整路径为x 0→x 1→…→x n,则总贡献f(x 0→x 1→…→x n)计算公式如下：
Figure PCTCN2022136882-appb-000001
其中，L 0(x 0→x) 1为声源特性，M(x i→x i+1)为x i→x i+1处的声介质特性，ρ(x i-1→x i→x i+1)为x i处直边衍射的径向函数表示；
(4)将所有路径累加的贡献至冲激响应函数：对于路径x 0→x 1→…→x n，已知其贡献为步骤(3)中的f(x 0→x 1→…→x n)，并记其生成概率为p(x 0→x 1→…→x n)，传播延迟为t S，则在原有冲激响应函数之上累加
Figure PCTCN2022136882-appb-000002
其中δ(x)为单位冲激函数；
(5)输出步骤(4)产生的冲激响应函数作为衍射模拟结果。
具体地，所述步骤(2)包括如下子步骤：
(2.1)初始化路径，使其包含的唯一节点为路径起点；
(2.2)设当前路径状态为x 0→x 1→…→x n,从x n出发随机射出一条射线；若x n不是路径起点，则使用与直边衍射径向函数相匹配的采样函数决定射线的出射角度；
(2.3)记射线与场景几何表面的首个交点为
Figure PCTCN2022136882-appb-000003
并称其为伪交点；若该点不存在，则求交失败，转到步骤(2.8)；
(2.4)在步骤(2.3)中伪交点
Figure PCTCN2022136882-appb-000004
所在三角面的三条边中随机选择一条；设三角面为Δabc，被选择的边为ab，该边所对的顶点为c，则用直线连接c与
Figure PCTCN2022136882-appb-000005
并延长，与ab交于一 点，记为x n+1，称为交点；
(2.5)若步骤(2.4)中的边ab不是凸直边，或者线段x n→x n+1在两端点以外的地方与场景几何相交，则求交失败，转到步骤(2.8)；否则将交点x n+1作为节点添加至路径x 0→x 1→…→x n末尾，构成新的路径x 0→x 1→…→x n+1；
(2.6)计算路径x 0→x 1→…→x n+1的生成概率：若已知路径x 0→x 1→…→x n的生成概率为p(x 0→x 1→…→x n)，则路径x 0→x 1→…→x n+1的生成概率计算方法为
Figure PCTCN2022136882-appb-000006
其中，
Figure PCTCN2022136882-appb-000007
为x n处入射路径方向为x n-1→x n时，生成出射射线方向为
Figure PCTCN2022136882-appb-000008
的概率，若n＝0，则须将这一项替换为射线源出射射线方向为
Figure PCTCN2022136882-appb-000009
的概率；
Figure PCTCN2022136882-appb-000010
为射线
Figure PCTCN2022136882-appb-000011
与场景几何交于
Figure PCTCN2022136882-appb-000012
的概率；p C(x n→x n+1)为概率修正项；
(2.7)若路径长度未满足用户要求，则转到步骤(2.2)；
(2.8)将当前路径作为场景中的衍射路径。
进一步地，所述步骤(2)中，从声源和听者两端同时生成路径，从声源出发的路径称为前向子路径，从听者出发的路径称为后向子路径；之后，在前后向子路径的节点之间随机建立连接，以构成从声源到听者的完整路径；设前向子路径从声源到被连接处的部分生成概率为p F，后向子路径从听者到被连接处的部分的生成概率为p B，建立连接的概率为p L，则完整路径的生成概率为p F·p B·p L。
进一步地，所述步骤(5)中，路径的贡献需乘以一用于抑制极端值的因子C s：
Figure PCTCN2022136882-appb-000013
其中
Figure PCTCN2022136882-appb-000014
为一由用户定义的变量，C o为路径中随机生成的各段对应的C o之积，每一路径段的C o为衡量包含该路径段的路径产生极端值的可能性的一参数。
进一步地，所述步骤(3)中，计算贡献需除以正反向路径连接节点对应位置的[F j,i]中元素之积；所述[F j,i]为一个M行N列的数组，其中N为前向与后向路径的最大节点数；所 述路径的节点存储于一个行列数相同的数组当中，其中的每一行对应一条单独的路径，令数组[F j,i]第一列各元素为1，节点数组第一列为路径起点；在第i-1列节点已经生成的情况下，数组中的第i列节点的生成过程包括如下子步骤：
(5.1)从节点数组每一行的第i-1个节点x i-2出发，随机射出一条射线；若i＞1，则使用与直边衍射径向函数相匹配的采样函数决定射线的出射角度；
(5.2)记射线与场景几何表面的首个交点为
Figure PCTCN2022136882-appb-000015
并称其为伪交点；若该点不存在，则求交失败，转到步骤(5.7)；
(5.3)在步骤(5.2)中伪交点
Figure PCTCN2022136882-appb-000016
所在三角面的三条边中随机选择一条；设三角面为Δabc，被选择的边为ab，该边所对的顶点为c，则用直线连接c与
Figure PCTCN2022136882-appb-000017
并延长，与ab交于一点，记为x i-1，称为交点；
(5.4)若步骤(5.3)中的边ab非凸直边，或者线段x i-2→x i-1在两端点以外的地方与场景几何相交，则求交失败，转到步骤(5.7)；否则将交点x i-1作为节点添加至当前行的路径x 0→x 1→…→x i-2末尾，构成新的路径x 0→x 1→…→x i-1；
(5.5)计算路径x 0→x 1→…→ i-x的生成概率：若已知路径x 0→x 1→…→ i-x的生成概率为
Figure PCTCN2022136882-appb-000018
则路径
Figure PCTCN2022136882-appb-000019
的生成概率计算方法为:
Figure PCTCN2022136882-appb-000020
其中，
Figure PCTCN2022136882-appb-000021
为x i-2处入射路径方向为x i-3→x i-2时，生成出射射线方向为
Figure PCTCN2022136882-appb-000022
的概率，若i＝1，则须将这一项替换为射线源出射射线方向为
Figure PCTCN2022136882-appb-000023
的概率；
Figure PCTCN2022136882-appb-000024
为射线
Figure PCTCN2022136882-appb-000025
与场景几何交于
Figure PCTCN2022136882-appb-000026
的概率；p C(x i-2→x i-1)为概率修正项；
(5.6)若路径长度未满足用户要求，则转到步骤(5.1)，准备生成下一列节点；
(5.7)对于第i-1列的每个节点x i-2，统计出其后继节点生成成功的行集合S +(x i-2)和失败的行集合S -(x i-2)；
(5.8)将S +(x i-2)中各行第i列的节点复制进S -(x i-2)中各行的第i列；若第j行节点被复制k次，则令F j,i＝(k+1)F j,i-1。
本发明的有益效果是：
与传统方法相比，本文方法计算效率更高，且在处理大规模，遮挡关系复杂的场景时无需进行预求解波动方程，几何简化等预处理计算。本文方法使用的衍射模型无基尔霍夫近似，无限长直边假设等近似操作，因此结果精度更高，对于较简单几何体，其结果与理论衍射解完全符合。此外，基于路径跟踪的衍射算法是基于路径跟踪的传统声传播算法的自然扩展，因此本发明的方法也较易于添加至现有的声传播模拟器之中。
附图说明
图1是本发明中直边坐标系的示意图；
图2是本发明直边衍射介面项函数的绝对值图像；
图3是本发明直边衍射介面项采样概率分布的示意图；
图4是一阶衍射测试场景的布置示意图；
图5是本发明与对照方法在一阶衍射测试场景布置下的计算结果示意图；
图6是二阶衍射测试场景的布置示意图；
图7是本发明与对照方法在二阶衍射测试场景布置下的计算结果示意图。
具体实施方式
本发明提供了一种基于路径跟踪的声波直边衍射计算方法，包括如下步骤：
(1)读取与衍射相关的声学环境状态，所述声学环境状态包括声源与听者位置、场景几何描述、以及几何体表面的介面特征；
(2)在场景中随机生成大量衍射路径，所述路径由场景中一系列节点与连接节点的线段构成，路径的起始节点为声源，末尾节点为听者，其余中间节点位于场景中几何体的凸直边上；除各构成路径的节点外，路径的其它部分与场景中的几何体表面不相交；
(3)计算每一路径的最终贡献，设有路径包含n+1个节点，各节点从声源到听者依次编号为x 0到x n，记节点之间线段为x i→x i+1，完整路径为x 0→x 1→…→x n,则总贡献f(x 0→x 1→…→x n)计算公式如下：
Figure PCTCN2022136882-appb-000027
其中，L 0(x 0→x 1)为声源特性，M(x i→x i+1)为x i→x i+1处的声介质特性，ρ(x i-1→x i→x i+1) 为x i处直边衍射的径向函数表示；
(4)将所有路径累加的贡献至冲激响应函数：对于路径x 0→x 1→…→x n，已知其贡献为步骤(3)中的f(x 0→x 1→…→x n)，并记其生成概率为p(x 0→x 1→…→x n)，传播延迟为t S，则在原有冲激响应函数之上累加
Figure PCTCN2022136882-appb-000028
其中δ(x)为单位冲激函数；
(5)输出步骤(4)产生的冲激响应函数作为衍射模拟结果。
所述基于路径跟踪的声波直边衍射计算方法，具体为：
路径跟踪：声传播模拟中路径跟踪的技术是声传播模拟使用的路径跟踪算法以及衍射表示在算法中的运用。
在声传播模拟中，路径跟踪算法对下面的渲染方程进行求解：
Figure PCTCN2022136882-appb-000029
式中的t表示时间。t＝0的时刻为声源开始发声的时刻；L(x′→x,t)为声源x′在位置x处激发的能量响应或冲激响应，即本方程需求解的目标函数；L 0(x′→x,t)代表直射成分，即声源不与场景几何交互，直接在x处激发出的响应；V(x″,x)为可见性函数，在x″与x两位置互相可见时为1，反之为0；ρ(x′→x″→x)为介面项，表示入射声方向为x′→x″时，介面上一点x″附近在方向x″→x上产生的出射声强度；M t(x″,x′,t)为介质项，表示x″与x′之间声音介质的影响，包括能量吸收和传播延迟。星号*为卷积符号；Ω表示对应场景几何体表面的点集，A Ω为该集合上的面积测度。在计算衍射时，对于均匀介质，使用下面的介质项：
Figure PCTCN2022136882-appb-000030
其中a为介质的声衰减系数，||x″-x′||为向量x″-x′的长度，δ(·)为单位冲激响应函数，c为介质中的声速。
记上面的积分为算子T，则渲染方程可简记为L＝L 0+TL；使用诺依曼诺依曼(Neumann)级数展开来表达该方程：
Figure PCTCN2022136882-appb-000031
算子T表示了声波与场景几何的一次交互。由于T本质上是积分，其值可以通过蒙特卡洛积分法算出。
在路径跟踪中，本发明将T nL 0表示为路径贡献的数学期望；一条路径由场景中一系列节点与连接节点的线段构成，其起始节点为声源，末尾节点为听者，中间节点位于场景几何体的表面上；一条含有n个中间节点(n+2个节点)的路径称为n阶路径，与级数中的T nL 0一项相对应。
由介质项的表达式不难得知，一条路径的总传播延迟即路径各段时间延迟(长度除以声速)之和。因此可以单独计算出来。忽略时间项，仅计算路径贡献的强度，使用无时间参数的介质项
Figure PCTCN2022136882-appb-000032
替换渲染方程中原本的介质项M t。此时一条n阶路径x 0→x 1→…→x n+1的贡献表达式如下：
Figure PCTCN2022136882-appb-000033
若上面的路径为随机生成的，则计算其贡献的数学期望时需要除以其生成概率。再考虑时间延迟，则贡献期望可计算如下：
Figure PCTCN2022136882-appb-000034
其中t S为路径的总时间延迟。使用多条路径计算期望时，期望值为各路径结果的平均值。
生成随机路径的具体方法有很多种；使用单向路径跟踪或使用双向路径跟踪均可。
直边衍射的介面项表达：使用路径跟踪计算衍射的关键在于构造对应的介面项。这里本发明将给出介面项的形式。如图1所示，给出了直边附近的坐标系与入射，出射路径方向的表示方法及相关符号。本发明称图中的n为法线，b为副法线，t为切线(t方向与直边方向重合，正负皆可)。三者构成了直边上一点处的一个标架。直边处二面角的一半记为φ，由于本文仅考虑凸直边，φ的值应小于π/2。入射方向向量v i和出射方向向量v o使用一系列角度表示：
Figure PCTCN2022136882-appb-000035
Figure PCTCN2022136882-appb-000036
在这样的坐标系下，写出直边衍射的界面项形式；首先定义以下辅助变量：
Figure PCTCN2022136882-appb-000037
Figure PCTCN2022136882-appb-000038
ω lu＝θ i-π
ω lr＝-θ i-π+2φ
ω ru＝θ i+π
ω rr＝-θ i+π-2φ
Figure PCTCN2022136882-appb-000039
此时，若直边上一点附近的几何表面满足Dirichlet边界条件，则该点处的介面项为：
Figure PCTCN2022136882-appb-000040
若几何表面满足Neumann边界条件，则介面项为：
Figure PCTCN2022136882-appb-000041
上述介面项的形式可由数学证明。在其它条件下，可考虑采用形如aρ D+bρ N的介面项，其中a,b为二常数。
在路径跟踪中，需要从直边上一点发射射线，使其发射方向的概率分布接近介面项函数的绝对值。这一概率分布需要专门设计。为描述这一分布，需要定义一些辅助变量。令
Figure PCTCN2022136882-appb-000042
Figure PCTCN2022136882-appb-000043
使用Dirichlet边界条件时，令D＝D D；使用Neumann边界条件时，令D＝D N。其余变 量定义如下：
Figure PCTCN2022136882-appb-000044
Figure PCTCN2022136882-appb-000045
Figure PCTCN2022136882-appb-000046
Figure PCTCN2022136882-appb-000047
Figure PCTCN2022136882-appb-000048
那么，射线出射方向的概率密度函数
Figure PCTCN2022136882-appb-000049
便可以表示如下：
Figure PCTCN2022136882-appb-000050
其中，p θ是θ o的概率密度函数，
Figure PCTCN2022136882-appb-000051
是
Figure PCTCN2022136882-appb-000052
的概率密度函数。这两个函数的表达式如下：
Figure PCTCN2022136882-appb-000053
Figure PCTCN2022136882-appb-000054
上面构造出的函数满足概率密度函数的非负性和归一性要求。
图2展示了在Dirichlet条件下，
Figure PCTCN2022136882-appb-000055
θ i＝0.4636，
Figure PCTCN2022136882-appb-000056
时的介面项绝对值，图3为本发明直边衍射介面项采样概率分布的示意图，展示了对应方向下的采样概率密度函数。可以看到这两图中函数的形态是非常相似的。
衍射路径的生成：在路径跟踪算法中，常见的路径生成方法是依次生成路径中的各个节点。一般是从当前路径的末尾节点射出一条方向随机的射线，若其与场景几何相交，则将交点作为原末尾节点的后继节点。然而在模拟衍射时，路径的中间节点必须位于几何体的边上，而一般的射线求交并不能保证交点满足这一要求。因此需要对其进行改动。
假设场景几何由三角网格描述。在射线求交成功之后，称射线与几何体的交点为伪交点。检查伪交点所在的三角面的各条边，并判断这些边是否为凸边。接下来在这些凸边中随机选择一条。设三角面为Δabc 1，被选择的边为ab，则记选中该边的概率为
Figure PCTCN2022136882-appb-000057
在此之后，使用直线连接伪交点与c 1点并延长，与ab交于一点。称其为伪交点对应的(实际)交点，并 将它作为原末尾节点的后继节点。
现计算这种方法下节点的生成概率。记原本的路径为x 0→x 1→…→x n，伪交点为
Figure PCTCN2022136882-appb-000058
对应的交点为x n+1。ab有两个相邻面，记其中与Δabc 1不同的另一个面为Δabc 2。那么，所有对应x n+1的伪交点均位于线段x n+1c 1和x n+1c 2之上。接下来，计算两个辅助变量S(x n,x n+1c 1)和S(x n,x n+1c 2)。S(x n,x n+1c 1)的计算方法如下：
(3.1)找出x n+1c 1对于x n可见(未被遮挡)的部分。由于场景由三角面片组成，可见部分必然是一系列不相交线段的集合；通过找出与场景几何中三角形Δx nx n+1c 1相交的三角面片来算得所有遭到遮挡的线段，然后从x n+1c 1中去除这些被遮挡的线段，余下的就是可见的线段集合。记这些可见线段为
Figure PCTCN2022136882-appb-000059
其中
Figure PCTCN2022136882-appb-000060
相比
Figure PCTCN2022136882-appb-000061
更靠近c 1一端。
(3.2)计算S(x n,x n+1c 1)。公式为
Figure PCTCN2022136882-appb-000062
S(x n,x n+1c 2)的计算方法与上面相同。
最终，给出节点生成概率的一个估计：
Figure PCTCN2022136882-appb-000063
其中，
Figure PCTCN2022136882-appb-000064
为x n处入射路径方向为x n-1→x n时，生成出射方向
Figure PCTCN2022136882-appb-000065
的概率，具体概率分布由算法实现者决定。若n＝0，则须将这一项替换为射线源出射射线方向为
Figure PCTCN2022136882-appb-000066
的概率。
Figure PCTCN2022136882-appb-000067
为出射射线与场景交于
Figure PCTCN2022136882-appb-000068
的概率，其计算公式如下：
Figure PCTCN2022136882-appb-000069
其中
Figure PCTCN2022136882-appb-000070
为射线
Figure PCTCN2022136882-appb-000071
与
Figure PCTCN2022136882-appb-000072
所在三角面法线的夹角。p C(x n→x n+1)为概率修正项，计算公式如下：
Figure PCTCN2022136882-appb-000073
其中
Figure PCTCN2022136882-appb-000074
为Δabc 1的面积，
Figure PCTCN2022136882-appb-000075
同理。由于x n+1对应的
Figure PCTCN2022136882-appb-000076
不止一个，上面的概率本身并 不直接等同于生成节点x n+1的生成概率，但其期望确实与生成x n+1的概率相等。
增加节点后的路径生成概率如下：
Figure PCTCN2022136882-appb-000077
所述双向路径跟踪中，对于双向路径跟踪需要从声源和听者两端同时生成路径，然后在两条路径的节点之间建立连接，生成完整路径。此时完整路径的生成概率为p F·p B·p L。其中p F为前向子路径从声源到被连接处的部分生成概率，p B为后向子路径从听者到被连接处的部分的生成概率，p L为建立连接的概率。
此外，本发明还可采用多重重要性采样，它为是一种提高样本质量的方法，在双向路径跟踪中通常会配合使用。对于一条路径，多重重要性采样要求给出这条路径在所有生成方式下的生成概率。假设通过双向路径跟踪产生了一条路径x 0→x 1→…→x n。它的连接段可能是路径中的任何一段，因此共有n-1中可能的生成方式。这些生成方式下路径生成概率的求法。在声波直边衍射计算基础之上结合多重重要性采样技术后，为计算多重重要性采样技术中各种策略的概率，须为每个非起点或终点的节点生成对应的前向和后向伪交点，其中前向伪交点对该节点的前趋节点可见，后向伪交点对该节点的后继节点可见；具体如下：
设路径x 0→x 1→…→x n由双向路径跟踪生成，其连接段为x i→x i+1，那么这条路径的生成概率为p(x 0→x 1→…→x i)p(x i+1←x i+2←…←x n)p L(x i→x i+1)。这里x i+1←x i+2←…←x n的含义等同于x n→x n-1→…→x i+1。将计算路径生成概率所涉及到的场景位置列出如下：
Figure PCTCN2022136882-appb-000078
这里类似于
Figure PCTCN2022136882-appb-000079
的点与
Figure PCTCN2022136882-appb-000080
一样，是伪交点，称之为“后向伪交点”。它们对应的是反向路径上的节点。由反向路径的生成过程可知，
Figure PCTCN2022136882-appb-000081
对x i+1可见，但不能保证对x i-1可见，与
Figure PCTCN2022136882-appb-000082
这样的伪交点(称为前向伪交点)正好相反。
改变上面图示中i的值可知，为了求得连接段为其他路径段时的生成概率，每一个中间节点都需要有对应的前向和后向伪交点。因此需要主动补齐图示中缺失的伪交点。这里将给出补齐前向伪交点的方法，补齐反向伪交点的方法与之相似。
假设x i+1没有对应的前向伪交点，首先寻找x i+1所在的边，此处记其为ab。与这条边相邻的两个邻面记为Δabc 1和Δabc 2。再找出线段x i+1c 1和x i+1c 2相对于x i可见的部分，并将其表示为线段的集合，每条线段使用
Figure PCTCN2022136882-appb-000083
表示。定义函数：
Figure PCTCN2022136882-appb-000084
Figure PCTCN2022136882-appb-000085
接下来，对每条线段
Figure PCTCN2022136882-appb-000086
若其为x i+1c 1的子集，则通过
Figure PCTCN2022136882-appb-000087
将其映射到集合[0,1]中；若其为x i+1c 2的子集，则通过
Figure PCTCN2022136882-appb-000088
将其映射到集合[0,1]中；其次，在映射后的线段集合中按照均匀分布随机选取一点；最后通过f lerp将选取的这一点映射回原本的线段集合，并以之作为伪交点。
补全伪交点后，便可算出所有生成方式下的路径生成概率。
此外，本发明也可实用其他优化：
当生成路径概率过小，或路径贡献相对生成概率过大时，结果中会出现可能影响用户体验的极端值。因此需要识别出这些极端值并去除之。
给定一条路径x 0→x 1→…→x n+1，定义辅助变量如下：
Figure PCTCN2022136882-appb-000089
其中p V与
Figure PCTCN2022136882-appb-000090
的定义与第三节中相同。若i＝0，p V也需要按照第三章中的方法替换。这一变量表示了包含该段的路径产生极端值的可能性。整条路径的C o值为路径各段的C o值之积。若路径由多段子路径连接而成(如双向路径跟踪)，则其C o值为各段子路径的C o值之积。接下来，在路径的最终贡献上乘以一个值C s:
Figure PCTCN2022136882-appb-000091
这里的
Figure PCTCN2022136882-appb-000092
为一由用户定义的变量(通常在100-1000之间即可)。上面的乘数将压低
Figure PCTCN2022136882-appb-000093
的所有路径的贡献，从而减少极端值的产生。
若生成路径后继节点方法的效率并不够高，在双向路径跟踪中，为了提高节点的利用效 率，可以使用下面的方法：首先假设模拟算法中一次生成的路径有M条，路径的最大阶数为N-2。因此这些路径共有M×N个节点。这些节点，无论实际是否被成功生成，均排列于一个M行N列的数组当中。对于数组中的每个节点，定义一个变量F j,i,j＝1,2, …M,i＝1,2,…N，称为拆分度，并令F j,1＝1。数组中的节点需要逐行生成。对于数组中的第i列节点，具体生成过程如下：尝试通过第三节中的方法生成第i列中的所有节点；对于第i-1列的每个节点x，统计出其后继节点生成成功的行集合S +(x)和失败的行集合S -(x)；将S +(x)中各行第i列的节点复制进S -(x)中各行的第i列。若第j行节点被复制k次，则令F j,i＝(k+1)F j,i-1。设路径到此节点的贡献期望为L j,i，建议各节点的复制次数与L j,i/F j,i-1成正比。
在建立正反向路径连接，生成最终路径时，应当让最终路径的贡献期望除以正反向路径连接节点对应的两个拆分度之积。通过这种方式，充分利用数组当中的空位，节省存储空间，提高节点连接效率。
实施实例
发明人在一台配备Intel i7 2.60GHz四核中央处理器，12GB内存的计算机上实现了上文所描述的算法，具体技术选择为结合多重重要性采样的双向路径跟踪，并使用了实施过程描述第五节中的所有优化技巧。本发明将其结果与使用UTD和BTM模型的模拟结果进行了比较。
图4与图6展示了用于比较的一阶和二阶衍射场景。图4中的场景包含一个直边长度为2m的直二面角，图6中的场景则包含两个平行的直二面角，从声源发出的路径必须经过两次衍射，才能到达听者。两场景几何均使用狄利克雷(Dirichlet)边界条件。图5与图7展示了场景中使用不同方法算得的冲激响应，此处使用的方法包括BTM、UTD、本文方法以及选取UTD中8个频点(位于62.5-8000Hz之间，以指数间隔方式分布)重构出的结果UTD方法计算二阶及以上衍射误差甚大，在图7中未展示。可以看到，本文算法在一至二阶衍射情况下的计算结果质量与已有方法没有明显的差异，且明显优于常见的简化UTD方案(仅计算8个频率上的结果)，由于BTM方法对于一阶衍射的解是精确的，这一比较也验证了本发明中方法的正确性。对于大规模复杂场景，本方法展现出了接近实时的计算效率。在使用12000条路径的情况下，对于包含22974三角面的模型，本发明的实现计算衍射冲激响应的时间开销为约140ms。对于包含279133面的模型，时间开销则为约205ms。由此可见，本发明方法适合在复杂场景中进行快速的衍射结果计算。
以上所述实施例表达了本发明的具体实施方式，其描述较为具体和详细，旨在用于帮助 理解本发明的方法及其核心思想，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (5)

  1. 一种基于路径跟踪的声波直边衍射计算方法，其特征在于，包括如下步骤：
    (1)读取与衍射相关的声学环境状态，所述声学环境状态包括声源与听者位置、场景几何描述、以及几何体表面的介面特征；
    (2)在场景中随机生成大量衍射路径，所述路径由场景中一系列节点与连接节点的线段构成，路径的起始节点为声源，末尾节点为听者，其余中间节点位于场景中几何体的凸直边上；除各构成路径的节点外，路径的其它部分与场景中的几何体表面不相交；
    (3)计算每一路径的最终贡献，设有路径包含n+1个节点，各节点从声源到听者依次编号为x 0到x n，记节点之间线段为x i→x i+1，完整路径为x 0→x 1→…→x n,则总贡献f(x 0→x 1→…→x n)计算公式如下：
    Figure PCTCN2022136882-appb-100001
    其中，L 0(x 0→x) 1为声源特性，M(x i→x i+1)为x i→x i+1处的声介质特性，ρ(x i-1→x i→x i+1)为x i处直边衍射的径向函数表示；
    (4)将所有路径累加的贡献至冲激响应函数：对于路径x 0→x 1→…→x n，已知其贡献为步骤(3)中的f(x 0→x 1→…→x n)，并记其生成概率为p(x 0→x 1→…→x n)，传播延迟为t S，则在原有冲激响应函数之上累加
    Figure PCTCN2022136882-appb-100002
    其中δ(x)为单位冲激函数；
    (5)输出步骤(4)产生的冲激响应函数作为衍射模拟结果。
  2. 根据权利要求1所述的一种基于路径跟踪的声波直边衍射计算方法，其特征在于，所述步骤(2)包括如下子步骤：
    (2.1)初始化路径，使其包含的唯一节点为路径起点；
    (2.2)设当前路径状态为x 0→x 1→…→x n,从x n出发随机射出一条射线；若x n不是路径起点，则使用与直边衍射径向函数相匹配的采样函数决定射线的出射角度；
    (2.3)记射线与场景几何表面的首个交点为
    Figure PCTCN2022136882-appb-100003
    并称其为伪交点；若该点不存在，则求交失败，转到步骤(2.8)；
    (2.4)在步骤(2.3)中伪交点
    Figure PCTCN2022136882-appb-100004
    所在三角面的三条边中随机选择一条；设三角面为Δabc，被选择的边为ab，该边所对的顶点为c，则用直线连接c与
    Figure PCTCN2022136882-appb-100005
    并延长，与ab交于一点，记为x n+1，称为交点；
    (2.5)若步骤(2.4)中的边ab不是凸直边，或者线段x n→x n+1在两端点以外的地方与场景几何相交，则求交失败，转到步骤(2.8)；否则将交点x n+1作为节点添加至路径x 0→x 1→…→x n末尾，构成新的路径x 0→x 1→…→x n+1；
    (2.6)计算路径x 0→x 1→…→x n+1的生成概率：若已知路径x 0→x 1→…→x n的生成概率为p(x 0→x 1→…→x n)，则路径x 0→x 1→…→x n+1的生成概率计算方法为
    Figure PCTCN2022136882-appb-100006
    其中，
    Figure PCTCN2022136882-appb-100007
    为x n处入射路径方向为x n-1→x n时，生成出射射线方向为
    Figure PCTCN2022136882-appb-100008
    的概率，若n＝0，则须将这一项替换为射线源出射射线方向为
    Figure PCTCN2022136882-appb-100009
    的概率；
    Figure PCTCN2022136882-appb-100010
    为射线
    Figure PCTCN2022136882-appb-100011
    与场景几何交于
    Figure PCTCN2022136882-appb-100012
    的概率；p C(x n→x n+1)为概率修正项；
    (2.7)若路径长度未满足用户要求，则转到步骤(2.2)；
    (2.8)将当前路径作为场景中的衍射路径。
  3. 根据权利要求1所述的一种基于路径跟踪的声波直边衍射计算方法，其特征在于，所述步骤(2)中，从声源和听者两端同时生成路径，从声源出发的路径称为前向子路径，从听者出发的路径称为后向子路径；之后，在前后向子路径的节点之间随机建立连接，以构成从声源到听者的完整路径；设前向子路径从声源到被连接处的部分生成概率为p F，后向子路径从听者到被连接处的部分的生成概率为p B，建立连接的概率为p L，则完整路径的生成概率为p F·p B·p L。
  4. 根据权利要求1所述的一种基于路径跟踪的声波直边衍射计算方法，其特征在于，所述步骤(5)中，路径的贡献需乘以一用于抑制极端值的因子C s：
    Figure PCTCN2022136882-appb-100013
    其中
    Figure PCTCN2022136882-appb-100014
    为一由用户定义的变量，C o为路径中随机生成的各段对应的C o之积，每一路径段的C o为衡量包含该路径段的路径产生极端值的可能性的一参数。
  5. 根据权利要求3所述的一种基于路径跟踪的声波直边衍射计算方法，其特征在于，所述步骤(3)中，计算贡献需除以正反向路径连接节点对应位置的[F j,i]中元素之积；所述[F j,i]为一个M行N列的数组，其中N为前向与后向路径的最大节点数；所述路径的节点存储于一个行列数相同的数组当中，其中的每一行对应一条单独的路径，令数组[F j,i]第一列各元素为1，节点数组第一列为路径起点；在第i-1列节点已经生成的情况下，数组中的第i列节点的生成过程包括如下子步骤：
    (5.1)从节点数组每一行的第i-1个节点x i-2出发，随机射出一条射线；若i＞1，则使用与直边衍射径向函数相匹配的采样函数决定射线的出射角度；
    (5.2)记射线与场景几何表面的首个交点为
    Figure PCTCN2022136882-appb-100015
    并称其为伪交点；若该点不存在，则求交失败，转到步骤(5.7)；
    (5.3)在步骤(5.2)中伪交点
    Figure PCTCN2022136882-appb-100016
    所在三角面的三条边中随机选择一条；设三角面为Δabc，被选择的边为ab，该边所对的顶点为c，则用直线连接c与
    Figure PCTCN2022136882-appb-100017
    并延长，与ab交于一点，记为x i-1，称为交点；
    (5.4)若步骤(5.3)中的边ab非凸直边，或者线段x i-2→x i-1在两端点以外的地方与场景几何相交，则求交失败，转到步骤(5.7)；否则将交点x i-1作为节点添加至当前行的路径x 0→x 1→…→x i-2末尾，构成新的路径x 0→x 1→…→x i-1；
    (5.5)计算路径
    Figure PCTCN2022136882-appb-100018
    的生成概率：若已知路径
    Figure PCTCN2022136882-appb-100019
    的生成概率为
    Figure PCTCN2022136882-appb-100020
    则路径
    Figure PCTCN2022136882-appb-100021
    的生成概率计算方法为:
    Figure PCTCN2022136882-appb-100022
    其中，
    Figure PCTCN2022136882-appb-100023
    为x i-2处入射路径方向为x i-3→x i-2时，生成出射射线方向为
    Figure PCTCN2022136882-appb-100024
    的概率，若i＝1，则须将这一项替换为射线源出射射线方向为
    Figure PCTCN2022136882-appb-100025
    的概率；
    Figure PCTCN2022136882-appb-100026
    为射线
    Figure PCTCN2022136882-appb-100027
    与场景几何交于
    Figure PCTCN2022136882-appb-100028
    的概率；p C(x i-2→x i-1)为概率修正项；
    (5.6)若路径长度未满足用户要求，则转到步骤(5.1)，准备生成下一列节点；
    (5.7)对于第i-1列的每个节点x i-2，统计出其后继节点生成成功的行集合S +(x i-2)和失败的行集合S -(x i-2)；
    (5.8)将S +(x i-2)中各行第i列的节点复制进S -(x i-2)中各行的第i列；若第j行节点被复制k次，则令F j,i＝(k+1)F j,i-1。
PCT/CN2022/136882 2022-12-06 2022-12-06 一种基于路径跟踪的声波直边衍射计算方法 Ceased WO2024119366A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/CN2022/136882 WO2024119366A1 (zh) 2022-12-06 2022-12-06 一种基于路径跟踪的声波直边衍射计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2022/136882 WO2024119366A1 (zh) 2022-12-06 2022-12-06 一种基于路径跟踪的声波直边衍射计算方法

Publications (1)

Publication Number Publication Date
WO2024119366A1 true WO2024119366A1 (zh) 2024-06-13

Family

ID=91378437

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2022/136882 Ceased WO2024119366A1 (zh) 2022-12-06 2022-12-06 一种基于路径跟踪的声波直边衍射计算方法

Country Status (1)

Country Link
WO (1) WO2024119366A1 (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006208273A (ja) * 2005-01-31 2006-08-10 Asahi Kasei Homes Kk 音場シミュレーションシステム
US20150378019A1 (en) * 2014-06-27 2015-12-31 The University Of North Carolina At Chapel Hill Methods, systems, and computer readable media for modeling interactive diffuse reflections and higher-order diffraction in virtual environment scenes
CN105893719A (zh) * 2016-06-16 2016-08-24 浙江大学 一种基于双向路径跟踪的实时声音传播模拟方法
CN114003855A (zh) * 2021-10-20 2022-02-01 国网重庆市电力公司电力科学研究院 基于声线跟踪理论的声波障碍物绕射模拟方法及存储介质

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006208273A (ja) * 2005-01-31 2006-08-10 Asahi Kasei Homes Kk 音場シミュレーションシステム
US20150378019A1 (en) * 2014-06-27 2015-12-31 The University Of North Carolina At Chapel Hill Methods, systems, and computer readable media for modeling interactive diffuse reflections and higher-order diffraction in virtual environment scenes
CN105893719A (zh) * 2016-06-16 2016-08-24 浙江大学 一种基于双向路径跟踪的实时声音传播模拟方法
CN114003855A (zh) * 2021-10-20 2022-02-01 国网重庆市电力公司电力科学研究院 基于声线跟踪理论的声波障碍物绕射模拟方法及存储介质

Similar Documents

Publication Publication Date Title
Aleksandrov et al. Determining approximate shortest paths on weighted polyhedral surfaces
Kulla et al. Importance sampling techniques for path tracing in participating media
US10679407B2 (en) Methods, systems, and computer readable media for modeling interactive diffuse reflections and higher-order diffraction in virtual environment scenes
Mehra et al. Wave-based sound propagation in large open scenes using an equivalent source formulation
Rong et al. GPU-assisted computation of centroidal Voronoi tessellation
Drumm et al. The adaptive beam-tracing algorithm
Melvær et al. Geodesic polar coordinates on polygonal meshes
CN115453753A (zh) 基于nurbs表面的光学系统高精度光线追迹方法及装置
CN104361624A (zh) 一种计算机三维模型中全局光照的渲染方法
Antonacci et al. Fast tracing of acoustic beams and paths through visibility lookup
Zou et al. An algorithm for triangulating multiple 3D polygons
CN115935673A (zh) 一种基于射线追踪的fdtd网格剖分方法及系统
Carr et al. Scalable contour tree computation by data parallel peak pruning
Zhao et al. Progressive discrete domains for implicit surface reconstruction
CN105893719B (zh) 一种基于双向路径跟踪的实时声音传播模拟方法
CN102722653B (zh) 一种基于MapReduce的射线跟踪加速算法
CN116385623A (zh) 具有深度信息的参与介质多重散射的绘制方法及系统
Meng et al. Point-based acoustic scattering for interactive sound propagation via surface encoding
US20260059258A1 (en) Specular reflection path generation and near-reflective diffraction in interactive acoustical simulations
CN117076827A (zh) 一种基于路径跟踪的声波直边衍射计算方法
Siltanen et al. Room acoustics modeling with acoustic radiance transfer
CN118247317A (zh) 一种基于新点到面和自适应鲁棒损失的点云精配准方法
Shapot et al. Solution building for arbitrary system of linear inequalities in an explicit form
Yue et al. Toward optimal space partitioning for unbiased, adaptive free path sampling of inhomogeneous participating media
Martinelli An application of the level set method to underwater acoustic propagation

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22967528

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 22967528

Country of ref document: EP

Kind code of ref document: A1