一种基于双向路径跟踪的实时声音传播模拟方法
【专利摘要】本发明公开了基于双向路径跟踪的实时声音传播模拟方法，该方法首先输入当前帧的声学环境状态，以各声源为起点随机生成正向子路径，以听者为起点生成反向子路径，在正向与反向子路径节点间随机建立连接，生成完整路径；并计算出每一路径的强度，生成概率和多重重要性采样权重，最后得到并输出能量响应曲线。本发明以双向路径跟踪的路径生成能力为基础，结合基于信噪比的质量标准，通过控制不同反射次数路径的数量对声能曲线各时间段质量做出调整，从而可以平衡各时段反射声的质量。本发明由于使用了双向路径跟踪，可以在相同的计算代价下产生数量更多，质量更高的路径，且路径质量较少受声源与听者位置的影响。
【专利说明】
一种基于双向路径跟踪的实时声音传播模拟方法
技术领域
[0001] 本发明涉及声学模拟技术领域，尤其涉及一种基于路径跟踪的几何声学模拟方 法。
【背景技术】
[0002] 早期的声学模拟技术主要应用于建筑声学设计，大多使用离线算法。当前虚拟现 实技术的发展使得实时声学模拟成为必要。当下的声学模拟算法可分为波动声学方法和几 何声学方法两大类。波动法通过求解波动方程来模拟声音传播，其结果精确性高，然而计算 代价不菲，难以达到实时模拟的要求。几何声学方法则将声音传播的规则简化为与光线相 似。这类方法模拟与声音波动性相关的现象(如衍射与干涉）时会遇到困难，然而其计算代 价远低于波动法。因此在传统建筑声学模拟领域，这一方法的应用极为广泛。近年来几何声 学方法的发展可参考SAVI0JA，L.，AND SVENSS0N，U.P.2015.0verview of geometrical room acoustic modeling techniques.The Journal of the Acoustical Society of America 138,2,708-730〇
[0003] 很多实时模拟算法利用对声学环境的假设来提高计算效率。若我们假设场景中物 体，声源或听者静止不动，我们便可通过预计算声音传播的方法实时地产生音效。这类技术 的代表可参考RAGHUVANSHI，N.，AND SNYDER，J.2014.Parametric wave field coding for precomputed sound propagation.ACM Transactions on Graphics(T0G)33,4,38〇另一些 方法利用声学环境在邻近帧之间变化较小的假设来简化计算，代表方法可参考SCHISSLER， C.,MEHRA,R.,AND MAN0CHA,D.2014.High-order diffraction and diffuse reflections for interactive sound propagation in large environments.ACM Transactions on Graphics(T0G)33,4,39。此外，一些方法通过简化模型，合并声源等方法，去除声学场景中 不重要的细节来加速计算。代表性的方法可参考ANTANI，L.，AND MAN0CHA，D. 2013.Aural proxies and directionally-varying reverberation for interactive sound propagation in virtual environments.Visualization and Computer Graphics,IEEE Transactions on 19,4,567-575.
[0004] 路径跟踪方法本质上是为模拟光能传播而特化出的蒙特卡罗积分法。由于几何声 学中的声音传播规则与光线相似，这一方法在声学模拟领域也得到广泛应用。传统的路径 跟踪法以声源为起点生成大量路径。当路径与听者相交，则将这一路径的能量贡献累加到 最终的声能曲线之上。由于代表听者的几何体(通常为人头大小的球体)相对场景而言过 小，这一方法的收敛速度极为缓慢。
[0005] 由于在实际应用中，场景中往往存在多个声源，但很少需要同时模拟多个听者， Schissler等人提出了反向路径跟踪法。具体可参考SCHISSLER，C.，MEHRA，R.，AND MAN0CHA,D.2014.High-〇rder diffraction and diffuse reflections for interactive sound propagation in large environments.ACM Transactions on Graphics(TOG)33, 4,39。该方法以听者为起点生成路径，而用几何体表示声源，声源越多，则几何体总表面积 越大，从而提高了路径的利用率。
[0006] 反向路径跟踪法的问题在于其效率与声源几何体的大小相关，且不能模拟理想点 声源。尽管我们可以缩小几何体来近似点声源，然而这也会降低路径的命中率。SchrSder提 出的路径雨（diffuse rain，DR)技术很好地弥补了这些问题。该方法不使用几何体检测路 径，而将路径上的节点直接与目标相连(对于正向路径跟踪，目标为听者，反向路径跟踪则 为声源）。具体方法可参考SCHR〇DEl?，D. 2011 .Physically based real-time auralization of interactive virtual environments, vol .11.Logos Verlag Berlin GmbH。对于完全 无遮挡的场景，这种方法可以利用到路径上的每一个节点，极大地提高了计算效率。然而， 这一技术的缺点也非常明显。由于路径雨建立的连接有一端为目标，目标的位置对连接质 量的影响是非常明显的。当目标处于不合适的位置时，这一方法的结果质量将会出现明显 的下降。在目标经常移动的场景中，为了获得质量稳定的结果，用户需要不断调节模拟程序 的参数。
[0007] 最后，上述所有基于路径跟踪的方法都有一个共同的问题，即不能对声能曲线各 时间段的质量进行精确的控制。实现这一点的障碍主要有二:其一，目前并没有一种可以对 声能曲线不同时间段质量分别评价的质量标准。其二，声音传播的延迟取决于其传播距离。 如果要调节声能曲线某一段时间的质量，我们必须控制算法产生的路径长度，而由于路径 生成的随机性，这并不是一件容易实现的事情。由于这一问题，声音传播模拟很容易出现早 期反射声或后期反射声质量过低，与其余部分难以配合的情况，对总体结果质量造成不利 的影响。
【发明内容】
[0008] 本发明目的在于针对现有技术的不足，提供了一种基于双向路径跟踪的实时声音 传播模拟方法。
[0009] 本发明的目的是通过以下技术方案来实现的：一种基于双向路径跟踪的实时声音 传播模拟方法，每一帧中的处理过程包括如下步骤：
[0010] (1)输入当前帧的声学环境状态，包括声源与听者位置、场景几何描述、声学材质 状况和声音介质特性；
[0011] (2)根据积分TnLo的米样概率xn,d，向各积分分配样本;T nLo代表直接声Lo在场景中 经过n次反射后对能量响应曲线的贡献，d表示当前帧的序号;d = l时，xn,i=l/N，N为N为积 分TnLo的个数，即最大反射次数；
[0012] (3)对于每一积分，在积分内的各连接策略间尽可能平均地分配样本；
[0013] (4)以各声源为起点随机生成正向子路径，以听者为起点生成反向子路径；
[0014] (5)根据先前确定的各连接策略的样本数，在正向与反向子路径节点间随机建立 连接，生成完整路径；
[0015] (6)计算出每一路径的强度，生成概率和多重重要性采样权重;对于通过连接正向 子路径第s个节点与反向子路径第t个节点得到的路径 XfXs+t+1，使用公式如下：
[0019] 其中f( ?)为路径强度，Ps,t( ?)为路径生成概率，Ws,t( ?)为多重重要性采样权 重；
[0020] Fi为节点Xl处的声材质特性，G1>1+1为路径中各段的几何传播因子，11 1>1+1为路径中 各段的声介质特性，I代表路径中各段的生成概率，cs,t代表在正向子路径第s个节点与反 向子路径第t个节点建立连接的概率，星号为卷积符号；
[0021] (7)使用步骤6得到的f( ? )，ps,t( ?)和ws,t( ?)计算能量响应曲线/j = i/%，TnLo 0 的计算公式为：
[0023]其中，E[ ?]表示数学期望，Sn为用于计算TnLo的样本总数，sn>1- n为连接正向子路径 第i个节点与反向子路径第n-i个节点得到的样本数，为使用这种连接方式得到的第 j个样本；
[0024] (8)计算本帧各积分在各信噪比窗口中的方差〇2,令其质量因子Q为计算积分使用 的样本数，对于每一方差估计〇f，使用下面的更新公式：
[0028] 其中，Gf为用户定义的方差质量标准，d=l时€^1与〇(1-1应取0;
[0029] (9)使用下面的公式计算出下一帧的各积分采样概率：
[0031]其中，M为信噪比窗口的数量;a为小于1的正数;为步骤8中更新后的，积分TnL〇 在信噪比窗口m中的方差估计；
[0032 ] (10)输出步骤7得到的能量响应曲线。
[0033] 本发明的有益效果是:与传统方法相比，双向路径跟踪生成路径的能力更强，路径 的平均质量也更好。本发明以双向路径跟踪的路径生成能力为基础，结合基于信噪比的质 量标准，通过控制不同反射次数路径的数量对声能曲线各时间段质量做出调整，从而可以 平衡各时段反射声的质量。本发明由于使用了双向路径跟踪，可以在相同的计算代价下产 生数量更多，质量更高的路径，且路径质量较少受声源与听者位置的影响。本发明的样本分 配方法则可以平衡声能曲线各部分的质量，使结果总体上更加令人满意。
【附图说明】
[0034] 图1是本发明中关于连接策略的示意图；
[0035] 图2是试验场景中，本发明与对照方法生成的声能曲线对比图；
[0036] 图3是试验场景中，本发明与对照方法的信噪比曲线对比图；
[0037] 图4是示例情况1的场景布置示意图；
[0038] 图5是示例情况2的场景布置示意图；
[0039] 图6是对照方法在情况1与情况2下的信噪比曲线对比图；
[0040] 图7是本发明在情况1与情况2下的信噪比曲线对比图。
【具体实施方式】
[0041]在基于几何声学原理的各种声学模拟方法中，路径跟踪由于其实现简单，场景适 应性强和最终估计无偏的特点而被广泛应用。然而之前的路径跟踪算法存在采样效率低， 结果质量差且不稳定，不能支持理想点声源等种种问题。本发明使用双向路径跟踪，即从声 源与听者两个方向同时生成部分路径，而后在其间建立连接以产生完整路径的方法解决了 上述问题。
[0042] 下面首先介绍本发明的理论基础。
[0043] -、几何声学与蒙特卡罗积分法基础
[0044] 几何声学模拟的最终结果用各声源的能量响应曲线E(x，co，t)表示。其中x代表听 者位置，w代表声音的传入方向，t代表声音传播的时间延迟。
[0045] 几何声学中，声能的传播可用下面的方程描述：
[0046] 乙（a' 4 a.，') = (Ay 4 a'，') + jL、（a' 4 x''）* M (a' e n
[0047] 其中：
[0048] Lk (x' -> x, /) = A (-> x', /) G (xf/ x') f (xf x)
[0050]这里我们解释一下式中符号的含义。我们一般用L(x，co，t)来代表随时间变化的 声辐射度，其参数与能量响应曲线E(x，co，t)的含义相同。为方便描述传播过程，我们用L (x'-x，!:)表示L(x'，〇 (x^4x)，t)。式中其余情况同理。Lo代表声源的出射福射度。f(x〃 - X' -x)为双向散射分布函数，代表X'处声学材料的性质。G -v)为几何因子，用于表不 声音的扩散与场景中的遮挡关系。其中当x与x'互相可见时为1，反之为0。星号* 代表卷积。表示 x〃与Y之间声音介质的影响。这种影响包括能量吸收和传播 延迟。对于均匀介质，我们有
[0052]其中a为介质的能量吸收系数，5( ?)为单位冲激响应函数，c为声速。最后，Q代表 场景中的几何体。A则是几何体上定义的面积测度。
[0053] 注意上面的模型不能描述衍射和干涉等波动现象。对于不同频率的声音，我们假 设它们互不相干，每一频率上都可算得一独立的声辐射度。
[0054] 我们可以将上面的积分方程写为算子形式：
[0055] L = Lo+TL
[0056] 其中
[0057] IE - ^ .y'/) ^ M ^ a-'.t)dA n
[0058]简单地说，算子T的含义相当于"一次反射"。我们可将这一方程展开为Neumann级 数：
[0059] Lz=YJ'，L^
[0060] 可以看到，最终的声辐射度等于声音每一次反射后的辐射度之和。由于算子T本质 上是一次积分，我们便可以用数值积分法求出级数中的各项，以获得最终结果。在实际应用 中，由于计算能力有限，且声能在场景中反射之后会不断衰减，我们设定一最大反射次数N， 并只考虑r^N的情况。
[0061] 接下来我们介绍蒙特卡罗积分法。设我们需计算积分'我们可以在I：上 v 构建概率空间（2，x，y2)，从中取出一样本X，然后求下面的期望：
[0063] 这里的p(X)是定义于两测度间的Radon-Nikodym导数。也经常不严格地称为"取到 样本X的概率"。对于我们的问题，p(X)的公式将在后文给出。
[0064]最后，能量响应曲线E(x，《，t)与声辐射度的L(x，《，t)形式一致。事实上，假设听 者x处有一无限小球体，则我们在其球面上有E(x, ?，t)=L(x, ?，t)。
[0065] 二、双向路径跟踪
[0066] 我们之前提到过，路径跟踪本质上是特化了的蒙特卡罗积分法。我们现在将路径 跟踪中的概念和上文中的公式一一对应。一条连接声源与听者，含有i个中间节点的路径可 以视为从Q 1中取出的一个样本，这里Q 1为i个集合Q的笛卡尔直积，称为i次反射的路径空 间。这条路径对应的f( ?)称为该路径的强度，P( ?)称为该路径的生成概率。使用蒙特卡罗 积分法时，这条路径将用于计算积分ru。
[0067] 在双向路径跟踪中，我们首先从声源和听者出发，随机生成大量的子路径(从声源 出发的子路径称为正向子路径，从听者出发则称为反向子路径）。然后在其节点间建立连接 以生成完整的路径。假设一条子路径为xo… Xs，那么这条子路径的生成概率可用如下公式计 算：
[0069] 这里Pg(XQ -XI)代表从声源或听者XQ出发，生成一条XQ -XI方向的射线的概率，Pf (Xi-l^Xi^Xi+l)代表Xi接收到一条Xi-l^Xi方向的射线时，生成一条Xi-Xi+l方向的反射射 线的概率。这些概率由具体的路径生成实现方法决定。为简单起见，我们用Vo代表p g(X04 Xl)，Vi代表Pf(Xi-1-Xi-Xi+l)，Gi,i+l代表G -Y/+| )。那么上面的公式可以简化为： S-1
[0070] 0
[0071] 现在，我们假设有一条子路径的第s个中间结点与另一条路径的第t个中间结点相 连，构成一条s+t次反射的完整路径X0…Xs+t+1。那么这条路径的生成概率为：
[0073] 然后我们考虑一下f( ?)。根据与当前路径相关的积分Ts+tLo的展开形式，我们有：
[0074] /(% …Xs+f+1) = fps (巧 4 7；-0：
[0075] 我们希望能够将上面的公式化为与路径生成概率相似的形式。于是，我们定义Fi = f(Xi-1-Xi-Xi+i)，M,_ (,y, G 丨）。特别地，我们令卩。=1。(叉。4叉1)，卩糾+1=1，于是有： S+/+I 州
[0076] /(~…~+1)=rT^r^w *??? ) i~0
[0077] 由于冲激响应间的卷积结果仍然是冲激响应，我们只需将它们的强度相乘，时间 延迟相加就可以得到结果。
[0078] 现在我们考虑如何计算积分。由于我们这里使用了多重重要性采样技术，我们的 积分计算公式比一般的蒙特卡罗积分法要复杂一些：
[0080]我们解释一下公式中符号的含义。我们知道，构成一条反射次数为n的路径有很多 种方法。只要连接正向子路径的第i个中间节点和反向子路径的第n_i个中间节点即可。公 式中Xy-u代表采用这种连接方式的第j个样本，m.M代表这种连接方式的实际样本数，Sn n 代表计算积分TnLo所用的总样本数，即乏代表使用这种连接方式的概率，可以是 r-0 Siti/Sn，也可以是分配样本所用的概率分布。Wwi为多重重要性采样权重。本算法中使用 的权重如下：
[0082] 这一权重称为平衡型权重。原理参见VEACH，E.，AND⑶IBAS，L.J.1995.0ptimally combining sampling techniques for monte carlo rendering.In Proceedings of the 22nd annual conference on Computer graphics and interactive techniques，ACM， 419-428。
[0083]本文中，i、j、k均为数学表示，无特殊意义，取值为自然数。
[0084]三、样本分配策略
[0085]在讨论样本分配策略之前。我们需要先介绍一下我们基于信噪比，用于评价声音 传播模拟算法的质量标准。我们将能量响应曲线按时间分为多个窗口，每一个窗口内的能 量积分视为一个随机变量。通过多次模拟的结果，我们可以算得各随机变量的均值与标准 差估计。我们定义随机变量X的信噪比如下：
[0087]各随机变量的信噪比按时间排列可组成一"信噪比曲线"。这一曲线反映了能量响 应曲线各部分的质量。我们后面的样本分配将以这一标准为优化目标。
[0088] 样本的总数应在算法开始执行前确定，每一次模拟之前的样本分配过程分为两 步:第一步在各个积分TnLo中分配样本，第二步在属于同一积分的各种连接策略(决定使用 前向子路径的第几个中间节点）中分配样本。产生样本时，我们随机取两条子路径，并按先 前确定的连接策略确定连接，计算积分。
[0089] 附图1给出了子路径最大反射次数为6的情况下各种连接策略的示意图。注意各种 策略能够提供的连接数是不相同的，其原因在于不同的子路径可能包含相同的起点。对于 直接声来说，只有一种可能的连接，就是将声源与听者直接相连。涉及子路径起点的连接策 略比不涉及起点的连接策略能够提供的连接数要少一些。标为"不建议使用"的连接策略容 易对最终结果造成不良影响，应当略过。
[0090] 对于同一积分的不同策略（在图1中从右上到左下排成一斜线），我们可以证明最 优的样本分配方案为平均分配，然而由于各策略最大连接数不同，有可能出现不能平均分 配的情况，此时应使实际的样本分配与平均分配尽可能接近。
[0091] 各积分之间的样本分配以最大化总信噪比之和为目标。假设待分配的总样本数为 S，能量响应曲线划分的窗口数为M，积分个数为N，积分rLo的采样概率为 Xl，那么总信噪比 之和为
[0093] 这里Xmn为积分TnL〇对窗口m的贡献。E |；X_与采样方式无关。然而；fXw贝1J _ n=l J \_n?=\. _ 是可以通过调整样本分布来控制的。设)(_互不相关，则有
[0094] O'2 YXmn
[0095] 其中是积分TnLo中单个样本对窗口 m贡献的方差。由上可知，为了最大化总信噪 比，我们需要求解下面的问题： .M N.. N.
[0096] min^5lg<r^^-~, s.t. ^.v,t =1. 1 m~\ n-l H-1
[0097] 为了在积分之间进行样本分配，我们需要做两件事情。一是求出二是解上面 的最优化问题。
[0098] 我们首先讨论〇i。在每一次模拟结束后，我们当然可以用模拟得出的结果求出 然而实时模拟使用的样本数一般不会太多，容易导致方差估计不稳定。为了保证方差 估计相对稳定，我们将用当前帧的新方差估计与之前的估计做加权平均。
[0099]设我们在当前帧使用Q个样本作出了新的方差估计〇2,我们将Q记为这次估计的 "质量因子"。设前一帧的最终方差估计为，质量因子为Qh，则当前帧的方差估计和质 量因子按下面的公式更新：
[0102]其中Cf为一预先制定的质量标准。y的计算公式如下：
[0104]注意对于Qd-iiQ#的特殊情况(实际中极为常见），上面的公式可以简化为
[0106]现在我们讨论前面的优化问题。为求解这一问题，我们给出了下面的迭代公式：
[0108] 为防止除数为0,这一迭代公式要求对于任意M和N至少有一个不为0。我们可以 略去对最终能量响应曲线贡献太小的积分和窗口来达到这一要求。另一种方法是设定一个 极小的方差估计值下限，当估计值过低时将其强制提高至下限以上。我们可以证明，当a取 小于1的正数时，上面的迭代公式局部线性收敛至优化问题的解，且收敛速度至少为线性。 实际应用中，a的取值范围在0.4到0.7之间效果较好。
[0109] 每一帧结束前，我们首先使用当前帧的数据更新方差估计。然后使用当前帧的各 积分采样概率作为初值带入迭代公式，计算出下一帧的概率。当场景静止时，这样的迭代将 收敛至最优概率分布。当场景出现变化时，由于变化的连续性，之前的概率分布仍然是一个 很不错的初始值。
[0110]最后需要注意的是，我们需要给采样概率为〇的积分稍微保留一点样本，以保证下 一帧的所有方差估计都可以计算出来。
[0111]四、方法实现
[0112] 模拟开始前，我们需设定一初始的各积分采样概率Xn(建议平均分配）。所有的方 差估计与质量因子Q应初始化为〇。下面给出本方法中每一帧的实际流程：
[0113] 1、首先根据各积分采样概率Xn，向各积分分配样本。注意对采样概率为0的积分需 保留少量样本用于本帧的方差估计；
[0114] 2、对于每一积分，在积分内的各策略间尽可能平均地分配样本；
[0115] 3、随机生成子路径；
[0116] 4、根据先前确定的各策略的样本数，在子路径节点间随机建立连接，生成完整路 径；
[0117] 5、使用第二章中的方法，计算出每一路径的f( ? )，p( ?)和多重重要性采样权重；
[0118] 6、根据第二章的积分计算公式算出声能曲线；
[0119] 7、根据本帧中产生的样本，使用第三章中的公式更新方差估计(Ti与质量因子Q;
[0120] 8、根据第三章中的迭代公式，计算出下一帧的各积分采样概率。
[0121] 实施实例
[0122] 发明人在一台配备Intel i7 3.50GHz四核中央处理器，32GB内存的计算机上实现 了上文所描述的算法，并实现了用于对照的另一算法(反向光线跟踪与路径雨相结合）。实 验结果表明，我们算法的效率在大量场景中高于传统方法，在空旷的大型室内场景中，声源 离地面较近的情况下，我们的方法与传统方法的效率差异尤其明显。附图中给出了一个较 大的室内场景中，两方法在相同计算时间下产生的能量响应曲线（图2)与其信噪比（图3)。 此外，对照算法结果质量受声源位置影响的问题在我们的算法中得到了明显的缓解。附图 给出了一组示例场景中的对比。图4和图5中，光源与摄像机的位置分别对应声源与听者的 位置。图6与图7展示了两种算法在情况1(图4)与情况2(图5)下的结果信噪比曲线。可以看 到对照方法的结果有明显的质量差异，而我们的方法质量则基本稳定。
【主权项】
1. 一种基于双向路径跟踪的实时声音传播模拟方法，其特征在于，每一帧中的处理过 程包括如下步骤： (1) 输入当前帧的声学环境状态，包括声源与听者位置、场景几何描述、声学材质状况 和声音介质特性等； (2) 根据积分TnLo的米样概率Xn,d，向各积分分配样本。TnLo代表直接声Lo在场景中经过η 次反射后对能量响应曲线的贡献，d表示当前帧的序号。d = 1时，Xn, i = 1/Ν，N为N为积分TnL0 的个数，即最大反射次数。 (3) 对于每一积分，在积分内的各连接策略间尽可能平均地分配样本。 (4) 以各声源为起点随机生成正向子路径，以听者为起点生成反向子路径。 (5) 根据各连接策略的样本数，在正向与反向子路径节点间随机建立连接，生成完整路 径。 (6) 计算出每一路径的强度，生成概率和多重重要性采样权重。对于通过连接正向子路 径第s个节点与反向子路径第t个节点得到的路径 XfXs+t+1，使用公式如下：其中f( ·)为路径强度，ps,t( ·)为路径生成概率，ws,t( ·)为多重重要性采样权重。F1为节点X1处的声材质特性，G1,1+1为路径中各段的几何传播因子，M 1,1+1为路径中各段 的声介质特性,V1代表路径中各段的生成概率，cs,t代表在正向子路径第s个节点与反向子 路径第t个节点建立连接的概率，星号为卷积符号。 (7) 使用步骤6得到的f( · )，ps,t( · MPws,t( ·)计算能量响应曲自 的计 算公式为：其中，E[ ·]表示数学期望，Sn为用于计算TnLo的样本总数，sn,μ为连接正向子路径第i 个节点与反向子路径第n-i个节点得到的样本数，X1^1,」为使用这种连接方式得到的第j个 样本。 (8) 计算本帧各积分在各信噪比窗口中的方差〇2,令其质量因子Q为计算积分使用的样 本数，对于每一方差估计Of，使用下面的更新公式：其中，Qt为用户定义的方差质量标准，d = 1时€^_1与〇(1-1应取O。 (9)使用下面的公式计算出下一帧的各积分采样概率：其中，M为信噪比窗□的数量。α为小于1的正数。为步骤8中更新后的，积分TnLo在信 噪比窗口m中的方差估计。 (I 〇)输出步骤7得到的能量响应曲线。
【文档编号】G01S11/14GK105893719SQ201610436309
【公开日】2016年8月24日
【申请日】2016年6月16日
【发明人】任重, 曹春晓, 周昆
【申请人】浙江大学