一种融合约束的基于位置动力学的实时雪模拟方法

本发明涉及计算机图形学和虚拟现实仿真领域，具体涉及一种融合约束的基于位置动力学的实时雪模拟方法。总体框架示意如图1所示，其中主要通过基于松弛因子的非渗透约束和摩擦约束矫正粒子位置，克服显式计算物理力导致的仿真效率低下的问题，通过融合希尔伯特曲线编码提升邻域搜索效率，最终得到实时的符合物理规律的粒子雪仿真效果。

背景技术：

1、雪是户外场景中不可或缺的重要大气自然现象，雪仿真(snow simulation)是计算机虚拟现实仿真领域的重要课题，其真实感模拟仿真有助于创造出令人信服的虚拟环境，对于增加沉浸感和体验感具有重要意义。雪的模拟可分为静态雪模拟和动态雪模拟两个方面。

2、静态雪模拟的重点是获得逼真的雪景图像效果。需要仔细模拟雪的分布状态、形态特征以及与灯光的相互作用，从而在视觉上传达令人信服、栩栩如生的静态雪景外观。这对提升虚拟场景的真实感和细节层次起到关键作用。但这种静态模拟往往基于良好的渲染框架，在目前成熟的渲染体系里已有良好的效果，其改进创新和拓展空间并不多。

3、相比之下，动态雪模拟面临更高的挑战。除了需要准确模拟雪的形成过程外，还必须自然地再现雪的运动轨迹、积雪逐渐堆积的动态变化过程。动态雪模拟不仅需要精确计算雪粒运动路径，还要考虑雪粒相互作用、与物体碰撞等复杂情况才能获得高度真实、连贯自然的雪景动画效果。

4、针对不同的应用场景，动态雪模拟技术在计算效率方面也有不同的要求。例如，对于电影特效，渲染雪的计算时间可能是离线的；对于游戏、虚拟模拟等实时应用，由于需要在保证保真度的同时仍满足现实条件，因此对雪模拟的效率要求更高的时间响应。为此，模拟技术必须在视觉质量和计算效率之间取得平衡。具体来说，动态雪的模拟方法大致可分为基于非物理的方法和基于物理的方法两大类。

5、在游戏视频等领域的传统的雪仿模拟仿真大多基于非物理的方法，这些方法可以仿真出达到实时速率要求并且满足一定视觉效果的结果。然而，其技术核心仍然集中于通过贴图进行纹理映射并操控网格达到一定程度上的交互，这些方法需要大量的美术人员参与其中，不仅效率低而且无法展现真正的雪的物理特性。与此同时，该方法无法满足日益提升的人们的审美需求，更无法提供一个令人信服的具有沉浸感的虚拟自然环境。

6、特别的，基于物理方法模拟雪与本发明的工作最为相关，对其相关研究背景依次介绍如下。

7、(1)基于材质点法的雪模拟的方法

8、材质点法(material point method，mpm)是一种基于混合欧拉-拉格朗日视角的方法，它的基本思路流程是：

9、1.根据粒子的分布更新欧拉网格(particle to grid，p2g)；

10、2.在欧拉网格中进行投影；

11、3.根据欧拉网格的信息更新粒子的状态(grid to particle，g2p)；

12、4.从拉格朗日粒子的角度进行步进。

13、p2g的过程是将粒子信息映射到网格上，亦或是从一组粒子重建整个场的信息。通常使用一个径向对称的核函数来确定一个粒子对周边空间的影响，即一种加权计算；其中粒子的影响是局部化的，如图2所示核函数通常选择有限支撑，即存在一个最大支撑距离。

14、g2p是在计算完网格上各种场的对应信息之后，将周围网格的信息叠加到原粒子上。然而，在g2p的过程中，粒子信息是由网格信息叠加而成的，存在大量的信息缺失(能量损耗)，主要是由于在旋转和切变的过程中角动量损耗及其严重。为了解决信息丢失，除了位置和速度，apic(affine particle-in-cell)方法为每个粒子增加了一个仿射(affine)矩阵，该矩阵记录粒子的散度和角速度。

15、2013年，stomakhin等人首次提出基于apic法为粒子额外加入形变梯度f、体积应变比j、和杨氏模量e等属性并应用于雪的模拟仿真。该方法的本构模型结合了压缩性、塑性和硬化性，自动处理了断裂和堆积雪的特性粘性效应。是一种利用用户可控弹塑性本构模型与欧拉/拉格朗日混合材料点法集成的新型雪模拟方法。该方法是基于连续体的，其混合性质允许使用常规的笛卡尔网格来自动处理自碰撞和断裂。虽然这种方法许多参数都是手动调整的，而且忽略了与空气的相互作用，难以模拟像粉雪和雪崩这样的材质，但却是基于物理方法的雪模拟极具意义的提出。(参考文件1：stomakhin a,schroeder c,chai l,etal.a material point method for snow simulation[j].acm transactions ongraphics(tog),2013,32(4):1-10.)

16、2015年，gast等人展示了如何在不牺牲效率的情况下，通过约束最小化将简单的碰撞直接合并到求解器中。还提出了针对自碰撞和脚本主体碰撞的新颖惩罚碰撞公式，这些公式专为满足此求解器的独特需求而设计。(参考文件2：gast t f,schroeder c,stomakhin a,et al.optimization integrator for large time steps[j].ieeetransactions on visualization and computer graphics,2015,21(10):1103-1115.)

17、2018年，胡渊明等人介绍了移动最小二乘材料点方法(mls-mpm)，模拟了mpm以前不支持的各种新现象，包括材料切割、动态开放边界以及与刚体的双向耦合。(参考文件3：hu y,fang y,ge z,et al.a moving least squares material point method withdisplacement discontinuity and two-way rigid body coupling[j].acmtransactions on graphics(tog),2018,37(4):1-14.)

18、2019年，han等人利用拉格朗日网格进行内力计算，利用欧拉网格进行自碰撞以及与外部材料耦合，在mpm模拟中显示了摩擦接触，尤其是雪与头发之间的摩擦接触。(参考文件4：han x,gast t f,guo q,et al.a hybrid material point method for frictionalcontact with diverse materials[j].proceedings of the acm on computer graphicsand interactive techniques,2019,2(2):1-24.)

19、(2)基于平滑粒子流体力学的雪模拟方法

20、在流体力学中，平滑粒子流体力学(smoothed particle hydrodynamics，sph)算法的基本设想是将连续的流体想象成一个个相互作用的微粒，这些粒子相互影响，共同形成了复杂的流体运动，对于每个单独的流体微粒，依旧遵循最基本的牛顿第二定律。流体中某点r，在光滑核半径h范围内有数个粒子，r0,r1,r2,…rj，则该处某项属性ai(r)在空间域ω上由sph方法的核函数w积分求解如式(1)所示：

21、ai(r)＝∫ωajw(r-rj,h)drj                   (1)

22、实际上，在计算机求解积分时，通常采用将被积函数离散在空间域内的不同子域vj上并通过累加的方式进行迭代计算，又由于v＝m/ρ，因此式(1)可以转换为如式(2)所示：

23、

24、其中mj和ρj是粒子j的质量和密度。核函数w是一种随着距离而衰减的函数，它描述粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，最大影响半径为核半径h。sph方法在问题域上求解的基本流程就是根据光滑核函数逐个推出流体中某点的属性ai(r)，例如密度、压力、速度相关的累加结果，进而推导出此处的加速度，从而模拟流体的运动趋势。

25、1996年，cg(computer graphics，计算机图形学)相关的研究者desbrun等人最早将sph法引入计算机模拟“软体”的物理效果。该团队使用平滑的颗粒代表样品点，可以实现介质内局部物理量的值和衍生物的近似值，还可以确保对描述材料物理行为的状态方程的有效稳定模拟。最终，模拟出具有多种刚度和粘度的非弹性物体的动画。(参考文件5：desbrun m,gascuel m p.smoothed particles:anew paradigm for animating highlydeformable bodies[c]//computer animation and simulation’96:proceedings of theeurographics workshop in poitiers,france,august 31–september 1,1996.vienna:springer vienna,1996:61-76.)

26、雪被归类为非牛顿流体是一种弹塑性材质，虽然不能由牛顿流体的原始n-s方程控制，但是可以通过引入其他控制项来改进。2012年，takahashi等人基于sph方法将雪模拟为流体，修改了n-s方程的黏度项来建模雪，是雪模拟仿真的一项突破。于此同时，他们通过将熔结项(sintering)附加在n-s方程中使得雪可以黏附在刚体上。然而由于他们采用黏度模型，其雪模拟仅对剪切速率作出反应，无法真实地捕捉复杂的物理雪行为。(参考文件6：takahashi t,fujishiro i.particle-based simulation of snow trampling takingsintering effect into account[m]//acm siggraph 2012posters.2012:1-1.)

27、2020年，gissler等人提出两个隐式求解器来计算雪的弹性形变。这两个求解器的组合允许模拟雪及雪与其他相(如流体)的相互作用。他们将单个雪片建模为单个粒子，使用sph方法进行离散化，从而简化了降雪和积雪的模拟。此外，作者还在基于粒子的边界处提出了隐式摩擦公式。不过，单个雪片或冰晶在第一次碰撞表面后往往会弹开，而gissler等人所提出方法的边界处的隐式摩擦公式还无法参数化单个粒子的反弹效应。(参考文件7：gissler c,henne a,band s,et al.an implicit compressible sph solver for snowsimulation[j].acm transactions on graphics(tog),2020,39(4):36:1-36:16.)

28、(3)基于离散元法的雪模拟方法

29、离散元法(discrete element method，dem)广泛用于模拟颗粒材料即离散粒子，可以假设颗粒材料是固体颗粒的集合。与sph的离散化思想不同的是，dem将求解空间离散为单元阵，主要计算两个离散单元之间的法向力和切向力。之后根据合力以及牛顿运动第二定律求得离散单元的加速度，在时间上积分进而得到单离散元的速度和位移。

30、2015年，hagenmuller等人使用dem对变形下的硬化雪进行了颗粒描述，来开发基于微观学捕获的完整3-d微结构的离散元件模型。该模型假设雪是由通过局部触点相互作用的刚性晶粒组成的，这些刚性颗粒通过局部接触相互作用，模拟了单轴约束压缩，同时使用多个变形阶段来确定雪的弹性或塑性行为。(参考文件8：hagenmuller p,chambon g,naaim m.microstructure-based modeling of snow mechanics:a discrete elementapproach[j].the cryosphere,2015,9(5):1969-1982.)

31、goswami等人也使用dem法进行离散模拟仿真雪的压缩、粘结行为并使用gpu加速。该研究通过迭代方案矫正雪物理求解器中的粒子位置，根据每个粒子大小或压缩级别计算每个粒子密度迭代求最小化粒子重叠。(参考文件9：goswami p,nordin a,nyléns.iterative discrete element solver for efficient snow simulation[c]//eurographics symposium on parallel graphics and visualization(egpgv),rome,june 13.2022.)

32、综上所述，目前大多数基于物理的雪仿真方法存在着过度关注物理力计算导致模拟效率低下的问题，而对于不显式计算物理力直接通过位置矫正的方法模拟雪却相对较少。本发明主要解决这一问题，该方法不仅可以较为准确地还原雪的物理特性，还可以在满足实时效率的情况下模拟大规模粒子组成的雪的运动行为。

技术实现思路

1、本发明旨在克服复杂的物理方程导致显式计算物理力模拟雪的物理行为效率低下的问题，通过设计基于松弛因子的非渗透约束并融合摩擦力约束来矫正雪粒子的位置来避免直接计算物理力，从而保证在减轻模拟计算压力的基础上真实地还原雪在下落时的运动规律。通过融合希尔伯特曲线的哈希编码提高邻域搜索的效率从而提高整体模拟帧率，整体思路如图1所示。

2、在基于位置动力学(position based dynamic，pbd)中，如图3所示，基本框架首先初始化所有粒子的位置、速度和质量，然后施加外力包括重力并对速度阻尼；接下来预测新的位置；最后在约束方程的约束下对预测的位置进行矫正，更新到最新的粒子位置同时更新速度，主要是防止速度震荡。

3、约束本质是为了模拟内力其稳定度估值是一个方程，入参是所有被该约束影响的粒子，返回值是一个浮点数表示稳定的程度。约束方程可以是非线性不等式方程也可以是等式方程，通常有两种表示：

4、(1)等式表示：

5、(2)不等式表示：

6、其中，j∈[1,…,m]代表有m个约束，i∈[1,…,n]代表有n个粒子。预测的位置并不一定能满足所有约束的条件。即本来整个系统稳定，或者说所有约束的估值和为0，但是当预测后系统开始不稳定，也即是估值函数返回值不为零。此时约束就要开始对粒子预测的位置进行矫正，目标是期望让整个系统稳定下来。用p代替所有被该约束影响的粒子，那么约束方程可以表示为c(p)。p是原状态也就是原位置，那么p+δp就是下一时刻的状态，约束方程可以用泰勒展开式近似如式(3)所示，这代表着稳定最快的方式，就是沿着梯度方向下降：

7、

8、其中是约束函数c对位置p的梯度，即δp的矫正应该为约束的梯度方向，对应到每个粒子的位置矫正如式(4)所示：

9、

10、这个条件隐式地保证了线性和角动量守恒，同时允许求解未确定的约束系统。将式(3)带入式(4)得到整合后的式(5)：

11、

12、回归到多个粒子系统的约束中，可以定义缩放因子s，如式(6)所示：

13、

14、对于粒子质量并不一样的情况，则需要将式(4)乘以一个矩阵质量的系数，即式(4)被改写为其中m＝diag(m0,m1……mn-1)，mi代表粒子i的质量，同时记

15、经过推导，对于粒子i的位置修正即如式(7)所示：

16、

17、在以上概念定义的基础上，接下来对本发明的主要内容进行详细介绍，具体包括以下步骤：

18、步骤一：设计基于松弛因子的非渗透约束

19、雪粒子具有一定的刚体特性，两个刚体粒子不应该过度共享同样的空间。因此非渗透约束作为几何约束能够有效防止刚体粒子发生不合理的相互穿透和重叠。

20、另外，雪整体具有一定的不透性和一定的连续性，即使由离散粒子构成也不应呈现出过于离散和分裂的状态。非渗透约束能够控制雪粒子之间的最小距离，使得相邻粒子保持合理的间距从而增强雪的整体连续感，有助于模拟出更加自然逼真的雪效果。

21、如式(8)所示，通过确保任意两个粒子之间的距离不小于两个粒子半径之和来模拟非渗透约束。通过这样的约束可以维持雪的整体形态，避免出现粒子过度的覆盖导致不自然的仿真结果。

22、c(p)＝||pi-pj||-(ri+rj)≥0           (8)

23、其中，ri和rj代表粒子半径，所以结合式(7)得到的最终基于非渗透计算的位置矫正方程如式(9)所示：

24、

25、式(9)更深层次的物理意义实际上模拟了两个刚体球形粒子在运动过程中发生碰撞接触的情况，它基于粒子半径判断重叠，首先检查任意两个粒子i和j之间的距离是否小于两个粒子的半径之和。如果小于该距离就认为两个粒子发生了重叠，即需要进行位置修正。也就是说，对于发生重叠的两个粒子i和j，非渗透约束会计算出一个修正向量δpi，使得两粒子被分离到刚好不重叠的最小距离。修正向量的计算方法是：先求出两粒子位置连线的单位向量n，然后让i粒子沿n方向移动修正量的一半，j粒子沿相反方向-n移动另一半修正量。

26、非渗透约束提供了一种简单高效的方式，可以基于粒子半径计算出所需的位置修正量，从而驱动模拟过程中粒子位置的更新。这种约束解决方案避免了复杂的碰撞检测和处理，可以高效地模拟大量粒子之间的相互作用。

27、macklin等人提出的xpbd框架主要解决了原始pbd的刚度系数的迭代次数依赖问题，其中的λ求解方程如式(10)所示：

28、

29、可以看到xpbd的核心在于融入了之前的缩放因子s当中新的系数(称为柔度系数)，而不是使用刚度系数k。显而易见的是，当的时候式(10)就退化回了pbd框架。(参考文件10：macklin m,müller m,chentanez n.xpbd:position-based simulation ofcompliant constrained dynamics[a].proceedings of the 9th internationalconference on motion in games[c].acm,2016:49-54.)

30、从数学角度看待，实际上是刚度系数矩阵的倒数也就是矩阵的逆，它和刚度系数一样代表约束的强弱，即如果柔度系数为0，则约束必须被满足；而柔度系数越大则代表这个约束强度较弱。由此可以理解到在pbd仿真框架里有一种直观的弊端，即约束是通过直接修改粒子位置来满足的，也就导致了系统刚性很大、约束响应不够真实自然。图4形象地解释了在位置矫正的过程中，会因为刚度较强的修正导致粒子位置的突变，从而在视觉上呈现出从b到d的突然抖动或突然不合理断裂的情况。

31、考虑到在xpbd里引入会增加求解系统的复杂性，本发明提出一个轻量级的概念称为松弛因子αnon-penetration，αnon-penetration以显示的形态存在于式(5)中，主要作用于非渗透约束求解λ时的方程，如式(11)所示：

32、

33、松弛因子方程主要是将原pbd约束投影公式替换为引入了αnon-penetration的新公式以削弱约束修正量。在迭代求解过程中，αnon-penetration的存在使得系统状态在严格满足约束条件和放松约束之间取得平衡，表现出“亚刚性”的特点。也就是说松弛因子允许粒子被不严格的纠正。引入松弛因子相当于在每个约束方向上施加了一个“阻尼”，使得迭代过程更加平稳而不易出现振荡或发散。

34、步骤二：设计融合约束的基于位置动力学的迭代框架

35、除了步骤一设计的基于松弛因子的非渗透约束，本发明融合了摩擦力约束和密度约束，作为辅助约束方程一同矫正粒子位置。

36、首先，考虑到雪材质呈现疏松多孔的特征，合理的密度约束是确保雪总体质地的关键。雪由无数小的冰晶颗粒组成，这些颗粒在一定程度上是松散的，但又不能过于分散。

37、因此本发明需要融合密度约束来维持雪的整体密实度。具体可以设置每个粒子周围的密度上下限，当密度过高时向外移动粒子位置，密度过低时向内移动粒子位置，使粒子密度保持在一个合理范围内。

38、考虑到雪具有一定的流动性，而基于位置的流体模拟(position based fluids,pbf)是通过在流体粒子上应用密度约束来克服需要很小的时间步长来防止过冲的问题。在pbf中，粒子i的密度ρi的计算公式如式(12)所示：

39、ρi＝∑jmjw(xij,h)            (12)

40、其中j是粒子i的邻居，xij是粒子i和j之间的距离。w是一个核函数，用于计算邻居粒子对当前粒子密度的贡献。考虑到雪的不可压缩性，现给出密度约束方程如式(13)所示：

41、

42、这个不等式指出，当粒子i的密度ρi大于静止密度ρ0时约束函数的值为正，需要通过调整粒子位置来使密度降低。通过融合这种密度约束，能够防止粒子聚集过密从而保持模拟的稳定性。联立式(7)和式(12)，得式(14)为最终基于密度计算的位置矫正公式：

43、

44、其次，考虑到雪具有一定的内部摩擦力，因此运动时会呈现出特有的阻尼效应。而传统框架里无法对静摩擦进行建模，由于摩擦阻尼力无法校正约束求解过程中已经发生的位置变化，所以粒子位置的变化会导致模型体迅速坍塌的情况。因此本发明还融合摩擦力模型，作为计算非渗约束之后的摩擦阻尼矫正。当非渗透约束计算完毕，粒子在切线方向上的位移计算如式(15)所示：

45、

46、库仑摩擦模型是描述两个接触表面间摩擦行为的经典模型。它的基本原理讲述了静止时的最大静摩擦力与施加在表面间的法向力成正比，比例系数为静摩擦系数μs；一旦开始相对运动，动摩擦力与速度和接触面面积无关，只与法向压力和接触面粗糙程度有关，但比例系数降低为动摩擦系数μk(μk<μs)。该模型体现了静止物体受到的最大静摩擦力，以及运动过程中动摩擦力与速度无关的特性。在式(15)的基础上，将连续的摩擦力转化为对粒子位置的位移矫正得到式(16)：

47、

48、当两粒子i和j的相对法向位移|δp⊥|小于最大静摩擦力n＝μsd时，只消除所有的切向移动来进行位置矫正；否则施加动摩擦矫正，其大小由动摩擦力决定且被限制在一个最大值以避免过度矫正。

49、融合摩擦力矫正也是模拟雪很重要的一个环节，它可以在每个时间步引入基于相对速度和接触面积的摩擦力矫正，对粒子的速度进行修正。合理的摩擦力矫正能使雪在运动时逐渐减速最终停止，并能维持雪堆的形态不会发生剧烈扩散。

50、基于以上，本发明的整体迭代算法框架如算法1所示：

51、

52、算法1的循环里首先进行粒子的初始化流程，赋予每个粒子初始位置和速度以及其他参数；4-6步与原始框架基本相同，进行位置预测；8-11步主要根据约束求解位置校正，分别包括前文给出的密度约束、非渗透约束以及摩擦力约束；最后根据得到的位置矫正进行最终的速度、位置更新。

53、为了得到更加真实自然的物理模拟效果，在基于pbd的算法的迭代计算框架中，每个时间步长内pbd算法需要求解由各种约束方程组成的一个大型线性方程组以确定粒子在下一时间步的新位置。

54、这个线性方程组往往是一个巨大的系统，包含成很多个未知数，而求解如此大规模的线性方程组是一个计算量很大的问题。因此，需要使用迭代算法来近似求解而不是直接采用高斯消元法等直接求解方法。也就是对应了算法1的第8-11步，应在每一帧计算时额外重复计算多次以确保各种约束条件得到充分满足。同时，在求解约束方程时所采用的迭代方法对求解的性能和精度也有很大影响。

55、gauss-seidel(高斯-赛戴尔)迭代法和jacobi(雅可比)迭代法就是两种常用的迭代求解线性方程组的格式。它们的思路是通过不断迭代修正未知数值，使得最终的解越来越逼近真实解，这两种迭代格式各有优缺点。传统的pbd方法用gauss-seidel迭代的迭代格式，如式(17)所示：

56、

57、gauss-seidel的基本思想是：在计算某个未知量的新估计值时，利用其他未知量的已计算出的新估计值，而不是使用它们在前一次迭代中的旧估计值。也就是说，gauss-seidel算法认为新计算出来的估计值会比上一次迭代的值更加精确，因此迭代过程中会尽量利用新算出的值进行更新和收敛。这种做法的优点是收敛较快，但缺点是必须等到所有其余分量都完成计算后才能开始。对并行计算有一定限制。

58、另一种常用的迭代格式是jacobi的迭代格式，如式(18)所示：

59、

60、与gauss-seidel迭代不同，jacobi迭代法在计算时,只使用前一次迭代的结果，并没有对新计算出分量的重复多次使用，虽然这种做法会使收敛速度比gauss-seidel迭代慢一些，但是各分量的计算是相互独立的，更加适合并行计算框架。

61、总的来说，gauss-seidel迭代法的收敛速度更快但并行效率较低；而jacobi迭代法虽然收敛慢一些但完全无需同步操作，并行效率很高。考虑到并行计算的优势，本发明采用jacobi迭代法求解约束方程，计算雪粒子的位置矫正。

62、步骤三：融合希尔伯特曲线的哈希编码方式

63、在前两步已经得到融合约束的基于位置动力学的迭代框架的基础上，本发明的第三步通过在邻域搜索的步骤中融合希尔伯特曲线的方式作为哈希编码来提高整个模拟的运行帧率。

64、在模拟仿真粒子行为过程中，识别每个粒子周围的邻域粒子是一项基本且关键的任务。最简单直接的方法是采用全局搜索算法，也称为穷举搜索算法。该算法的基本思想非常简单易于理解和实现。具体来说，对于待处理的当前粒子，全局搜索算法需要遍历整个粒子系统中的所有其他粒子。对于每一对粒子，计算它们之间的欧几里得距离并与给定的邻域半径值进行比较。如果两粒子间距离小于等于邻域半径，则将该粒子确定为当前粒子的邻域粒子并存储。重复上述步骤，直至遍历完所有粒子对。

65、虽然全局搜索算法思路简单直白，但其算法复杂度为o(n2)。这意味着，当粒子数量逐渐增大时算法所需的计算时间将呈现平方级的快速增长。这种计算开销的急剧增加，使得全局搜索算法在处理大规模粒子系统时，执行效率将严重下降，难以满足实时性仿真的需求。因此，在粒子数量较大的情况下有必要探究更高效的邻域搜索算法。

66、希尔伯特曲线(hilbert curve)是以德国数学家大卫·希尔伯特命名的一种连续空间填充曲线。其空间填充特性和良好的局部性，使其在数学和计算机图形学中等许多领域都有重要的应用。它是一种分形曲线具有自相似性，即曲线的任意部分与整体结构相似。

67、可以将其理解为一个映射hn，曲线的构造始于一条连接正方形四个顶点的曲线，起点和终点是正方形的左下角和右下角:hn＝(0,0),hn＝(1,0)。然后将正方形划分为四个相等的子正方形，并在每个子正方形中递归地绘制相同的曲线，但方向可能会有所调整。随着递归的进行，曲线会逐渐填满整个正方形。如图5所示是二维平面上希尔伯特曲线的构造过程。

68、值得一提的是，hn满足一个重要的过程，即把每个的小区间映射到4个区域内。也就是说，把映射到左下角把映射到左上角等。并且，这个特点可以在迭代中保持。例如，被映射到被映射到等。

69、类似地，希尔伯特曲线也可以扩展到三维或更高维度的空间。希尔伯特曲线可以无限递归下去，随着递归的深入曲线会逐渐填满整个二维平面或者更高维度的空间。对于三维情况的曲线绘制如图6所示。

70、具体来说，图6中的希尔伯特曲线是一条连续不交叉的分形曲线，能够以自相似的方式覆盖整个三维空间。通过一系列的迭代过程，曲线被细分为越来越多的线段，最终密集地填充了整个目标空间域。在这个过程中，目标空间被划分为大量的离散单元格，每个单元格将被希尔伯特曲线一次且仅一次遍历经过。

71、希尔伯特曲线最显著的特点是能够将多维空间映射到一维，同时保持原有的空间局部性，这一独特的性质使希尔伯特曲线在多维数据的线性化编码、空间细分等领域有着广泛的应用。由于希尔伯特曲线是一条连续曲线，所以可以将每个离散单元格的位置转化为曲线上的一个唯一线性序号，即该单元格的希尔伯特编码。

72、对于一个8*8的二维网格，希尔伯特曲线编码结果如图7所示。可以看到，由于希尔伯特曲线在映射过程中能够很好地保持原始空间的局部性，即在希尔伯特曲线上相邻的点，它们在原始空间中也一定是相邻的。这种编码具有保留原空间中局部性质的特点，在高维场景同样适用：在空间域中位置相邻的单元格，其对应的希尔伯特编码值也必定相近。利用这一性质，可以将原本任意维度的空间映射为一个线性序列，而不会打乱元素间的空间局部拓扑结构。算法2给出了计算希尔伯特曲线编号的算法：

73、

74、算法2首先计算网格维度的最大值，并通过对数运算得到该值所需的最大比特位数hilbert_order。其次计算网格维度的最小值，并通过对数运算得到所需的最小比特位数，从而推导出用于邻域编码的比特位数neighbor_bits。在主循环中，按比特级别从高到低遍历：

75、(1)根据当前比特级别，从gridpos的xyz坐标中分别提取出该级别的比特值，组合成一个3比特子索引subindex。

76、(2)将subindex对应的值(基于希尔伯特曲线的映射规则)左移到相应的位置，叠加到hilbertindex上。

77、(3)根据subindex的值，对gridpos的xyz坐标进行旋转、翻转等变换，为下一级别做准备。

78、(4)将currentlevel右移一位进入下一循环层级。

79、完成循环后将hilbertindex向右移动neighbor_bits位，得到最终的希尔伯特曲线编码值hilbertidentifier。

80、本发明将采用基于希尔伯特曲线计算的网格id作为邻居网格的编号进行邻域搜索，后通过与传统的哈希编码方式对比来探究模拟帧率的提升效果。

81、总体来说，本发明的主要贡献是设计了融合约束的基于位置动力学的粒子雪模拟方法，属于计算机图形学和虚拟现实仿真领域比较重点的研究课题。本发明从约束方程的本质出发，考虑到雪的内在物理特性，设计了基于松弛因子的非渗透约束来使得雪粒子平稳地还原刚性特点，融合了密度约束和摩擦约束使得雪粒子还原了其具有一定松散度和粗糙度的特点，即不显式计算物理力的情况下还原了雪的物理特点。

82、除此之外，还通过融合希尔伯特曲线的哈希编码方式提高邻域搜索的效率，这也是本发明的重要贡献。实验结果证明，本发明所提出的模拟方案可以真实地还原雪在运动过程中的物理规律，并且在大规模的粒子数量下仍然可以满足实时的效果。