GPU上KD树的射线追踪加速系统及KD树输出方法

本发明属于声学技术领域，涉及水声学快速声场强度计算的仿真技术，具体是一种GPU上KD树的射线追踪加速系统及KD树输出方法，用于对于水下目标的快速预报。

当今，准确计算水下目标在各种姿态下的声场特性，对等水下兵器的目标检测，以及国防建设有着重要的意义。随着和声呐技术的发展，要求对目标的声场强度特性计算的精度更高。目前分析水下目标回波特性的方法主要有基于Kirchhhoff 近似法的板块元方法，以及引入了Gordon积分法的改进板块元方法。而这两种方法在目标表面面元两两遮挡的情况下，会进行复杂的遮挡判断导致计算量大且效率低。

范军，汤渭霖等在2012年发表在流体与结构声学上的非专利文献“声呐目标回声特性预报的板块元方法”中提出了将计算雷达散射截面的板块元方法引入到声呐目标强度计算中的方法。在应用Kirchhoff近似计算水中目标散射声场时,用一组平面板块元近似复杂形状目标曲面,所有板块元的散射声场的和就是总散射声场的近似值。通过将单个板块元的积分化为代数和避免了面积分运算,使得板块元法比通常的面元积分法计算速度提高许多倍。针对声呐的情况将此方法推广到目标近场和非刚性表面回声特性预报。但该方法在板块元空间坐标划分时会出现积分分母为零的情况,因而由面元计算得出的目标强度可能存在奇异值,从而导致计算不稳定。

孙乃葳，李建辰等2016年发表在技术上的非专利文献“基于改进板块元法的潜艇目标强度预报仿真”中提出了一种在计算目标强度时可以避免产生奇异值的方法。该方法针对板块元法计算潜艇目标强度时积分分母可能为零,从而导致计算结果不稳定的问题,将Gordon积分算法用于潜艇目标强度的预报仿真,并针对复杂目标简化面元遮挡判断流程。因此该方法的计算结果不会出现奇异值，并且更加稳定。但是由于对遮挡面元的判断，导致该方法在计算遮挡较多的目标时，计算结果不够精确。

现有技术中，对于水下目标的预报计算精度不够高，计算时间不能满足实时性要求。

为了克服上述现有技术的不足，本发明提供了一种计算精度高，计算速度快的GPU上KD树的射线追踪加速系统及KD树输出方法。

本发明首先是一种基于GPU的KD树的射线追踪加速系统，依信号连接方向包括和连接有数据预处理模块，加速搜索生成模块，虚拟孔径面生成模块，射线追踪模块，积分模块，界面显示模块，其特征在于，在系统的最初设有显存预申请模块，所述加速搜索生成模块为KD树加速搜索生成模块，各模块相互连接共同构成基于GPU 以KD树加速搜索实现射线追踪，对各模块分述如下：

显存预申请模块：该模块的输入为三角面元模型的三角面元的数量，在该模块内部根据三角面元的数量计算出三角面元包围盒模型数组所占显存大小，在CPU端预先一次性申请GPU上的显存空间，该显存空间至少足以无冲突生成KD树，该模块的输出是预申请的显存的起始地址；

数据预处理模块：在GPU上，使用显存预申请模块申请的显存，对输入该模块的三角面元模型的每个三角面元并行构建包围盒，并初始化当前入射角度，输出包围盒模型数组；

KD树加速搜索生成模块：该模块是在GPU上并行生成的，在KD树上加速射线追踪过程，以地址访问的方式动态使用显存预申请模块申请的显存，加速KD树的生成，数据预处理模块中输出的包围盒模型数组为该模块的输入，该模块的两路同步输出，一路为对应于该包围盒模型数组的KD树，另一路为重新组织过的包围盒模型数组，该两路同步输出直接作为射线追踪模块的输入；

虚拟孔径面生成模块：该模块的输入为数据预处理模块输出的包围盒模型数组及当前入射角度，在该模块内计算包围盒模型的边界值，并在GPU上并行生成声线管束，该模块的输出为若干个相同尺寸的声线管束；

射线追踪模块：该模块的输入分别为虚拟孔径面生成模块输出的声线管束和加速搜索生成模块的两路同步输出，该三路为同步输入，在GPU上并行加速声线管束的射线追踪，射线追踪过程是通过遍历KD树完成的，该模块的输出为声线管束在遍历KD 树后的反射信息；

Gordon积分模块：Gordon积分模块接收追踪模块输出的反射信息，该模块使用Gordon积分公式在GPU上对声线管束进行并行积分，使用归约操作对所有声线管束的积分结果进行并行求和，该模块的输出为当前角度下所有声线管束积分值的代数和；

界面显示模块：Gordon积分模块输出的积分值及相关信息均显示在界面上，便于信息查看与调试。

本发明还是一种KD树叶节点内包围盒数组统一输出方法，其特征在于，包括有如下步骤：

1)输入数据并初始化：输入三角面元包围盒模型数组，以及两个辅助数组，辅助数组的长度和KD树每层的包围盒数组等长，用于辅助划分当前节点内的包围盒数组，初始化根节点作为当前节点；

2)终止剖分判断：判断当前节点是否满足终止剖分条件，当前节点地址位置中包含的包围盒数量是否小于终止剖分阈值，若满足，执行步骤7)，若不满足，执行步骤 3)；

3)计算最佳剖分面位置：用中分法或SAH法计算最佳剖分面位置和剖分轴；

4)辅助数组赋值：利用两个辅助数组辅助划分，根据最佳剖分面位置与剖分轴，位于最佳剖分面左侧的包围盒，在第一个辅助数组中将对应位置标为1；位于最佳剖分面右侧的包围盒，在第二个辅助数组中将对应位置标为1；被最佳剖分面切分的包围盒，在两个辅助数组中将对应位置均标为1；

5)扫描辅助数组：计算当前节点中的数组复制到左、右子节点时的起始存储地址，并得到数组中每个元素在左、右子节点对应地址中的位置索引；

6)剖分当前节点：将当前节点剖分为左、右两个子节点，根据位置索引，将当前节点内包围盒数组复制到左、右节点存储地址的相应位置中，更新左、右子节点相关信息，并将左、右子节点分别作为当前节点进入步骤2)，再次进入KD树构建循环判断；

7)形成输出数组：将当前节点标记为叶节点，将叶节点存储的地址内的数组附加到输出数组尾部，并将叶节点内的数组在输出数组中的地址记录在叶节点内，留给射线追踪模块使用；

8)终止剖分过程：当所有的节点都终止剖分之后，得到一棵构建完成的高质量 KD树，以及重新组织后的包围盒模型数组。

本发明解决了现有技术中对于表面具有较多遮挡的凹面目标，被遮挡面元需要计算复杂的遮挡关系，进而导致计算量大量增加且计算效率较低，为水下复杂目标声场强度的快速准确计算提供了技术手段。

与现有技术相比，本发明的技术优势：

计算精度高：相对于传统板块元及改进板块元计算方对被遮挡面元采取抛弃策略而导致的精度下降，本发明对于表面具有较多遮挡的凹面目标，使用射线追踪技术获得被遮挡面元反射的声波，因此使用本发明的系统可以获得较高的计算精度；

KD树构建速度快：现有技术中只将KD树构建过程的SAH计算部分使用GPU加速，并没有做到完全的GPU加速，本发明将父节点剖分成子节点的过程也使用GPU 进行加速，减少了GPU上KD树构建的时间，本发明在构建KD树之前一次性申请了足够无冲突构建KD树的显存空间，以地址访问的形式代替了建树过程中较为耗时的动态申请，因此本发明中的建树方法可以加快KD树的构建速度；

KD树结构所占显存小：由于KD树节点内部不保存包围盒数组，只保存节点内包含包围盒数组在显存中的起始地址，在使用时根据地址访问，因此，本发明构建的KD 树所占显存较小。

图1为本发明中基于GPU的KD树的射线追踪加速系统的结构示意图；

图2为一种KD树叶节点内包围盒数组统一输出方法流程图；

图3KD树加速搜索生成模块内部构成示意图；

图4为当前节点内的包围盒数组划分至左、右子节点并终止剖分示意图；

图5为实施例7中使用的三角面元模型示意图；

图6为使用本发明对图5所示模型仿真的声场强度曲线图；

图7为使用本系统对图5所示模型仿真的射线追踪模块计算用时及生成的声线管束数量曲线图。

下面结合附图对本发明详细说明

实施例1：

射线追踪技术被广泛应用于图片渲染，电磁场强度计算等领域，射线追踪技术还可应用于对水下复杂目标声场强度的快速准确预报，在传统的水下目标声场预报仿真系统中，使用的都是板块元及改进板块元方法对水下目标进行仿真预报。传统技术对于表面遮挡较多的目标，需要进行大量的面元遮挡判断，并将被遮挡的面元抛弃，从而导致计算精度有一定的损失。随着GPU技术的不断发展，在GPU上应用一些更为快速的算法成为可能。本发明针对这种现状展开了研究，提出一种基于GPU的KD树的射线追踪加速系统，依信号连接方向包括和连接有数据预处理模块，加速搜索生成模块，虚拟孔径面生成模块，射线追踪模块，积分模块，界面显示模块，参见图1，本发明在数据预处理模块之前设有显存预申请模块，本发明中加速搜索生成模块为KD 树加速搜索生成模块，系统各模块相互连接共同构成基于GPU以KD树加速搜索实现射线追踪，对各模块分述如下：

显存预申请模块：该模块的输入为三角面元模型的三角面元的数量，在该模块内部根据三角面元的数量计算出三角面元包围盒模型数组所占显存大小，在CPU端预先一次性申请GPU上的显存空间，该显存空间至少足以无冲突生成KD树，该模块的输出是预申请的显存的起始地址。

数据预处理模块：在GPU上，使用显存预申请模块申请的显存，对输入该模块的三角面元模型的每个三角面元并行构建包围盒，并初始化当前入射角度，输出包围盒模型数组。

KD树加速搜索生成模块：本发明中该模块是在GPU上并行生成的，在KD树上加速射线追踪过程，以地址访问的方式动态使用显存预申请模块申请的显存，加速KD 树的生成，数据预处理模块中输出的包围盒模型数组为该模块的输入，该模块的两路同步输出，一路为对应于该包围盒模型数组的KD树，另一路为重新组织过的包围盒模型数组，该两路同步输出直接作为射线追踪模块的输入。

虚拟孔径面生成模块：该模块的输入为数据预处理模块输出的包围盒模型数组及当前入射角度，在该模块内计算包围盒模型的边界值，并在GPU上并行生成声线管束，该模块的输出为若干个相同尺寸的声线管束。

射线追踪模块：本发明中追踪模块为射线追踪模块，该模块的输入分别为虚拟孔径面生成模块输出的声线管束和加速搜索生成模块的两路同步输出，该三路为同步输入，该模块在GPU上并行加速声线管束的射线追踪，射线追踪过程是通过使用虚拟孔径面模块生成的声线管束遍历KD树加速搜索生成模块生成的KD树完成的，该模块的输出为声线管束在遍历KD树模型后的反射信息。

Gordon积分模块：Gordon积分模块接收追踪模块输出的反射信息，该模块使用Gordon积分公式在GPU上对声线管束进行并行积分，使用归约操作对所有声线管束的积分结果进行并行求和，该模块的输出为当前角度下所有声线管束积分值的代数和。

界面显示模块：Gordon积分模块输出的积分值及相关信息均显示在界面上，便于信息查看与调试。

本发明改变入射角度，重复调用虚拟孔径面生成模块、射线追踪模块、积分模块、界面显示模块，完成对三角面元模型0-180度的积分计算及显示。

由于在水下目标声场强度仿真预报系统中，对水下目标声场强度仿真预报的速度及精度有较高的要求，本发明使用射线追踪技术提高了传统方法的预报精度，使用KD 树结构加速了射线追踪的搜索过程，同时将KD树的构建以及射线追踪过程使用GPU 进行加速，进一步提高了射线追踪的速度，从而提高了对水下目标声场强度仿真预报的速度与精度。在KD树的构建中，将叶节点内的数组统一组织到输出数组内，既方便了射线追踪过程中对数据的合并访问，也节省了KD树的显存占用。

实施例2：

基于GPU的KD树的射线追踪加速系统同实施例1，本发明中显存预申请模块预先申请的显存空间大小为 (D·P·I)+2F+(P·I)

其中，D为预估的最大树深，P为膨胀系数，I为根节点中包围盒数组所占显存的大小，F为一个辅助数组所占显存的大小，该模块为KD树的每层都分配了大小为 (P·I)的显存空间，同时为输出数组申请了大小为(P·I)的显存空间。本发明在剖分父节点时，根据剖分面及剖分轴，将父节点的包围盒模型数组拷贝至属于子节点的显存块中，因此在并行拷贝的过程中，不会发生显存访问冲突；此外，膨胀系数的设定也保证了被剖分面分割的三角面元包围盒可以被安全的放置到两个子节点中。

在现有技术中，KD树的构建过程中是动态申请显存空间的，此过程较为耗时。本发明的显存预申请模块在构建KD树之前一次性申请了足够无冲突构建KD树的显存空间，以地址访问的形式代替了建树过程中较为耗时的动态申请，因此本发明可以加快 KD树的构建速度。

实施例3：

基于GPU的KD树的射线追踪加速系统同实施例1-2，

在现有技术中，在KD树构建完成后，在KD树的叶节点内存储了数据的索引，由于SAH算法在节点内包围盒数量大于剖分阈值时会提前停止剖分，因此节点内部数据的索引需要的内存是不确定的。本发明中，只在叶节点内保存数组在输出数组中的地址段信息，同时各叶节点内的数组在输出数组中是连续的，在射线追踪过程中可以合并访问数据，加快射线追踪的过程。

本发明KD树加速搜索生成模块生成KD树的过程是在GPU上并行完成的，本发明KD树加速搜索生成模块利用KD树叶节点内包围盒数组统一输出方法生成KD树，参见图3，对KD树加速搜索生成模块中设有的子模块分述如下：

条件判别子模块：该子模块首先接收数据预处理模块输出的数组，根据元素数量判断决定当前节点数据是进入节点剖分子模块还是终止剖分子模块，如果进入终止剖分子模块，则停止当前节点的剖分。如果进入节点剖分子模块，则之后依次进入子节点内包围盒数组生成子模块、子节点信息更新子模块，最后再返回到条件判别子模块对子节点进行再次判断。之后循环接收节点剖分子模块反馈的数据，反复循环判断，逐层构建KD树结构，直至所有节点都进入终止剖分子模块。

节点剖分面计算子模块：该子模块的输入为条件判别子模块的输出，该子模块根据接收到的数据计算当前节点的剖分面位置和剖分轴，将一个父节点剖分成两个子节点，子节点中包含有父节点的部分信息，并将计算结果和该子模块接收到的数据输出到子节点内包围盒数组生成子模块。

大节点剖分单元：根据节点剖分子模块的输入使用中分法计算当前节点的剖分面位置及剖分轴。

小节点剖分单元：根据节点剖分子模块的使用SAH法计算当前节点的剖分面位置及剖分轴。

子节点内包围盒数组生成子模块：该子模块的输入为节点剖分面计算子模块的输出，根据输入，在GPU中将当前节点内包围盒数组划分至当前节点的两个子节点中，并行输出生成的两个子节点分别至两个子节点信息更新子模块。

子节点信息更新子模块：该子模块为两个结构相同并列的子模块，各自的输入为子节点内包围盒数组生成子模块输出的一个子数组，当前节点的子节点先继承当前节点的线索值，然后根据接收到的子数组更新子节点的相关信息，并将子节点作为当前节点与子数组输出至条件判别子模块进行重新判断；再次进入KD树构建循环过程。

剖分终止子模块：该子模块的输入为条件判别子模块输出的数据，该子模块终止当前节点内包围盒数组的剖分，并将当前节点包含的子数组输出至输出数组中，等待全部节点停止剖分后一次性输出。

本发明中KD树加速搜索生成模块的输入为数据预处理模块输出的包围盒模型数组，KD树加速搜索生成模块的输入首先进入该模块内部的条件判别子模块，在条件判别子模块内判断将条件判别子模块的输入输出到节点剖分子模块还是剖分终止子模块；若条件判别子模块的输出进入剖分终止子模块，则停止当前节点的剖分；若条件判别子模块的输出进入节点剖分子模块，则先在节点剖分子模块中进行计算，再将节点剖分子模块的输出作为子节点内包围盒数组生成子模块的输入，子节点内包围盒数组生成子模块在GPU上加速生成子数组的过程，并输出两个子数组分别作为条件判别子模块的单独输入，并行进行节点剖分，直到进入剖分终止子模块。

实施例4：

基于GPU的KD树的射线追踪加速系统同实施例1-3，

现有技术中，在KD树父节点剖分成子节点的过程中，需要先在父节点内部数组在剖分轴的方向上进行排序，排序后再根据最佳剖分面位置将数组分别复制到左、右子节点中。由于排序算法的时间复杂度较高，会影响到整体的构建时间。为缩短KD 树构建时间，本发明没有对父节点内部的数组进行排序，而是利用GPU的并行加速特性，根据最佳剖分面位置，仅判断父节点内的数组是复制到左子节点？还是复制到右子节点？或者同时复制到左、右两个子节点中。大大缩短了数据移动的时间。

所述子节点内包围盒数组生成子模块是在GPU上执行的，对子节点内包围盒数组生成子模块内设有的单元分述如下：

扫描单元：该单元的输入为节点剖分子模块的输出，在GPU上使用包含扫描操作，并行计算出两个子节点内包围盒数组元素在当前节点内包围盒数组中的索引，并将索引保存在显存预申请模块申请的辅助数组对应的位置上并输出给数组索引单元。

数组索引单元：该单元的输入为扫描单元输出的索引，在GPU上并行索引，将当前节点存储的地址中的数组按照索引分别复制到两个子节点存储的地址中，将当前节点存储地址内包含的数组划分至两个子节点存储的地址内，然后将两个子节点分别输出至条件判别子模块进行重新判断。

现有技术中只将KD树构建过程的SAH计算部分使用GPU加速，并没有对数组从父节点复制到子节点的过程在GPU上加速。本发明在子节点剖分子模块中，通过设置扫描单元，预先计算出父节点内包围盒数组复制到左、右子节点时的位置索引，再使用数组索引单元根据索引将数组中的数据分别复制到左、右子节点内，在GPU上并行实现了父节点内包围盒数组复制到左、右子节点的过程，加速了KD树的构建。

实施例5：

本发明还是一种KD树叶节点内包围盒数组统一输出方法，是在基于GPU的KD 树的射线追踪加速系统中的KD树加速搜索生成模块中实现叶节点内包围盒的统一输出，基于GPU的KD树的射线追踪加速系统同实施例1-4，

参见图2，包括有如下步骤：

1)输入数据并初始化：输入三角面元包围盒模型数组，以及两个辅助数组，辅助数组的长度和KD树每层的包围盒数组等长，用于辅助划分当前节点内的包围盒数组，初始化根节点作为当前节点。

2)终止剖分判断：判断当前节点是否满足终止剖分条件，当前节点地址位置中包含的包围盒数量是否小于终止剖分阈值，若满足，执行步骤7)，若不满足，执行步骤 3)。

3)计算最佳剖分面位置：用中分法或表面积启发式算法(SAH)计算最佳剖分面和剖分轴。

4)辅助数组赋值：利用两个辅助数组辅助划分，根据最佳剖分面位置与剖分轴，位于最佳剖分面左侧的包围盒，在第一个辅助数组中将对应位置标为1；位于最佳剖分面右侧的包围盒，在第二个辅助数组中将对应位置标为1；被最佳剖分面切分的包围盒，在两个辅助数组中将对应位置均标为1。

5)扫描辅助数组：计算当前节点中的数组复制到左、右子节点时的起始存储地址，并得到数组中每个元素在左、右子节点对应地址中的位置索引。

6)剖分当前节点：将当前节点剖分为左、右两个子节点，根据步骤5)中得到的位置索引，将当前节点内包围盒数组复制到左、右节点存储地址的相应位置中，更新左、右子节点相关信息，并将左、右子节点分别作为当前节点进入步骤2)，再次进入 KD树构建循环判断。

7)形成输出数组：将当前节点标记为叶节点，将叶节点存储的地址内的数组附加到输出数组尾部，并将叶节点内的数组在输出数组中的地址记录在叶节点内，留给射线追踪模块使用。

8)终止剖分过程：当所有的节点都终止剖分之后，得到一棵构建完成的高质量 KD树，以及重新组织后的包围盒模型数组。

参见图4，本发明中，在GPU上使用扫描单元与数组索引单元加速了当前节点内数组复制到左、右子节点内的过程；当子节点无法再剖分时，使用KD树叶节点内包围盒统一输出方法，将节点内数据统一输出。由于仅在叶节点内保存数据在输出数组中的地址段信息，同时各叶节点内的数组在输出数组中是连续的，因此在射线追踪过程中可以实现数据的合并访问，加速射线追踪的过程。

在建树过程中，本发明不再在节点内存储数据的索引，改为将数组统一存储到显存预申请模块申请的显存中，只在节点内存储该节点所包含数组在显存中的地址段。当节点为叶节点时，将叶节点内的数组复制到输出数组中，并在叶节点内记录叶节点内的数组在输出数组中的地址段。本发明可以节省KD树在GPU中所占显存的大小；可以在节点内的包围盒数组中的元素数量大于终止剖分阈值时停止剖分过程，并且不会出现数据溢出的情况；使得后续射线追踪模块的射线求交过程可以合并访问叶节点包含的数据，从而加速射线追踪的过程。

本发明在较短的时间内输出一棵质量较高KD树。将KD树结构线性保存在显存中，将KD树的根节点和重新组织过后的输出数组作为射线追踪模块的两路输入，以查表的方式使用。

下面给出一个更加详尽的例子，对本发明进一步说明：

基于GPU的KD树的射线追踪加速系统及KD树叶节点内包围盒统一输出方法同实施例1-5，参见图1：

显存预申请模块：该模块的输入为三角面元模型的三角面元的数量，在该模块内部根据三角面元的数量计算出三角面元包围盒模型数组所占显存大小，在CPU端预先一次性申请GPU上的显存空间，该显存空间至少足以无冲突生成KD树，该模块的输出是预申请的显存的起始地址。

数据预处理模块：在GPU上，使用显存预申请模块申请的显存，对输入该模块的三角面元模型的每个三角面元并行构建包围盒对于每一个三角面元，分配一个线程，计算该三角面元的轴对称(AABB)包围盒，以下统称包围盒。并初始化当前入射角度，输出包围盒模型数组及当前入射角度。

KD树加速搜索生成模块：该模块是在GPU上并行生成的，KD树构建完成后，在 KD树上遍历可以大大加速射线追踪过程，在KD树生成的过程中，以地址访问的方式动态使用显存预申请模块申请的显存，加速了KD树的生成，数据预处理模块中输出的包围盒模型数组为该模块的输入，该模块的两路同步输出，一路为对应于该包围盒模型数组的KD树，另一路为重新组织过的包围盒模型数组，该两路同步输出直接作为射线追踪模块的输入。

虚拟孔径面生成模块：该模块的输入为数据预处理模块输出的包围盒模型数组及当前入射角度，在该模块内计算包围盒模型的边界值，并在GPU上并行生成声线管束，该模块的输出为若干个相同尺寸的声线管束。

射线追踪模块：该模块的输入分别为虚拟孔径面生成模块输出的声线管束和加速搜索生成模块的两路同步输出，该三路为同步输入，在GPU上并行加速声线管束的射线追踪，射线追踪过程是通过遍历KD树完成的，该模块的输出为声线管束在遍历KD 树后的反射信息。

Gordon积分模块：Gordon积分模块接收追踪模块输出的反射信息，使用Gordon 积分公式在GPU上对声线管束进行并行积分，使用归约操作对所有声线管束的积分结果进行并行求和，该模块的输出为当前角度下所有声线管束积分值的代数和。

界面显示模块：Gordon积分模块输出的积分值及相关信息均显示在界面上，便于信息查看与调试。

显存预申请模块中，本例中的三角面元模型是潜艇的模型，如图5所示，潜艇模型的尺寸为：62m×7.5m×11m，模型由147427个节点，294850个三角面元组成。其中三角面元尺寸为2λT(本例中λT＝5cm)；预申请的显存空间大小为 (D·P·I)+2F+(P·I)

其中，D为预估的最大树深,本例中D取20，P为膨胀系数，本例中M取10，I 为根节点中包围盒数组所占显存的大小，F为一个辅助数组所占显存的大小，该模块为KD树的每层都分配了大小为(P·I)的显存空间，同时为输出数组申请了大小为 (P·I)的显存空间；在剖分父节点时，根据剖分面及剖分轴，将父节点的包围盒模型数组拷贝至属于子节点的显存块中，因此在并行拷贝的过程中，不会发生显存访问冲突；此外，膨胀系数的设定也保证了被剖分面分割的三角面元包围盒可以被安全的放置到两个子节点中。

数据预处理模块中，根据潜艇模型的三角面元数量，为每个三角面元都开辟了一个线程，用来计算个三角面元在三维空间内的坐标极值(Xmax,Ymax,Zmax)和 (Xmin,Ymin,Zmin)，其中Xmax,Ymax,Zmax分别为三角面元的三个顶点在X,Y,Z轴上的最大值，Xmin,Ymin,Zmin分别为三角面元的三个顶点在X,Y,Z轴上的最小值。由每个三角面元的两个坐标极值，可以确定每个三角面元的紧致包围盒的大小和在空间中的位置；初始化当前入射角度为0度。

KD树加速搜索生成模块是利用KD树叶节点内包围盒数组统一输出方法生成的 KD树，对KD树加速搜索生成模块中设有的子模块分述如下：

条件判别子模块首先接收数据预处理模块输出的数据，根据当前节点内包围盒的数量Box_Num判断决定当前节点数据是进入节点剖分子模块还是终止剖分子模块；当 Box_Num小于等于停止剖分阈值Stop_Num时，进入终止剖分子模块，本例Stop_Num 取32。如果进入终止剖分子模块，则停止当前节点的剖分。当Box_Num小于等于停止剖分阈值Stop_Num时，进入节点剖分子模块。如果进入节点剖分子模块，则之后依次进入子节点内包围盒数组生成子模块、子节点信息更新子模块，最后再返回到条件判别子模块对子节点进行再次判断。条件判别子模块循环接收节点剖分子模块反馈的数据，反复循环判断，逐层构建KD树结构，直至所有节点都进入终止剖分子模块。

节点剖分面计算子模块的输入为条件判别子模块的输出，该子模块根据接收到的数据计算当前节点的剖分面位置和剖分轴，将一个父节点剖分成两个子节点，子节点中包含有父节点的部分信息，并将计算结果和该子模块接收到的数据输出到子节点内包围盒数组生成子模块。

节点剖分面计算子模块设有两个单元，根据输入该子模块的面元数量N决定使用哪个单元进行节点剖分。

当节点内面元数量大于阈值Big_Node时，使用大节点剖分单元对当前节点进行剖分，本例阈值Big_Node取256，该单元使用中分法确定剖分面位置及剖分轴，剖分轴为节点内包围盒最长边所平行的坐标轴，剖分面位置坐标为节点内包围盒最长边的中点坐标值，由于该过程需要求得节点内包围盒的边界，因此在GPU端使用归约操作求边界值可以加快求解过程。

当节点内面元数量小于等于阈值Big_Node时使用小节点剖分单元对当前节点进行剖分，该单元使用SAH算法确定剖分面位置及剖分轴。在SAH算法中需要SAH代价函数计算每个候选剖分面的代价值SAHcost，并选取使代价值SAHcost最小的那个候选剖分面作为最佳剖分面。

SAH代价函数为：

其中，ctraversal表示射线遍历到当前节点的代价，SA(V)表示当前节点内包围盒的表面积，SA(VL)表示当前候选剖分面下左子节点内包围盒的表面积，NL表示在当前候选剖分面下被划分到左子节点的面元数量，SA(VR)表示右子节点内包围盒的表面积，NR表示被划分到右子节点的面元数量，Chit表示射线与三角面元求交的代价值。

从公式可以看出，每个候选剖分面的SAH值计算都是独立进行的，所以SAH值的计算可以利用GPU进行加速，为每个候选剖分面分配一个block，用来计算该候选剖分面的SAH值，并在GPU端通过归约算法到最小的SAH值及相应的候选剖分面。当剖分面确定之后，剖分面所垂直的坐标轴即为剖分轴。

子节点内包围盒数组生成子模块的输入为节点剖分面计算子模块的输出，根据输入，在GPU中将当前节点内包围盒数组划分为两个子数组，并行输出生成的两个子数组分别至两个子节点信息更新子模块；对子节点内包围盒数组生成子模块内设有的单元分述如下：

扫描单元：该单元的输入为节点剖分子模块的输出，在GPU上使用包含扫描操作，并行计算出两个子节点内包围盒数组元素在当前节点内包围盒数组中的索引，并将索引保存在显存预申请模块申请的辅助数组对应的位置上。

数组索引单元：该单元的输入为扫描单元输出的索引，在GPU上并行索引，将当前节点存储的地址中的数组按照索引分别复制到两个子节点存储的地址中，将当前节点存储地址内包含的数组划分至两个子节点存储的地址内，然后将两个子节点分别输出至条件判别子模块进行重新判断。

子节点信息更新子模块：该子模块为两个结构相同并列的子模块，各自的输入为子节点内包围盒数组生成子模块输出的一个子数组，当前节点的子节点先继承当前节点的相关信息，然后根据接收到的子数组更新子节点的相关信息，并将子节点作为当前节点与子数组输出至条件判别子模块进行重新判断；再次进入KD树构建循环过程。

剖分终止子模块：该子模块的输入为条件判别子模块输出的数据，该子模块终止当前节点内包围盒数组的剖分，并将当前节点包含的子数组输出至输出数组中，等待全部节点停止剖分后一次性输出。

KD树加速搜索生成模块最终输出一棵构建完成的KD树以及经过重新组织的包围盒数组。经过重新组织的包围盒数组是使用KD树叶节点内包围盒数组统一输出方法进行重新组织并输出的，包括有如下步骤：

1)输入数据并初始化：输入三角面元包围盒模型数组，以及两个辅助数组，辅助数组的长度和KD树每层的包围盒数组等长，用于辅助划分当前节点内的包围盒数组，初始化根节点作为当前节点。

2)终止剖分判断：判断当前节点是否满足终止剖分条件，当前节点地址位置中包含的包围盒数量是否小于终止剖分阈值，若满足，执行步骤7)，若不满足，执行步骤 3)。

3)计算最佳剖分面位置：用中分法或SAH法计算最佳剖分面和剖分轴。

4)辅助数组赋值：利用两个辅助数组辅助划分，先将辅助数组清零，然后根据最佳剖分面位置与剖分轴，位于最佳剖分面左侧的包围盒，在第一个辅助数组中将对应位置标为1；位于最佳剖分面右侧的包围盒，在第二个辅助数组中将对应位置标为1；被最佳剖分面切分的包围盒，在两个辅助数组中将对应位置均标为1。

5)扫描辅助数组：计算当前节点中的数组复制到左、右子节点时的起始存储地址，并得到数组中每个元素在左、右子节点对应地址中的位置索引。

6)剖分当前节点：将当前节点剖分为左、右两个子节点，根据5)中得到的位置索引，将当前节点内包围盒数组复制到左、右节点存储地址的相应位置中，更新左、右子节点相关信息，并将左、右子节点分别作为当前节点进入步骤2)，再次进入KD 树构建循环判断。

7)形成输出数组：将当前节点标记为叶节点，将叶节点存储的地址内的数组附加到输出数组尾部，并将叶节点内的数组在输出数组中的地址记录在叶节点内，留给射线追踪模块使用。

8)终止剖分过程：当所有的节点都终止剖分之后，得到一棵构建完成的高质量 KD树，以及重新组织后的包围盒模型数组。

建树过程中，不再在节点内存储数组中元素的索引，改为将数据统一存储到显存预申请模块为该层的KD树节点预申请的显存中，只在节点内存储所包含数据在显存中地址的起点和终点。当节点为叶节点时，将显存预申请模块为该层的KD树节点预申请的显存中的数据复制到输出数组中，并在叶节点内存储数据在输出数组中的地址。因此也能节省KD树在GPU中所占显存的大小；在当前节点内包围盒数量大于终止剖分阈值时提前停止剖分过程，不会出现数组溢出的情况；后续射线追踪模块可以合并访问叶节点包含的数据，加速了射线追踪过程。

本发明在较短的时间内输出一棵质量较高KD树。将KD树结构线性保存在显存中，将KD树的根节点和重新组织过后的输出数组作为射线追踪模块的两路输入，以查表的方式使用。

步骤1)生成包围盒的过程是在GPU上并行完成的，在GPU上为每个三角面元开辟一个线程计算各自的包围盒，并按照原三角面元数组的顺序排列组织成一个包围盒数组。

步骤2)中若当前节点内包围盒的数量小于等于Stop_Num，则进入步骤7)停止当前节点的剖分；若当前节点内的包围盒数量大于Stop_Num，则进入步骤3)，将当前节点作为父节点剖分成左右子节点。

步骤3)中若当前节点内面元数量大于阈值N(本例中N为256)，则视当前节点为大节点，使用中分法确定剖分面位置及剖分轴，剖分轴为节点内包围盒最长边所平行的坐标轴，剖分面位置坐标为节点内包围盒最长边的中点坐标值，由于该过程需要求得节点内包围盒的边界，因此在GPU端使用归约操作求边界值可以加快求解过程。

若当前节点内面元数量小于阈值N，则视当前节点为小节点，使用SAH算法确定剖分面位置及剖分轴。在SAHcost SAH算法中我们需要SAH代价函数计算每个候选剖分面的代价值，选取使代价值SAHcost最小的那个候选剖分面作为剖分面。

计算SAH代价的函数为：

其中，ctraversal表示射线遍历到当前节点的代价，SA(V)表示当前节点内包围盒的表面积，SA(VL)表示当前候选剖分面下左子节点内包围盒的表面积，NL表示在当前候选剖分面下被划分到左子节点的面元数量，SA(VR)表示右子节点内包围盒的表面积，NR表示被划分到右子节点的面元数量，Chit表示射线与三角面元求交的代价值。

从公式可以看出，每个候选剖分面的SAH值计算都是独立进行的，所以SAH值的计算可以利用GPU进行加速，为每个候选剖分面分配一个线程，用来计算该候选剖分面的SAH值，并在GPU端通过归约算法到最小的SAH值及相应的候选剖分面。

当剖分面确定之后，剖分面所垂直的坐标轴即为剖分轴。

步骤4)两个辅助数组是在显存预申请模块中预申请的两个长度相等的数组，每个辅助数组的长度为三角面元模型中三角面元数量的P倍；因此在KD树的逐层构建过程中，每层的节点只需要通过地址访问辅助数组的对应位置，便可实现节点间的无冲突使用；根据最佳剖分面+位置与剖分轴，位于最佳剖分面左侧的包围盒，在第一个辅助数组中将对应位置标为1；位于最佳剖分面右侧的包围盒，在第二个辅助数组中将对应位置标为1；被最佳剖分面切分的包围盒，在两个辅助数组中将对应位置均标为1；通过辅助数组，便可知道包围盒数组中对应位置上的包围盒归属于哪个子节点；

步骤5)中，对辅助数组进行包含扫描操作，包含扫描操作会将原辅助数组变成一个新辅助数组，新辅助数组中的每一个位置上的新元素都是原辅助数组中从下标为0 的位置到该位置所有元素的和，因此，通过计算新辅助数组便可得到当前节点的包围盒复制到左、右子节点内包围盒数组时的对应下标，方便并行复制数据；

步骤6)中，根据步骤5)计算出的新辅助数组存储的索引下标，将当前节点内的包围盒数组分别复制到左、右子节点的对应索引位置上，并将当前节点剖分为左、右两个子节点，根据划分至左、右子节点内的包围盒数组，更新左、右子节点相关信息，并将左、右子节点分别作为当前节点进入步骤2)，再次进入KD树构建循环判断，直至所有的当前节点都进入步骤7)，成为叶节点。

步骤7)中，终止当前节点的剖分，将当前节点标记为叶节点，并将该叶节点内存储的地址内的数据附加到输出数组的末尾，并将数据在输出数组中的地址段记录在该叶节点内，待所有当前节点都被标记为叶节点时，进入步骤8)，完成KD树的构建。

步骤8)中，当所有的节点都终止剖分之后，得到一棵构建完成的高质量KD树，以及所有叶节点内的包围盒数组拼接在一起的输出数组，保存KD树的结构与输出数组在显存中，在射线追踪模块中作为查表使用。更新子节点的线索值，子节点内包围盒的边界值，以及子节点所包含的面元，实现了叶节点内包围盒数组的统一输出。

虚拟孔径面生成模块中，虚拟孔径面为benchmark在等相位面上的投影，虚拟孔径面由大量孔径相等的矩形声线管束组成，声线管束的孔径λs，本例中λs取0.0625

λT，可以根据实际需要设定大小，虚拟孔径面是在GPU端并行生成的，具体生成步骤如下：

首先要根据当前的入射角度、与模型的距离及输入的包围盒模型数组计算出虚拟孔径面的法向量及尺寸长L及宽W；

根据虚拟孔径面的尺寸L*W以及声线管束的尺寸λs，计算出声线管束的数量其中除法和只保留整数位；为GPU核函数开辟Snum个线程，每个线程生成一个声线管束，并将生成的声线管束输出至射线追踪模块，根据线程的索引将生成的声线管束存入相应地址的显存中，等待射线追踪模块的读取。

射线追踪模块中，该模块的输入分别为虚拟孔径面生成模块输出的声线管束和加速搜索生成模块的两路同步输出，该三路为同步输入，在GPU上为各个声线管束的每个射线开辟一个线程，并行执行射线追踪核函数，射线追踪过程即是以遍历二叉树的形式遍历KD树，模拟出光线照射到三角面元模型上并产生反射的物理过程，该模块的输出为虚拟孔径面中每根射线与三角面元模型的交点坐标及射线的传播路径。

在射线追踪的过程中，射线之间是相互独立的，在GPU上为每根射线分配一个线程，独立执行射线追踪模块，并将追踪结果存入对应显存地址中，等待后续积分使用。在串行的代码中，我们可以使用递归的方法实现射线追踪，但是GPU不支持递归，因此我们设置了一个最大追踪次数Trace_Num，当射线追踪函数执行了 Trace_Num次之后，当前射线就会自动停止追踪，防止陷入死循环中。

射线追踪模块中，需要对声线管束判断是否发生射线分裂，具体有以下过程：

判断声线管束是否满足射线分裂条件的原则是：声线管束的矩形的四个顶点是否都照射到同一个三角面元上。当四个顶点照射到同一个三角面元上时，不需要射线分裂；当四个顶点没有照射到同一个三角面元上时，声线管束发生分裂，声线管束需要进行均匀四等分，即由一个声线管束等分为四个子声线管束。重新追踪子声线管束，若声线管束分裂三次后，仍满足分裂条件，则将仍然满足分裂条件的子声线管束抛弃，不再进行射线追踪与后续计算。对于不再发生分裂的声线管束以及子声线管束，根据反射信息判断是有效的声线管束。声线管束有效与否的判断原则是：声线管束的四个顶点发出的射线是否照射到同一个三角面元上。若是，则标记为有效声线管束，返回反射信息；若不是，则标记为无效声线管束，不返回反射信息。

本发明的Gordon积分模块接收追踪模块输出的相关信息，使用Gordon积分公式在GPU上对声线管束进行并行积分，使用归约操作对所有声线管束的积分结果进行并行求和，该模块的输出为当前角度下所有声线管束积分值的代数和，具体过程如下：

本发明为每个有效的声线管束分配一个线程进行一次积分运算，并将积分结果作为中间值存储在显存中，当所有的有效声线管束积分完毕之后，使用归约操作对积分结果进行求和运算，得到当前入射角度下的TS值。时域物理声学应用于弹跳射线法法的声场积分公式为；

其中，i，j表示声线管束在虚拟孔径面上的第i行第j列，k为波数k＝2π/λ， K＝2k，Δa'n表示四边形各条边的向量，Δρ'n表示积分区域内任意一点到参考点的矢量Δρ'n＝(a'n+a'n+1)/2，T表示r'i,j在积分区域上的投影长度，RM为源点位置矢量，r'i,x表示接收点到参考点的单位向量

r'i,j＝u'x'+v'y'+w'z'

其中x',y',z'是对原始空间坐标系{O,X,Y,Z}进行投影变换后得到的三维平面坐标系{M,X',Y',Z'}，从三维平面坐标系{M,X',Y',Z'}到原始坐标系{O,X,Y,Z}的转换矩阵Q

其中，M表示原始坐标系下的四边形中心点同时作为三维平面坐标系下的坐标原点，X'表示点M到四边形任意顶点的单位向量，Z'表示四边形平面的单位法向量 n0，Y'＝X'×Z'。

将原始坐标系下的任意点乘以转换矩阵Q的逆矩阵Qinv，即完成了该点从原始坐标系到三维平面坐标系下的坐标投影变换。

将所有射线的反射参数代入所述公式，根据高频近似原理，计算所有三角面元积分结果的代数和即可得到当前入射角度下声场强度的TS值。

本发明的界面显示模块中，将当前入射角度下Gordon积分模块输出的积分值、benchmark模型，三角面元模型中三角面元的数量，虚拟孔径面尺寸，虚拟孔径面射线数量在界面上显示出来，便于信息查看与调试。

以benchmark模型的坐标原点作为顶点，以声源位置和坐标原点的连线为轴逆时针转动θ(本例θ取0.5)度，改变当前入射角度；顺时针或逆时针转动均可；每次转动的角度θ可以根据需要设定，最终完成对benchmark0-180度的积分计算。

下面通过仿真对本发明的技术效果再做说明：

实施例7：

基于GPU的KD树的射线追踪加速系统及KD树叶节点内包围盒数组统一输出方法同实施例1-6，

使用基于GPU的KD树的射线追踪加速系统来仿真声线管束在水下目标上的传播过程，并将追踪结果进行积分求和，被遮挡面元也能参与回波强度计算，因此可以获得较传统方法更为准确的仿真预报结果。

图6是使用本发明对图5进行仿真的结果，参见图6，其中横坐标是入射角度，纵坐标是当前入射角度下的声场强度值(单位为分贝)；由图6可以看出，声场强度值随着入射角度从0-180度的改变是先增大再减小的过程。在0度和180度附近时，由于潜艇头部与尾部横截面积较小，因此生成的虚拟孔径面也比较小，从而导致返回的声波强度值较低；随着入射角度的增大，声源逐渐位于潜艇的艇身侧面，此时的虚拟孔径面较大，产生的声线管束较多，声场强度值就比较大，由于潜艇表面反射的关系，因此声场强度达到峰值的角度并不是90度，而是在110度左右。

相对于传统板块元及改进板块元计算方法对被遮挡面元采取抛弃策略而导致的精度下降，本发明使用射线追踪技术获得被遮挡面元反射的声波，对于表面具有较多遮挡的凹面目标，被遮挡面元也能参与计算，因此使用本发明的系统可以获得较高的计算精度。

图7是使用本发明在十块GPU计算卡上对图5进行仿真产生的相关数据(GPU型号为Nvidia P100)，参见图7，其中横坐标是入射角度，左侧的纵坐标是射线追踪模块在当前入射角度下的耗时(单位为毫秒)，右侧的纵坐标是虚拟孔径面生成模块在当前角度下生成的声线管束的数量，图7中虚线为本发明的射线追踪模块在不同入射角度下的计算用时，实线为本发明的虚拟孔径面生成模块在不同入射角度下生成的声线管束的数量。从图7中可以看出，本发明中射线追踪模块的用时和虚拟孔径面生成模块生成声线管束的数量在大体趋势上是线性的关系，说明射线追踪模块耗时仅与声线管束的数量有关，与入射角度无关，证明本发明构建KD树时既有较快的速度，又有较好的质量。相对于传统的板块元需要先提取无遮挡面元再进行射线追踪的过程，本发明省略了无遮挡面元的提取过程，且加速了射线追踪的过程，所以本发明对目标声场强度仿真预报的速度有了极大的提高。

简而言之，本发明的基于GPU的KD树的射线追踪加速系统及KD树叶节点内包围盒数组统一输出方法，解决了对水下复杂目标的声场强度快速预报问题，系统按照信号连接方向依次包含：显存预申请模块，数据预处理模块，KD树加速搜索生成模块，虚拟孔径面生成模块，射线追踪模块，积分模块，界面显示模块。本发明增设了显存预申请模块，并在KD树加速搜索生成模块中加速了节点剖分的过程，本发明还提出了一种KD树叶节点内包围盒数组统一输出方法，包含以下步骤：输入数据并初始化，终止剖分判断，计算最佳剖分面位置，辅助数组赋值，扫描辅助数组，剖分当前节点，形成输出数组，终止剖分过程。本发明能够对给定的任意模型的三角面元模型进行声场强度计算。本发明使用基于GPU的KD树的射线追踪加速系统仿真声线的传播过程，无需对遮挡面元进行判断，具有计算速度快，精度高，对不同目标的适应性强的特点，可有效用于对水下目标强度的快速预报仿真。