集群计算机模拟电磁波传播的方法 【技术领域】
本发明涉及电磁波传播，特别是一种集群计算机模拟电磁波传播的方法。
背景技术
利用数值方法对电磁波传播进行精确模拟，对了解电磁波与各种复杂媒质作用的规律和缩短电磁产品的开发周期至关重要。现有的模拟方法主要有矩量法，电磁场时域有限差分方法以及电磁场时域伪谱方法等。
在上面所提及的各种方法中，电磁场时域伪谱方法于1997年由Q.H.Liu提出(Microwave and optical technology letters，卷：15，页：158-165，1997)，比较适用于大规模电磁波传播的模拟。这种方法首先选取合适的网格长度将电磁波传播的空间划分为三维网格，为模拟电磁波在无限空间的传播，三维网格的最外面nabsorb层为电磁波吸收层。空间中电场和磁场位于各网格单元的顶点，然后通过对麦克斯韦方程组的时间迭代得出时间t后电磁场在空间的分布。最后通过近场外推远场程序获得远场散射场的分布。在模拟过程中，麦克斯韦方程组中的空间导数由一种高精度的傅立叶伪谱方法计算，根据奈奎斯特定理，在这种方法中空间网格长度通常只需略小于电磁波在空间各媒质中传播时最短波长的一半，因而可极大地减少模拟计算中变量的数目。
然而现有的电磁场时域伪谱方法在应用时存在以下几个方面的缺陷：
(1)该方法中电场和磁场在空间位置上重合，在这一构造方法下所采用的空间导数算符作用于不连续函数时会产生吉布斯现象。因而现有方法的构造仅能很好地模拟扩展源激励的电磁波传播问题。对于点源激励的电磁波传播模拟，通常需采用空间扩展源来近似代替点源激励，这一做法不能确切反映真实模拟环境。
(2)在很多电磁波传播模拟中(例如光波在生物组织中传播的模拟)，变量的数目相当巨大，远远超出单个处理器的内存容量，而且模拟时间也会超出可接受的范围。为了有效地解决这个问题，通常可以利用集群计算机进行并行处理，将整个计算区域分解到集群计算机的各个计算节点，然后各个计算节点分别计算各自所在区域内电磁场的传播情况。然而由于电磁场时域伪谱方法中所用的空间导数算符是一个全局算符，它需要利用到整个区域内的电磁场数据，这就不可避免地在整个集群计算机的各个计算节点之间交换大量数据。这一事实使得现有的电磁场时域伪谱方法无法在分布式内存集群计算机上高效地并行处理以模拟大规模电磁波传播。
(3)电磁场时域伪谱方法的时间微分计算采用二阶精度的差分格式，当时间步长不是足够小时会引起较大的数值色散误差。这一点对于大规模的电磁场模拟问题尤为严重，因为数值色散误差带来的相位误差会随着电磁波传播距离的增大而不断积累。通常需在计算中需选取远小于稳定性条件的时间步长，这无疑增加了模拟所需要的迭代步数。
(4)电磁场时域伪谱方法中所采用的格点粗糙，因而在获取远场电磁场分布的近场外推远场程序中对等效面电流在封闭面上进行数值积分时误差较大。
因此，为了适应不同的电磁波传播问题的模拟，必须对现有的方法进行改进，特别需要一种高效可行的并行方案以解决大规模电磁波传播问题。
【发明内容】
本发明的主要目的是提供一种集群计算机模拟电磁波传播的方法，使其能高效地模拟大规模电磁波传播问题。本发明还旨在拓展电磁场时域伪谱方法的适用范围，提高远场计算结果的精度。
本发明的技术方案如下：
一种集群计算机模拟电磁波传播的方法，其特征在于，该方法包括如下步骤：
(1)根据需要模拟的电磁波传播问题建立模型并创建模型数据文件：
①需要模拟的电磁波传播问题为某种形式的电磁波入射到空间区域Ω后在该区域的传播，区域Ω之外为某种单一媒质；任意选取三个相互正交的方向和满足和存在一个包含区域Ω的最小的长方体A，A的各边平行于或以A地中心为原点，和分别为x轴、y轴和z轴的正方向建立三维笛卡尔坐标系；
②选取网格单元长度Δx、Δy和Δz；
③根据所选的网格单元长度将空间划分为三维网格，在x、y和z方向网格单元的数目分别为N′x、N′y和N′z，保证长方体A所占空间都在三维网格内；确定符合下列式子且符合能使快速傅立叶变换算法高效的特征长度的最小正整数Lx、Ly和Lz：
Lx≥((N′x+2nabsorb+20)+2noverlap(nx-1))/nx
Ly≥((N′y+2nabsorb+20)+2noverlap(ny-1))/ny
Lz≥((N′z+2nabsorb+20)+2noverlap(nz-1))/nz
其中noverlap为模型中重叠区域的厚度，nabsorb为电磁波吸收层的厚度，nx、ny和nz分别为集群计算机三维拓扑结构在x、y和z方向计算节点的数目；计算
Nx＝2(Lx-noverlap)+(nx-2)(Lx-2noverlap)
Ny＝2(Ly-noverlap)+(ny-2)(Ly-2noverlap)
Nz＝2(Lz-noverlap)+(nz-2)(Lz-2noverlap)；
④将已有的三维网格在x轴正方向和x轴负方向分别向外延伸nx1和nx2层网格单元，得到一个新的三维网格，x、y和z方向网格单元的数目分别为Nx、N′y和N′z，其中nx2＝Nx-N′x-nx1，函数ceil(s)取大于等于自变量s的最小正整数；
⑤将第④步得到的新的三维网格在y轴正方向和y轴负方向分别向外延伸ny1和ny2层网格单元，得到一个新的三维网格，x、y和z方向网格单元的数目分别为Nx、Ny和N′z，其中ny2＝Ny-N′y-ny1；
⑥将第⑤步得到的新的三维网格在z轴正方向和z轴负方向分别向外延伸nz1和nz2层网格单元，得到最终的三维网格，x、y和z方向网格单元的数目分别为Nx、Ny和Nz，其中nz2＝Nz-N′z-nz1；该三维网格所占据空间即为电磁波传播模拟的整个空间，最外面nabsorb层网格单元为电磁波吸收层；
⑦按照Yee网格的方法配置电场和磁场在三维网格各网格单元中的位置；去除位于三维网格各方向的正方向一侧端面上的所有电磁场分量，网格中各电磁场分量分别对应一个三维矩阵，这些矩阵具有相同的形状：x方向的长度均为Nx，y方向的长度均为Ny，z方向的长度均为Nz；
⑧利用重叠区域分解方法将模拟的整个空间划分为块，x、y和z方向的块的数目分别等于集群计算机三维拓扑结构中该方向计算节点的数目，各块中各电磁场分量的范围由下面的方法确定：将各电磁场分量组成的矩阵分别分割成连续的子矩阵，各方向子矩阵的数目分别等于该方向的块的数目，各子矩阵在方向ξ上元素的数目为Lξ-noverlap或Lξ-2noverlap，当该子矩阵为ξ方向的首或尾子矩阵时为Lξ-noverlap，其余情况为Lξ-2noverlap，其中ξ为x、y或z；将所得的各个子矩阵的区域延伸，原则是若该子矩阵在某个边界有相邻子矩阵且该子矩阵中对应的电磁场分量方向与该边界相切，则将该子矩阵在此边界处与此边界垂直方向的坐标向外延伸noverlap；上述各个电磁场分量的各个子矩阵一一对应于各个块，且延伸之前的区域为该块中的非重叠区域，后来延伸的区域为重叠区域；
⑨选取时间步长Δt为其中T为电磁波的周期，vmax为电磁波在空间各媒质中传播的最大速率；
⑩确定迭代次数N和记录电磁场分布的时间步n1和n2；
将空间中各媒质的介电常数和磁导率分别修改为εh′＝αεh和μh′＝αμh，其中为电磁波的角频率，εh和μh分别为第h类媒质实际的介电常数和磁导率，εh′和μh′分别为在模拟中使用的第h类媒质介电常数和磁导率；
创建模型数据文件：所述的重叠区域分解方法分割得到的每一个块，分别建立模型数据文件，且命名为“model_ib_jb_kb”，其中ib、jb和kb分别为该块在x、y和z方向的序号，将下列内容写入各块所对应的模型数据文件：该块中各电磁场分量在各方向重叠区域和非重叠区域的坐标范围、该块中非重叠区域中各电磁场分量位置处媒质修改后的介电常数和磁导率、网格单元在各方向的长度、模拟中的时间步长、迭代次数、时间步n1和n2、以及入射场的参数；
(2)集群计算机进行数据处理：
读取模型数据文件，用并行电磁场时域伪谱方法并行模拟电磁波的传播，输出记录时间步n1和n2时空间电磁场分布的结果文件；
①集群计算机的各个计算节点分别读取步骤(1)中所得到的对应的模型数据文件，得到模拟所需的参数，集群计算机三维拓扑结构中坐标为(ip，jp，kp)的计算节点对应的模型数据文件为“model_ip_jp_kp”；各个计算节点在内存中分配存储各电磁场分量的空间；
②利用并行电磁场时域伪谱方法模拟空间中的电磁波传播；
③输出结果文件“result1”和“result2”，该结果文件“result1”和“result2”分别记录了第n1步和第n2步时空间各电磁场分量的数值；
(3)利用近场外推远场程序计算远场的散射场分布：
①读取步骤(2)中输出的结果文件“result1”和“result2”，得到近场中各电磁场分量的复数场；
②在三维网格中，设置一个虚拟的包围区域Ω的封闭面，由六个切面x＝x1、x＝x2、y＝y1、y＝y2、z＝z1和z＝z2围成；
③插值得到上述各切面上与该切面法线方向垂直的各磁场分量；
④在各切面上分别对各等效面电流和各等效面磁流作二维傅立叶变换，将各切面上的各等效面电流和各等效面磁流展开为级数和的形式；
⑤解析积分得到外推方向的远场时所需的中间量和
⑥由和计算得到方向的远场散射场和
所述的并行电磁场时域伪谱方法是通过时间迭代模拟电磁波的传播，在将空间中电磁场由当前时间步迭代至下一时间步的过程中各个计算节点依次完成下列四步：
(1)利用电场的时间迭代式求下一时间步本计算节点中非重叠区域的电场；其中电场的时间迭代式中涉及到的各磁场分量的空间导数用本节点内的傅立叶伪谱方法计算得到，即仅针对本计算节点的磁场数据进行操作，包括其重叠区域和非重叠区域，且先在重叠区域的磁场上乘上一个平滑函数使其从与非重叠区域连接的一端到另一端平滑地下降为0；
(2)接收各个相邻计算节点发送过来的重叠区域的电场数据，并向各个相邻计算节点发送自身非重叠区域与该相邻计算节点的重叠区域部分重合的电场数据；
(3)利用磁场的时间迭代式求下一时间步本计算节点中非重叠区域的磁场；其中磁场的时间迭代式中涉及到的各电场分量的空间导数用本节点内的傅立叶伪谱方法计算得到，即仅针对本计算节点的电场数据进行操作，包括其重叠区域和非重叠区域，且先在重叠区域的电场上乘上一个平滑函数使其从与非重叠区域连接的一端到另一端平滑地下降为0；
(4)接收各个相邻计算节点发送过来的重叠区域的磁场数据，并向各个相邻计算节点发送自身非重叠区域与该相邻计算节点的重叠区域部分重合的磁场数据。
本发明的技术效果：
(1)在电磁场时域伪谱方法中，将三维网格中电场和磁场的位置在空间各个方向上错开半个网格长度，基于此构造方法的空间导数算符作用于点源时不会产生吉布斯现象。改进后的构造方法既可用于精确模拟扩展源激励的电磁波传播又可用于精确模拟点源激励的电磁波传播。
(2)采用重叠区域分解方法对电磁场时域伪谱方法进行并行处理。将整个计算区域分割成块并分配到集群计算机的各个计算节点中去，每个计算节点的数据分为重叠区域和非重叠区域。在每一时间步的计算过程中，各个节点仅计算各自非重叠区域的电磁场，而重叠区域的电磁场在每一时间步的计算后通过数据交换由相邻节点获得。尤其重要的是，在计算非重叠区域的电磁场时，仅需利用各计算节点内存中的数据来计算空间电磁场的导数。为了在此过程中不引起误差，将本地数据中的重叠区域部分乘上平滑函数以形成锥形尾端。这一并行处理方法避免了在集群计算机的各计算节点中交换大量数据，数据交换仅限于相邻计算节点之间少量的重叠区域，因而能够高效地模拟大规模电磁波传播。
(3)修改模拟中所采用的空间各媒质的介电常数和磁导率。通过在模拟中修改空间各媒质的介电常数和磁导率，本发明消除了在选取较大的时间步长时所带来的较大的数值色散误差，从而成倍地增大模拟中可采用的时间步长，减少模拟所需的迭代步数，节约计算资源。
(4)对近场外推远场计算中等效面电流和等效面磁流在封闭面上的积分方法进行了改进。首先将外推面上的等效面电流和等效面磁流作二维傅立叶变换获得傅立叶频谱，从而可以将它们展开为级数和的形式，对于空间快速变化的相位延迟部分即可解析积分。这一改进克服了电磁场时域伪谱方法中格点粗糙的缺点，远场计算结果更为精确。
总之，本发明的优点是可用于模拟点源激励的电磁波传播；可高效地并行模拟大规模电磁波传播；减少模拟中迭代的次数；远场结果更为精确。
【附图说明】
图1是本发明空间网格划分示意图
图2是本发明空间网格中每个网格单元上的电场和磁场位置示意图
图3是本发明集群计算机读取模型数据文件、处理模型并输出结果文件的流程示意图
图4是本发明并行电磁场时域伪谱方法中的平滑函数的示意图
图5是本发明并行电磁场时域伪谱方法中各个计算节点求电磁场空间导数的示意图
图6是本发明重叠区域分解方法的示意图
图7是本发明实施例1中平面波入射生物组织模型的示意图
图8是本发明实施例1的流程示意图
图9是本发明实施例1的模拟结果
图10是本发明实施例2的流程示意图
图11是本发明实施例2的模拟结果
【具体实施方式】
下面将结合附图和实施例对本发明作进一步说明，但不应以此限制本发明的保护范围。
本发明中集群计算机模拟电磁波传播的方法，包括三个步骤：
一、根据模拟的电磁波传播问题建立模型并创建模型数据文件：
①需要模拟的电磁波传播问题为某种形式的电磁波入射到空间区域Ω后在该区域的传播，区域Ω之外为某种单一媒质；任意选取三个相互正交的方向和满足和存在一个包含区域Ω的最小的长方体A，A的各边平行于或以A的中心为原点，和分别为x轴、y轴和z轴的正方向建立三维笛卡尔坐标系；
②选取网格单元长度Δx、Δy和Δz；各方向网格单元长度的选择方法与现有电磁场时域伪谱方法一致，必须小于电磁波在空间各媒质中最短波长的1/2，在此基础上不得大于该方向电磁性质不均匀的最小尺度；
③根据所选的网格单元长度将空间划分为三维网格，在x、y和z方向网格单元的数目分别为N′x、N′y和N′z，保证长方体A所占空间都在三维网格内；确定符合下列式子且符合能使快速傅立叶变换算法高效的特征长度的最小正整数Lx、Ly和Lz：
Lx≥((N′x+2nabsorb+20)+2noverlap(nx-1))/nx
Ly≥((N′y+2nabsorb+20)+2noverlap(ny-1))/ny
Lz≥((N′z+2nabsorb+20)+2noverlap(nz-1))/nz
其中noverlap为模型中重叠区域的厚度，nabsorb为吸收层的厚度，nx、ny和nz分别为集群计算机三维拓扑结构在x、y和z方向计算节点的数目，能使快速傅立叶变换高效的特征长度l满足形式l＝2m·3n·5d(m、n和d均为非负整数)；计算
Nx＝2(Lx-noverlap)+(nx-2)(Lx-2noverlap)
Ny＝2(Ly-noverlap)+(ny-2)(Ly-2noverlap)
Nz＝2(Lz-noverlap)+(nz-2)(Lz-2noverlap)；
④将已有的三维网格在x轴正方向和x轴负方向分别向外延伸nx1和nx2层网格单元，得到一个新的三维网格，x、y和z方向网格单元的数目分别为Nx、N′y和N′z，其中nx2＝Nx-N′x-nx1，函数ceil(s)取大于或等于自变量s的最小正整数；
⑤将第④步得到的新的三维网格在y轴正方向和y轴负方向分别向外延伸ny1和ny2层网格单元，得到一个新的三维网格，x、y和z方向网格单元的数目分别为Nx、Ny和N′z，其中ny2＝Ny-N′y-ny1；
⑥将第⑤步得到的新的三维网格在z轴正方向和z轴负方向分别向外延伸nz1和nz2层网格单元，得到如图1所示的最终的三维网格，x、y和z方向网格单元的数目分别为Nx、Ny和Nz，其中nz2＝Nz-N′z-nz1；该三维网格所占据空间即为电磁波传播模拟的整个空间，最外面nabsorb层网格单元为电磁波吸收层；
⑦按照Yee网格的方法配置电场和磁场在三维网格中的位置，如图2所示：电场位于网格单元的各棱边的中点，电场方向沿着所在棱边的正方向；磁场位于各网格单元面的中心，磁场方向沿着所在面的垂线的正方向；去除位于三维网格各方向的正方向一侧端面上的所有电磁场分量，网格中各电磁场分量分别对应一个三维矩阵，这些矩阵具有相同的形状，x方向的长度均为Nx，y方向的长度均为Ny，z方向的长度均为Nz；
⑧利用重叠区域分解方法将模拟的整个空间划分为块，x、y和z方向的块的数目分别等于集群计算机三维拓扑结构中该方向计算节点的数目，各块中各电磁场分量的范围由下面的方法确定：将各电磁场分量组成的矩阵分别分割成连续的子矩阵，各方向子矩阵的数目分别等于该方向上块的数目，各子矩阵在方向ξ上元素的数目为Lξ-noverlap或Lξ-2noverlap，当该子矩阵为ξ方向的首或尾子矩阵时为前面一个值，其余情况为后面一个值，其中ξ为x、y或z；将所得的各个子矩阵的区域延伸，方法是若该子矩阵在某个边界有相邻子矩阵且该子矩阵中对应的电磁场分量方向与该边界相切，则将该子矩阵在此边界处与此边界垂直方向的坐标向外延伸noverlap；上述得到的各电磁场分量的各个子矩阵一一对应于各个块，且延伸之前的区域为该块中的非重叠区域，后来延伸的区域为重叠区域；由上述方法可给出各块中各电磁场分量的坐标范围，以Ex为例给出表达式：由于Ex沿着x方向，因此在x方向上Ex没有重叠，而在y和z方向与相邻计算节点之间存在重叠。x方向上第i个块中Ex的数据区间[ai，bi]可由下式决定：
[ai,bi]=[1,Lx-noverlap]i=1[bi-1+1,ai+Lx-2noverlap-1]i=2,...,nx-1[bi-1+1,Nx]i=nx]]>
y方向上第j块中Ex的数据区间[cj，dj]由下式决定
[cj,dj]=[1,Ly]j=1[dj-1-2noverlap+1,cj+Ly-1]j=2,...,ny-1[dj-1-2noverlap+1,Ny]j=ny]]>
z方向上第k块中Ex的数据区间[ek，fk]由下式决定
[ek,fk]=[1,Lz]k=1[fk-1-2noverlap+1,ek+Lz-1]k=2,...,nz-1[fk-1-2noverlap+1,Nz]k=nz]]>
其余各电磁场分量的坐标范围与此类似；
⑨选取时间步长Δt为其中T为电磁波的周期，vmax为电磁波在空间各媒质中传播的最大速率；
⑩确定模拟中的迭代次数其中tsimulation为模拟中电磁波传播的时间，也即在空间中电磁场分布达到稳态所需要的时间，选取为电磁波往返整个模拟空间所需时间的两倍；
将记录空间中电磁场分布的时间步n1和n2分别设置为，n1＝N-m和n2＝N；
将空间中各媒质的介电常数和磁导率分别修改为εh′＝αεh和μh′＝αμh，其中为电磁波的角频率，εh和μh分别为第h类媒质实际的介电常数和磁导率，εh′和μh′分别为在模拟中使用的第h类媒质介电常数和磁导率；
创建模型数据文件：对于重叠区域分解方法分割得到的每一个小块，建立一个模型数据文件，且命名为“model_ib_jb_kb”，其中ib、jb和kb分别为该块在x、y和z方向的序号；将下列内容写入各块所对应的模型数据文件：该块中各电磁场分量在各方向重叠区域和非重叠区域的坐标范围、该块中非重叠区域中各电磁场分量位置处媒质修改后的介电常数和磁导率、网格单元在各方向的长度、模拟中的时间步长、迭代次数、时间步n1和n2以及入射波的参数；入射波的参数包括入射波形状、方向和频率等确定入射波在空间中任何时刻任意点的电磁场各分量值所需的所有参数；
二、集群计算机进行数据处理：
读取模型数据文件，用并行电磁场时域伪谱方法并行模拟电磁波的传播，输出记录时间步n1和n2时空间电磁场分布的结果文件；
图3展示了该步骤的流程图：
(1)调用消息传递接口程序，启动MPI(Message Passing Interface)进程；
启动p个进程(p＝nx×ny×nz)，分别对应于集群计算机中的p个计算节点；将这p个进程映射至三维拓扑结构，x、y和z方向的数目分别为nx、ny和nz；
(2)各个计算节点读取所对应的模型文件；
各个计算节点分别读取对应的模型文件(假设某个计算节点在三维拓扑结构中的坐标为(ip，jp，kp)，则其所对应的模型文件为“model_ip_jp_kp”)，得到如下参数：网格单元各方向的长度Δx、Δy和Δz，模拟中的时间步长Δt、总的迭代次数N、记录空间中电磁场分布的时间步n1和n2、各电磁场分量非重叠区域以及重叠区域的坐标范围、各电磁场分量在非重叠区域各位置处的电磁参数以及入射波的参数；分配存储各电磁场分量的内存空间；
(3)开始模拟，设置时间步n＝1；
(4)各计算节点利用电场的时间迭代式计算非重叠区域中下一步的电场；
电场的时间迭代式分为两个区域：
非吸收层，
Ex|i+1/2,j,kn+1=Ex|i+1/2,j,kn+Δtϵ′|i+1/2,j,k[(∂Hz∂y-∂Hy∂z)|i+1/2,j,kn+1/2-(ϵ′|i+1/2,j,k-ϵ0′)·(∂Ex,inc∂t)|i+1/2,j,kn+1/2]Ey|i,j+1/2,kn+1=Ey|i,j+1/2,kn+Δtϵ′|i,j+1/2,k[(∂Hx∂z-∂Hz∂x)|i,j+1/2,kn+1/2-(ϵ′|i,j+1/2,k-ϵ0′)·(∂Ey,inc∂t)|i,j+1/2,kn+1/2]Ez|i,j,k+1/2n+1=Ez|i,j,k+1/2n+Δtϵ′|i,j,k+1/2[(∂Hy∂x-∂Hx∂y)|i,j,k+1/2n+1/2-(ϵ′|i,j,k+1/2-ϵ0′)·(∂Ez,inc∂t)|i,j,k+1/2n+1/2]]]>
吸收层，
Px|i+1/2,j,kn+1=2-Δt·by|j2+Δt·by|jPx|i+1/2,j,kn+2Δt(2+Δt·by|j)ϵ0′(∂Hz∂y-∂Hy∂z)|i+1/2,j,kn+1/2Py|i,j+1/2,kn+1=2-Δt·bz|k2+Δt·bz|kPy|i,j+1/2,kn+2Δt(2+Δt·bz|k)ϵ0′(∂Hx∂z-∂Hz∂x)|i,j+1/2,kn+1/2Pz|i,j,k+1/2n+1=2-Δt·bx|i2+Δt·bx|iPz|i,j,k+1/2n+2Δt(2+Δt·bx|i)ϵ0′(∂Hy∂x-∂Hx∂y)|i,j,k+1/2n+1/2]]>
Ex|i+1/2,j,kn+1=2-Δt·bz|k2+Δt·bz|kEx|i+1/2,j,kn+2+Δt·bx|i+1/22+Δt·bz|kPx|i+1/2,j,kn+1-2-Δt·bx|i+1/22+Δt·bz|kPx|i+1/2,j,knEy|i,j+1/2,kn+1=2-Δt·bx|i2+Δt·bx|iEy|i,j+1/2,kn+2+Δt·by|j+1/22+Δt·bx|iPy|i,j+1/2,kn+1-2-Δt·by|j+1/22+Δt·bx|iPy|i,j+1/2,knEz|i,j,k+1/2n+1=2-Δt·by|j2+Δt·by|jEz|i,j,k+1/2n+2+Δt·bz|k+1/22+Δt·by|jPz|i,j,k+1/2n+1-2-Δt·bz|k+1/22+Δt·by|jPz|i,j,k+1/2n]]>
Einc为入射场，Px，Py，Pz为中间变量，bx，by，bz为吸收层的吸收参数可由下式确定：
bx|i=bxmax·(nabsorb-inabsorb)3i≤nabsorb0nabsorb<i<Nx-nabsorbbxmax·(i-(Nx-nabsorb)nabsorb)3i≥Nx-nabsorb]]>
by|j=bymax·(nabsorb-jnabsorb)3j≤nabsorb0nabsorb<j<Ny-nabsorbbymax·(j-(Ny-nabsorb)nabsorb)3j≥Ny-nabsorb]]>
bz|k=bzmax·(nabsorb-knabsorb)3k≤nabsorb0nabsorb<k<Nz-nabsorbbzmax·(k-(Nz-nabsorb)nabsorb)3k≥Nz-nabsorb]]>
其中c为真空中的光速，在本发明的实施例中选取nabsorb＝8；
上述时间迭代式中涉及到的各磁场分量的空间导数用本节点内的傅立叶伪谱方法计算得到，即仅针对本计算节点的磁场数据进行操作，包括磁场的重叠区域和非重叠区域。在用傅立叶伪谱方法求磁场的空间导数之前先将重叠区域的磁场乘以平滑函数w(l)，其中1≤l≤noverlap，w(l)前半段的值为1，后半段的值单调光滑地从1下降至0(定性见图4)，要求如此处理后的磁场能满足下面三个条件：①非重叠区域内的数据不变；②本节点上的数据满足周期性边界条件，即傅立叶变换操作的数据在左右两个端点为0；③磁场数据的一阶导数跨分界面连续。多次数值试验表明结果的精度对平滑函数的精确形式不敏感，仅须满足函数上端光滑连接到1，下端光滑下降到0。在本发明的实施例中选取平滑函数为(为程序编写方便，设noverlap为偶数)
w(l)=11≤l≤noverlap/2cos[π(l-noverlap/2)/noverlap]noverlap/2<l≤noverlap;]]>
重叠区域的第l个磁场数据乘上函数值w(l)，使重叠区域的磁场数据从与非重叠区域连接的一端到另一端光滑地下降为0。例如图5为某个计算节点中方向ξ(ξ可为x、y或z)上求空间导数的示意图，箭头代表磁场的位置，圆点代表电场的位置，电场位置处磁场的空间导数可按如下方法计算：
在本发明的实施例中进一步选取noverlap＝10，平滑函数的形式简化为
w(l)=11≤l≤5cos[π(l-5)/10]6≤l≤10]]>
将平滑函数作用到磁场的重叠区域上，作用后的磁场为
H~(i)=H(i)·w(i1+10-i)i1≤i≤i1+9H(i)i1+10≤i≤i2-10H(i)·w(i-i2+10)i2-9≤i≤i2,]]>
利用下面表达式求得磁场在电场位置处的空间导数
∂H∂ξ=F-1[jkξejkξΔξ/2F(H~)],]]>
其中j为虚数单位即j2＝-1；
(5)各计算节点接收各个相邻计算节点发送过来的重叠区域的电场数据，并向各个相邻计算节点发送自身非重叠区域与该相邻计算节点的重叠区域部分重合的电场数据；
(6)各计算节点利用磁场的时间迭代式计算非重叠区域中下一步的磁场；
磁场的时间迭代式分为两个区域：
非吸收层，
Hx|i,j+1/2,k+1/2n+3/2=Hx|i,j+1/2,k+1/2n+1/2+Δtμ′|i,j+1/2,k+1/2[(∂Ey∂z-∂Ez∂y)|i,j+1/2,k+1/2n+1-(μ′|i,j+1/2,k+1/2-μ0′)·(∂Hx,inc∂t)|i,j+1/2,k+1/2n+1]Hy|i+1/2,j,k+1/2n+3/2=Hy|i+1/2,j,k+1/2n+1/2+Δtμ′|i+1/2,j,k+1/2[(∂Ez∂x-∂Ex∂z)|i+1/2,j,k+1/2n+1-(μ′|i+1/2,j,k+1/2-μ0′)·(∂Hy,inc∂t)|i+1/2,j,k+1/2n+1]Hz|i+1/2,j+1/2,kn+3/2=Hz|i+1/2,j+1/2,kn+1/2+Δtμ′|i+1/2,j+1/2,k[(∂Ex∂y-∂Ey∂x)|i+1/2,j+1/2,kn+1-(μ′|i+1/2,j+1/2,k-μ0′)·(∂Hz,inc∂t)|i+1/2,j+1/2,kn+1]]]>
吸收层，
Bx|i,j+1/2,k+1/2n+3/2=2-Δt·by|j+1/22+Δt·by|j+1/2Bx|i,j+1/2,k+1/2n+1/2+2Δt(2+Δt·by|j+1/2)μ0′(∂Ey∂z-∂Ez∂y)|i,j+1/2,k+1/2n+1By|i+1/2,j,k+1/2n+3/2=2-Δt·bz|k+1/22+Δt·bz|k+1/2By|i+1/2,j,k+1/2n+1/2+2Δt(2+Δt·bz|k+1/2)μ0′(∂Ez∂x-∂Ex∂z)|i+1/2,j,k+1/2n+1Bz|i+1/2,j+1/2,kn+3/2=2-Δt·bx|i+1/22+Δt·bx|i+1/2Bz|i+1/2,j+1/2,kn+1/2+2Δt(2+Δt·bx|i+1/2)μ0′(∂Ex∂y-∂Ey∂x)|i+1/2,j+1/2,kn+1]]>
Hx|i,j+1/2,k+1/2n+3/2=2-Δt·bz|k+1/22+Δt·bz|k+1/2Hx|i,j+1/2,k+1/2n+1/2+2+Δt·bx|i2+Δt·bz|k+1/2Bx|i,j+1/2,k+1/2n+3/2-2-Δt·bx|i2+Δt·bz|k+1/2Bx|i,j+1/2,k+1/2n+1/2Hy|i+1/2,i,k+1/2n+3/2=2-Δt·bx|i+1/22+Δt·bx|i+1/2Hy|i+1/2,j,k+1/2n+1/2+2+Δt·by|j2+Δt·bx|i+1/2By|i+1/2,j,k+1/2n+3/2-2-Δt·by|j2+Δt·bx|i+1/2By|i+1/2,j,k+1/2n+1/2Hz|i+1/2,j+1/2,kn+3/2=2-Δt·by|j+1/22+Δt·by|j+1/2Hz|i+1/2,j+1/2,kn+1/2+2+Δt·bz|k2+Δt·by|j+1/2Bz|i+1/2,j+1/2,kn+3/2-2-Δt·bz|k2+Δt·by|j+1/2Bz|i+1/2,j+1/2,kn+1/2]]>
Hinc为入射场分布，Bx，By，Bz为中间变量，bx，by，bz与上述(4)中参数相同；
上述时间迭代式中涉及到的各电场分量的空间导数用本节点内的傅立叶伪谱方法计算得到，即仅针对本计算节点的电场数据进行操作，包括电场的重叠区域和非重叠区域。在用傅立叶伪谱方法求电场的空间导数之前先将重叠区域的电场乘以平滑函数，为保证电场和磁场空间导数计算的一致性，这里采用与上述(4)中相同的平滑函数w(l)。重叠区域的第l个电场数据乘上函数值w(l)，使重叠区域的电场数据从与非重叠区域连接的一端到另一端光滑地下降为0。图5中磁场位置处电场的空间导数可按如下方法计算：
将平滑函数作用到电场的重叠区域上，作用后的电场为
E~(i)=E(i)·w(i1+10-i)i1≤i≤i1+9E(i)i1+10≤i≤i2-10E(i)·w(i-i2+10)i2-9≤i≤i2,]]>
利用下面表达式求得电场在磁场位置处的空间导数
∂E∂ξ=F-1[jkξejkξΔξ/2F(E~)];]]>
(7)接收各个相邻计算节点发送过来的重叠区域的磁场数据，并向各个相邻计算节点发送自身非重叠区域与该相邻计算节点的重叠区域部分重合的磁场数据；
(8)判断时间步n是否为n1或n2：若是，则记录此时各计算节点内非重叠区域内的电磁场；若否，则进入步骤(9)；
(9)判断n＜N是否成立：若是，则n＝n+1并返回步骤(4)；若否，则进入步骤(10)；
(10)除了第一个计算节点外的其余计算节点将第n1步和第n2步时记录的电磁场发送给第一个计算节点；第一个计算节点接收从其余各计算节点中发送过来的数据，并与其自身记录的数据一起合并成整个模拟空间内的电磁场分布，将第n1步空间的电磁场分布写入输出文件“result1”，将第n2步空间的电磁场分布写入输出文件“result2”；
(11)终止所述的MPI进程。
三、利用近场外推远场程序计算远场的散射场分布；
①读取集群计算机的输出文件“result1”和“result2”，得到空间中近场电磁场在第n1步和第n2步时的分布和
近场中散射场的复数形式由下式获得：
E→=E→n1-jE→n2]]>
H→=H→n1-jH→n2;]]>
②根据标准的近场外推远场程序的要求，设置一个长方形的封闭面(由六个切面x＝x1，x＝x2，y＝y1，y＝y2，z＝z1，z＝z2围成，且能包围住区域Ω)；
③插值得到上述各切面上与该切面法线方向垂直的磁场；如为获得面x＝x1上y方向的磁场分量Hy，首先对位于面x＝x1法线上的所有Hy作一维傅立叶变换得到其傅立叶频谱，然后在各频谱分量上乘上移位因子ejkΔx/2，k为对应的离散波矢，然后对频谱作反傅立叶变换的结果，即
Hysca|x=x1=F-1(ejkΔx/2F(Hysca))|x=x1-Δx/2;]]>
④切面上的等效面电流和等效面磁流分别定义为和其中为是该切面上向外的单位法向矢量；分别对各切面上各等效面电流分量和各等效面磁流分量作二维傅立叶变换，根据得到的频谱将它们在所在切面上的分布展开成级数和的形式；例如在z＝z1面上等效面电流x方向的分量为
Jx(x,y,z1)=Σm,nJx~(m,n)ej(kx,mx+ky,ny)]]>
为对整个z＝z1切面上的等效面电流Jx作二维傅立叶变换时的频谱分量，kx，m和ky，n分别为对应的离散波矢；
⑤为获取方向的远场需先计算中间量和其中
N→(r^)=∫∫s(Jxx^+Jyy^+Jzz^)ejkr′cosψds′L→(r^)=∫∫s(Mxx^+Myy^+Mzz^)ejkr′cosψds′,]]>
ejkr′cosψ为向方向传播时积分面元处的相对相位延迟量，由于已经在各面上将等效面电流和等效面磁流展开成明确的级数和的形式，上式可解析积分；例如对z＝z1面上的x方向分量Jx积分为
∫x1x2∫y1y2Jx(x,y,z1)ej(kxx+kyy+kzz1)dxdy]]>
=∫x1x2∫y1y2Σm,nJx~(m,n)ej(kx,mx+ky,ny)ej(kxx+kyy+kzz1)dxdy]]>
=ejkzz1Σm,n(Jx~(m,n)ej(kx,m+kx)x2-ej(kx,m+kx)x1j(kx,m+kx)ej(ky,n+ky)y2-ej(ky,n+ky)y1j(ky,n+ky))]]>
其中kx、ky和kz分别为电磁波向方向传播时的波矢在x、y和z方向的分量。
⑥最后由下式计算得到方向的远场散射场
Er=0Eθ=-jke-jkr4πr(Lφ+η0Nθ)Eφ=+jke-jkr4πr(Lθ-η0Nφ)]]>
Hr=0Hθ=+jke-jkr4πr(Nφ-Lθη0)Hφ=-jke-jkr4πr(Nθ+Lφη0).]]>
为进一步阐释本发明重叠区域分解并行处理电磁场时域伪谱方法的思想，在图6中给出了重叠区域分解方法的示意图，以ξ方向区域的分解为例作说明。如图6-1所示，区域在ξ轴上的区间为[a，k]，宽度为Nξ。现要将其划分到集群计算机的4个计算节点(计算节点A、计算节点B、计算节点C和计算节点D)，每个节点中数据的宽度为Lξ，重叠区域厚度为noverlap。如图6-2所示，首先将[a，k]分割成10段[a，b]、[b，c]、[c，d]、[d，e]、[e，f]、[f，g]、[g，h]、[h，i]、[i，j]和[j，k]，各段宽度分别为Lξ-2noverlap、noverlap、noverlap、Lξ-4noverlap、noverlap、noverlap、Lξ-4noverlap、noverlap、noverlap和Lξ-2noverlap。然后如图6-3所示将区域[a，d]、[b，g]、[e，j]和[h，k]依次分配至计算节点A、计算节点B、计算节点C和计算节点D。在计算节点A中，区域[a，c]的电磁场用电磁场时域伪谱方法计算，区域[c，d]接收由计算节点B发送的对应段数据；在计算节点B中，区域[c，f]的电磁场用电磁场时域伪谱方法计算，区域[b，c]接收由计算节点A发送的对应段数据，区域[f，g]接收由计算节点C发送的对应段数据；在计算节点C中，区域[f，i]的电磁场用电磁场时域伪谱方法计算，区域[e，f]接收计算节点B发送的对应段数据，区域[i，j]接收计算节点D发送的对应段数据；在计算节点D中，区域[i，k]的电磁场用电磁场时域伪谱方法计算，区域[h，i]接收计算节点C发送的对应段数据。图6-4依次给出了各个计算节点中求电磁场空间导数时在自身数据上所乘因子的示意图：在非重叠区因子恒为1，电磁场数据不变；在重叠区域因子由非重叠区域和重叠区域边界向外平滑地降为0，作用后的电磁场数据也平滑地下降为0。
实施例1：
本实施例研究偏振光入射生物组织模型时后向散射光的偏振特性，从而为光学检测生物组织中的偏振门控技术提供参数。生物组织模型采用离散散射体模型，如图7所示，即在一定空间体积范围(lx×ly×lz)内随机分布大小均一的介质微球。微球的半径为1μm，折射率为1.59；微球以外空间是水，折射率为1.33。入射光为沿Z方向输入的x方向线偏振的平面波，其在自由空间中的波长为785nm。
在本实施例中共模拟了四组不同尺寸的组织模型，各组模型的空间尺寸、模型中微球的数目、在集群计算机模拟时网格的宽度及各计算节点中数据长度列于下表：
模 型 微球 数目 lx (μm) ly (μm) lz (μm) Nx Ny Nz Lx Ly Lz 1 1735 100 100 25 1092 1092 300 288 288 160 2 3470 100 100 50 1092 1092 556 288 288 288 3 5205 100 100 75 1092 1092 876 288 288 448 4 6940 100 100 100 1092 1092 1132 288 288 576
此外在各组模型中选取了相同的网格单元长度(Δx＝Δy＝Δz＝0.098μm)，所有的模拟利用集群计算机的32个计算节点，组成4×4×2的三维拓扑结构。
本实施例1中网格包括最内层的生物组织模型、最外层的吸收层以及介于上述两者的中间层。中间层的厚度可调节以使参数Lx、Ly和Lz的大小符合快速傅立叶变换算法的特征长度。
图8是本实施例1的流程图：
(1)设置i＝1；
(2)设置j＝1；
(3)生成第i组模型的第j个实现：在第i组模型限定的空间范围内随机产生微球位置，保证所有的微球互不重叠且在空间中分布均匀；
(4)建立第i组模型的第j个实现的模型并生成模型文件；
(5)将步骤(4)中得到的模型输入给集群计算机的各个计算节点；集群计算机利用改进了的电磁场时域伪谱方法模拟平面波在该生物组织模型中的传播情况；集群计算机输出记录空间中电磁场分布稳定后的电磁场分布的数据文件；
(6)读取步骤(5)得到的数据文件，用近场外推远场程序计算得到远场散射场分布和
(7)分别计算散射角度在175°≤θ≤180°，0°≤φ≤360°内与入射光偏振方向平行的总光强和与入射光偏振方向垂直的总光强，将两者相减得到偏振差分信号ΔIi，j：
将步骤(6)中得到的后向散射场的电场分解成与入射场的电场平行和垂直的两部分，磁场分解成与入射场的磁场平行和垂直的两部分
E→i,j(θ,φ)=E→||,i,j(θ,φ)+E→⊥,i,j(θ,φ)H→i,j(θ,φ)=H→||,i,j(θ,φ)+H→⊥,i,j(θ,φ)]]>
从而得到
ΔIi,j=12∫∫(|E→||,i,j(θ,φ)×H→*||,i,j(θ,φ)|-|E→⊥,i,j(θ,φ)×H→*⊥,i,j(θ,φ)|)sinθdθdφ.]]>
(8)判断j＜5是否成立：若是，j＝j+1并返回步骤(3)；若否，则进入步骤(9)；
(9)计算第i组模型的平均差分信号。
ΔIi=15Σj=15ΔIi,j]]>
(10)判断i＜4是否成立：若是，i＝i+1并返回步骤(2)；若否，结束整个模拟。
图9为本实施例的模拟结果，横坐标为各组生物组织模型对应的光学厚度，由Mie散射理论计算得到，纵坐标为所对应的模型后向散射差分信号强度。从图9可以看出，差分信号在光学厚度大于4时已经饱和。本实施例说明了通过检测后向散射差分信号可以将检测范围限定在生物组织的几个散射长度的表层组织内。
实施例2：
本实施例2模拟计算光波入射生物组织模型时的背向散射增强现象。本实施例中的生物组织模型构造方法与实施例1中的生物组织模型构造方法一致，且各参数与实施例1中的第2组模型相同。本实施例为抑制相干光与生物组织模型产生的散斑，对不同频率的光入射的结果作平均。
图10是本实施例2的流程图：
(1)生成组织模型：在限定的空间范围内随机产生微球位置，保证所有的微球互不重叠且在空间中分布均匀；
(2)设置i＝1；
(3)计算第i个模拟的入射光频率：fi＝(0.95+0.005(i-1))f0，其中f0为自由空间中波长为785nm光波的波长；
(4)建立生物组织模型在频率为fi的平行光入射时的模型并生成模型文件；
(5)将步骤(4)中得到的模型输入给集群计算机的各个计算节点；集群计算机利用改进了的电磁场时域伪谱方法模拟平面波在该生物组织模型中的传播情况；集群计算机输出记录空间中电磁场分布稳定后的电磁场分布的数据文件；
(6)读取步骤(5)得到的数据文件，用近场外推远场程序计算得到后向远场散射场分布和
(7)根据步骤(6)的结果计算不同散射角的后向散射强度Ii(θ)
Ii(θ)=12∫|Ei→(θ,φ)×Hi→*(θ,φ)|dφ;]]>
(8)判断i＜21是否成立：若是，i＝i+1，并返回步骤(3)；若否，进入步骤(9)；
(7)取不同频率入射时后向散射场的平均
图11为本实施例2的模拟结果：后向散射增强现象明显，且得到了来源于组织浅层后向散射增强角度范围(～2°)。这一模拟对将后向散射增强现象应用到光学检测生物组织中有重要指导意义。
值得指出的是，本发明中集群计算机模拟电磁波传播的方法不仅适用于上述两个实施例中光在生物组织模型中的传播模拟的，还适用于任意电磁波在任意复杂空间媒质中的传播。