提出了一种分频段电磁环境绘制算法.以雷达设备的电磁辐射功率密度表征电磁环境, 计算数据时不区分频段, 获得所有频段联合作用的数据场, 根据雷达设备的分布及作用情况将数据场划分为对应不同频段的子区域.使用光线投射算法绘制电磁环境数据场, 在设计传递函数时生成多频段颜色查找表, 以颜色区分频段, 并绘制数据场.在所有颜色查找表中保持前部若干颜色相同, 绘制所有频段的联合分布.最后, 使用统一计算设备架构加速实现了实时分频段绘制.实验结果表明, 绘制结果显示了电磁环境中的频段分布信息, 又反映了所有频段的联合作用情况, 方便用户直观理解电磁环境.
A band rendering method for electromagnetic environment was presented. The electromagnetic environment is represented by electromagnetic radiation power density of radars, when the data band computing is carried, a volumetric data for all bands is got. The volumetric data was then divided into sub-regions for each band based on each radar effect. Ray casting algorithm was used to render electromagnetic environment in which the transfer function was expressed as a set of color lookup tables for each band, different band has different color. Among all color lookup tables, there are some same colors for rendering the combined effects of all bands. A real-time rendering method was carried out using compute unified device architecture. Experiments show that the rendering images can display the distribution of bands and the combined effects of all bands for electromagnetic environment, and are helpful for understanding the electromagnetic environment intuitively.
电磁环境已经成为现代信息化战场的重要组成部分，将电磁环境展示给指挥员是战场环境可视化研究中的一项重要内容.电磁环境可视化主要有二维表现形式和三维表现形式[1].三维表现形式较常使用，等值面提取[2]和体绘制方法[3-4]是能展现电磁环境整体信息分布的三维可视化方法.
但是，这些已有的方法中均未考虑电磁环境的频段信息.针对此情况，笔者研究了区分频段的电磁环境绘制方法，使用效果最好的光线投射算法[5]进行绘制.光线投射算法计算量大，基于图形硬件加速的光线投射算法[6]较常见，特别是统一计算设备架构(CUDA, compute unified device architecture)提出后，使用图形硬件加速光线投射算法更加普遍.
光线投射算法的关键环节是传递函数设计[7]，即设计数据到颜色和不透明度值的转换规则.研究中，通过对数据场的频段划分和光线投射算法中传递函数的分频段设计，最终在CUDA下实现了电磁环境的实时分频段绘制.
1 算法概述算法主要包括数据预处理和分频段绘制2个阶段，其流程如图 1所示.数据预处理阶段主要计算生成包含频段信息的电磁环境数据场；分频段绘制阶段主要使用光线投射算法绘制电磁环境数据场，通过设计传递函数将不同频段的数据映射为不同颜色，在绘制结果中区分出不同频段.
|
图 1 算法流程 |
本研究使用公式计算生成电磁环境数据场，其无噪声，便于分析数据场进而设计高效的绘制算法.
2.1 数据场采样计算电磁设备的辐射信号叠加在一起构成电磁环境，但用频设备类型广、数量多，电磁环境非常复杂且其精确计算方法仍处于研究当中.研究中重点关注战场环境中最重要的电磁设备——雷达，使用多部雷达的联合电磁辐射功率密度空间分布表征雷达作用产生的电磁环境.电磁环境区域中某点A的数据值SA由区域中的n部雷达设备对此点处的辐射功率密度Si(i=1, 2，…, n)叠加得到，计算式为
|
(1) |
其中：Pi为雷达i的发射功率，Gi为雷达i的天线增益，Ri为雷达i与点A间的距离；iωA表示在地形的影响下雷达i能照射到A点，
只研究覆盖一定频率范围的电磁环境，整体考虑这一频率范围的同时对其中各个子频段在绘制效果上进行区分.因此，数据场计算时不区分频段：只要雷达设备能照射到A点即参与A点的电磁环境值计算.如此简化了数据计算和存储，方便绘制操作.分频绘制将在频段划分和绘制算法中体现.
研究有限区域的电磁环境，在区域内离散采样获得数据场，使用离散的数据采样点值表示区域电磁环境.电磁环境分布于地理空间中，但球坐标系下等坐标间隔对应的间隔距离不等，影响采样点的均匀分布.因此，在直角坐标系下对区域进行离散采样.离散采样如图 2所示，将立方体区域分割为若干个相等的正方体，每个正方体的所有顶点都对应有采样点.
|
图 2 区域离散采样示意图 |
计算得到的电磁环境数据场中，每个采样点均包含各个频段的数据.因此，不能依据采样点的数据值区分其所属频段.根据采样点在数据场中的位置判定其所属频段，不同频段对应数据场中不同的子区域，进而绘制出频段信息的空间分布.
区域中的一个数据点必然属于某个频段对应的子区域，将此点处辐射功率密度值最大的频段设定为其所属频段.由式(1) 可知，计算各点处的辐射功率密度值时只有距离是变化的，因此某点处某个频段的辐射功率密度值与距离的关系为
|
(2) |
其中：Sn为频段n的辐射功率密度值，k为属于频段n的设备下标，max{PG}表示所有设备发射功率与天线增益乘积的最大值.
将式(2) 中的βk命名为“距离增强系数”，可在数据计算阶段计算出全部设备的β值.由式(2) 可知，判断数据点的频段时，只需计算其与所有设备的距离，然后计算出对应频段的
根据距离判断数据点的频段归属可在数据计算阶段完成：计算得到数据场后，遍历所有数据点计算其频段归属，并将频段信息存储到数组中.执行绘制操作时只需查询数组便可得到频段信息，提高绘制时的频段划分效率，减少绘制操作的时间开销.
3 分频段显示绘制使用光线投射算法对电磁环境数据场进行绘制，算法中的传递函数经过特殊设计，能在绘制结果中区分显示电磁环境数据中的频段分布情况.
3.1 传递函数设计光线投射算法中，传递函数决定了数据属性到颜色属性的映射过程.研究中对传递函数进行设计，在生成的绘制结果中以颜色区分不同频段.
直接将数据采样点值映射为颜色和不透明度值，映射规则：将数据的值域[dmin, dmax]等分成n段，其中第i段[dmini, dmaxi)中所有数据值都对应于颜色Ci和不透明度αi，最大数据值dmax对应于颜色Cn和不透明度αn.将传递函数使用的n个颜色和不透明度值顺序存储于颜色查找表中，绘制操作执行时直接查找使用.
使用线性分段函数生成所使用的n个颜色和不透明度值.以颜色值为例，指定第1个颜色值为C1、最后1个颜色值为Cn以及之间的m个颜色值，对其中任意2个相邻指定颜色值Ca和Cb之间的颜色值Ct由式(3) 计算得到.参与式(3) 运算的每个颜色值C都含有r、g、b3个分量.生成n个不透明度值{αi|i=1, …, n}的过程与颜色值的生成过程相同.
|
(3) |
为了在绘制结果中区分不同频段，为每个频段生成各自的颜色查找表，确定采样点的所属频段后，在此频段对应的颜色表中查找其对应的颜色值.
若不同频段映射为完全不同的颜色，在频段区域相接的地方会出现颜色变化不连续的情况.将所有频段区域边缘附近的数据映射为相同颜色可以避免出现此情况.由式(1) 可知，数据点值跟其与设备距离的二次方成反比，而频段区域边缘相对区域中心处数据值较小.因此，将所有频段较小的数据值映射为相近或相同的颜色，即将所有颜色表的前k个颜色值设置为相同，从第k+1个颜色值开始逐渐向不同的颜色过渡. 图 3为2个频段使用的颜色查找表示意图，从绿色分别过渡到蓝色(频段1) 和红色(频段2)，其中前2个示例颜色相同.实际使用中，用户可以交互设置相同颜色的个数.
|
图 3 两频段颜色查找表示意图 |
在各个颜色查找表中，相同编号的颜色对应相同的不透明度值，从而不同频段的相同数据值对应不同的颜色值和相同的不透明度值，便于不同频段数据之间的对照.
所有颜色查找表生成后，即确定了数据采样点到颜色和不透明度值的映射规则.将1条光线上所有数据采样点对应的颜色和不透明度值累积为1个像素值即完成了从数据点到像素值的转换.对1条光线上的采样点序列{Li}，使用CiL表示点Li对应于颜色查找表中的颜色值，初始化Cend为0用于存储此光线对应像素点的颜色值，依次遍历序列{Li}，累积得到Cend的规则为
|
(4) |
其中：C.rgb为颜色的RGB分量，C.α为颜色的不透明度分量.
3.2 并行绘制策略使用CUDA实现并行光线投射算法，对电磁环境数据场进行分频段显示绘制.进行并行光线投射前，需要将所有绘制需要的数据存储到GPU中，以便执行绘制操作时供GPU高效地调度和使用.将数据场绑定到三维纹理，使用硬件支持的纹理采样器根据坐标对数据进行高效采样.
光线投射算法中光线相互独立，为每条光线开启一个线程，每个线程对应生成一个结果图像上的像素点，所有线程并行执行完成后即获得结果图像.光线投射线程执行流程的伪代码如下.
内核函数：光线生成像素点
pixX=blockIdx.x*blockDim.x+threadIdx.x; //确定像素点位置
pixY=blockIdx.y*blockDim.y+threadIdx.y; //确定像素点位置
if (pixX, pixY) is in the Image then
theRay=getRayforPix (pixX, pixY); //生成对应光线
pixValue=castRay (theRay); //光线穿越数据场获取像素值
Image [pixX, pixY]= pixValue;
其中getRayforPix()函数根据视点位置和像素坐标生成对应的光线，光线theRay的参数主要包括其穿越数据场时的入射点坐标、出射点坐标及传播方向；castRay()函数完成光线穿越数据场时对数据的采样、颜色映射、颜色累积等操作，最终生成光线对应的像素值.
castRay()函数是分频段绘制的关键操作，其具体执行分为以下几个步骤.
步骤1 初始化颜色和不透明度的累积值为零；
步骤2 计算光线当前对应采样点的纹理坐标；
步骤3 判断坐标是否在纹理内部，是则执行步骤4，否则返回颜色和不透明度的累积值；
步骤4 纹理采样获取采样点数据值；
步骤5 获取采样点所在频段信息，使用相应频段颜色查找表进行颜色映射；
步骤6 使用式(4) 对颜色、不透明度进行累积；
步骤7 判断不透明度是否达到1.0，是则返回颜色和不透明度的累积值，否则光线沿传播方向增加一个步长继续执行步骤2.
在castRay()函数的第3和第7步中，返回颜色和不透明度的累积值后函数执行结束，返回值为对应像素值.在第7步中，当不透明度累积达到1.0后光线提前截止，能减少不必要的后续投射计算，节省执行时间[8].
4 实验结果与分析算法在Microsoft Visual Studio 2005中编程实现，其中为支持CUDA编程集成安装了NVIDIA GPU Computing SDK 4.2.算法执行环境的硬件配置为Intel(R) Core(TM) i5-2400 CPU 3.10 GHz、4 GB内存、NVIDIA GeForce GT 530显卡.
实验中设置了18台分属2个频段的雷达设备，对应的2个颜色查找表都从绿色开始，一个渐变到蓝色，另一个渐变到红色.数据场的采样点规模为128×128×64，光线采样步长与采样点沿坐标轴方向的间隔相同，生成结果图像分辨率为1 200×900.
图 4所示为2个频段数据场在不同视角下的绘制结果，其中最左图显示了不同频段数据与颜色间的映射关系. 图 4所示绘制结果在显示出2个频段分布情况的同时又直观展现了各个雷达设备的相对覆盖范围，达到了分频段绘制显示电磁环境的目的，其中，绿颜色部分显示出了2个频段的联合电磁环境分布情况，这是各个频段颜色查找表前部若干颜色相同产生的结果.充分说明了2.2节频段划分和3.1节传递函数设计的作用，不仅对不同频段进行了区分绘制，同时还绘制了不同频段联合作用的电磁环境，能直观描述电磁态势.此外，图 4中每幅子图下方标示的时间是从数据场生成图像并显示出来这一过程的算法执行时间，可以看出完全达到了实时绘制电磁环境数据场的要求.
|
图 4 2频段电磁环境绘制效果 |
当所有雷达设备分布于3个频段时，使用3个颜色查找表. 图 5所示为3个频段的数据场绘制结果，其中最左图在图 4的基础上增加了一个新的频段颜色查找表. 图 5所示的绘制效果展示了3个频段在电磁环境中的分布情况，达到了分频段绘制电磁环境的目的.
|
图 5 3频段电磁环境绘制效果 |
由于频段划分在数据处理阶段完成，3频段绘制效率与2频段绘制效率相同，图 5下方的标示时间与图 4中标示的时间充分证明了这一情况，具体时间有较小波动是因为不同的视角下光线投射穿越数据场的情况不同.当频段继续增加时，只需要按照规则继续增加相应的颜色查找表即可，说明本算法在分频段电磁环境绘制方面具有很好的适应性.
在CUDA下实现并行光线投射算法，绘制效率不仅达到了实时，同时还支持实时改变对数据场的观察视角.用户可以通过鼠标交互调整视角，GPU获取用户输入后将新视角参数传入GPU中改变光线对于数据场的投射方向，从而生成新视角下的绘制图像. 图 6对连续视角变动情况下的绘制时间进行了统计，图 6(a)为视角的变化情况，大致沿一个圆周旋转两圈对数据场进行观察；图 6(b)是对视角变化时连续146帧图像绘制生成时间的统计.
|
图 6 不同视角下绘制时间的统计 |
图 6(b)中的绘制时间最大不超过27.5 ms，绘制速度大于36帧/s，不仅达到了实时绘制的要求，还为实时绘制时进行更多操作留有一定的时间宽裕度(大于12.5 ms).说明本算法不仅实现了对电磁环境的分频段实时绘制，较高的绘制效率还为绘制结果的实时应用提供了支持.同时，如此高效的绘制是因为在CUDA下的并行实现，经实验在CPU中串行实现本算法时，绘制图 6所示的146帧图像的每帧用时均在1 s以上，远不能达到实时绘制的要求.说明本算法的并行执行相对于串行执行的加速比很高，算法适于并行执行.
5 结束语针对电磁环境绘制研究中未考虑频段信息的情况，提出了一种分频段显示绘制算法.通过对颜色查找表的特殊设计和数据频段的划分，基于CUDA实现了分频段的实时光线投射算法，将电磁环境中的频段分布情况绘制到图像上.绘制结果在区分频段的同时还显示了所有频段共同作用的结果，丰富了分频段绘制的内容，能更好地帮助用户理解和掌握电磁环境的态势分布.
但是，本研究只侧重表现了电磁环境的频段信息，并未严格地从空域、时域、频域和能域等方面来讨论电磁环境.分析更多的电磁环境特征并将其可视化是后续研究的主要方向.
| [1] |
吴迎年, 张霖, 张利芳, 等. 电磁环境仿真与可视化研究综述[J]. 系统仿真学报, 2009, 21(20): 6332–6338.
Wu Yingnian, Zhang Lin, Zhang Lifang, et al. Survey on electromagnetic environment simulation and visualization[J].Journal of System Simulation, 2009, 21(20): 6332–6338. |
| [2] |
杨超, 徐江斌, 赵健, 等. 基于多层等值面的电磁环境三维可视化研究[J]. 系统工程与电子技术, 2009, 31(11): 2767–2772.
Yang Chao, Xu Jiangbin, Zhao Jian, et al. Research on 3D visualization of electromagnetic environment based on multi-isosurface[J].Systems Engineering and Electronics, 2009, 31(11): 2767–2772. |
| [3] | Yang Chao, Wu Lingda. GPU-based volume rendering for 3D electromagnetic environment on virtual globe[J].I. J. Image, Graphics and Signal Processing, 2010, 1: 53–60. |
| [4] |
杨超, 徐江斌, 吴玲达. 硬件加速的虚拟电磁环境体可视化[J]. 北京邮电大学学报, 2011, 34(1): 55–59.
Yang Chao, Xu Jiangbin, Wu Lingda. Hardware accelerated volume visualization in virtual electromagnetic environment[J].Journal of Beijing University of Posts and Telecommunications, 2011, 34(1): 55–59. |
| [5] |
马千里, 李思昆, 白晓征, 等. CFD非结构化网格格心格式数据高质量体绘制方法[J]. 计算机学报, 2011, 34(3): 508–516.
Ma Qianli, Li Sikun, Bai Xiaozheng, et al. High-quality volume rendering of unstructured-grid cell-centered data in CFD[J].Chinese Journal of Computers, 2011, 34(3): 508–516. |
| [6] | Schubert N, Scholl I. Comparing GPU-based multi-volume ray casting techniques[J].Comput Sci Res Dev, 2011, 26(1): 39–50. |
| [7] |
郭翰琦, 袁晓如. 体数据可视化传递函数研究[J]. 计算机辅助设计与图形学学报, 2012, 24(10): 1249–1258.
Guo Hanqi, Yuan Xiaoru. Survey on transfer functions in volume visualization[J].Journal of Computer-Aided Design and Computer Graphics, 2012, 24(10): 1249–1258. doi: 10.3969/j.issn.1003-9775.2012.10.001 |
| [8] | Moloney B, Ament M, Weiskopf D, et al. Sort first parallel volume rendering[J].IEEE Transactions on Visualization and Computer Graphics, 2011, 17(8): 1164–1177. doi: 10.1109/TVCG.2010.116 |