地下三维目标电磁散射的矩量法计算
Scattering from 3-D Targets in the Subsurface Using MOM
-
摘要: 该文给出了一种用层状介质中的混合势积分方程(MPIE)和基于RWG基函数的矩量法计算地下三维目标电磁散射的精确快速实施方法。对MPIE的RWG矩量法的开发、计算性能做了研究,尤其是对其中的多个不同形式的Sommerfeld积分的快速全波数值离散复镜像计算方法做了仔细研究。该文的实施方法退化到自由空间后的计算结果与解析解Mie Series吻合的很好,而且地下平板的计算结果也与以往公布结果吻合得很好,证实了该文实施方法的可行、高效、精确。除此以外,该文还计算了其它形状目标在不同大小、不同埋藏深度、以及不同地层媒质下的电磁散射特征。Abstract: In this paper, an accurate and efficient implementation approach of the Method of Moments (MOM) is developed to compute scattering from 3-D targets in the subsurface, which is based on the Mixed Potential Integral Equation (MPIE) in layered media and the RWG basis function. The implementation skill and computing performance of the proposed approach are studied in detail. In particular, how to efficiently calculate various different Sommerfeld integrals (SIs) is numerically investigated with the discrete complex images technique. The computed results of the targets in free space are in agreement with Mie Series results and the computed results of subsurface flat are also in agreement with the published data, which shows the reliability, efficiency and accuracy of the developed approach. Furthermore, the EM scattering of other shaped targets is computed at different size, buried depth, and subsurface medium.
-
1. 引言
空间填充曲线可将高维空间中的对象映射到1维空间,从而将复杂的高维问题降低为简单的1维问题,因此其在图像处理、空间数据管理、数据分区与负载均衡[1-6]、无线传感器网络[7,8]、Cache局部保序[9]等诸多领域得到了广泛的应用。
目前,常见的空间填充曲线包括Z曲线[10]和Hilbert曲线[11]等,大量理论分析与应用表明:Hilbert曲线相对其它曲线具有跳变性小,聚合性高等优点[12]。
尽管Hilbert曲线具有上述优点,但因其映射规则复杂,编码效率低等不足,从而限制了基于Hilbert应用的进一步推广,故许多应用转向编码效率更高的Z曲线[13,14]。Hilbert曲线编码和解码是基于Hilbert应用的前提,提高编解码效率具有重要的意义,故许多研究工作都致力于解决这一问题。
经典的Hilbert编解码算法可被划分为递归和迭代两类。对递归的算法,Butz等人[11]做出了开创性的工作,设计了一种以字节为导向的Hilbert曲线生成算法。在Butz工作的基础上,多种算法被提出,但其复杂度多为O(n2),其中n为编码的阶数,效率较低。Fisher[15]指出多数递归的算法不适用于Hilbert编码或效率较低,并提出了基于双查找表的迭代编码算法,故后续的研究多基于迭代的方式展开。
在基于迭代方式的实现中,基于状态视图的编解码算法因其具有简单、直观、高效、易实现等特性,从而得到了深入地研究。李绍俊等人[16]提出基于状态转移矩阵的、复杂度为O(n)的迭代编码算法,其设计了精简高效的状态转移表并采用位域共同体等优化手段来提高算法效率。考虑到基于状态视图的算法多需要针对编码和解码设计不同的状态视图,为解决这一问题,文献[17]通过编码前的象限映射和解码后的象限逆映射设计了基于一致状态视图的编解码算法,使设计的一个状态视图可同时适用于编解码过程。
除上述算法外,也有一些不依赖状态视图的实现算法。文献[18]通过计算替代查找表,可降低查找表的空间开销,但因其需要额外的计算开销,所以影响算法效率的提升;Moore[19]在Butz经典算法的基础上,设计了基于位操作的高效非递归编码和解码算法,在低维和高维编解码上均具有较好的效果,因此被作为比较的基准算法[20]并在空间范围查询等领域中得到了广泛的应用[21];曹雪峰等人[20]在考虑维度分布差异的基础上,设计了紧致的Hilbert编码,适用于各维度数据分布不均匀的场合,Zhang等人针对3维[22]、Liu等人[23,24]针对高维提出了相应的编解码算法。
考虑到状态视图的优势,本文同样基于状态视图进行编解码算法的研究。通过对现有多种高效的编解码算法对比后发现,现有的高效算法多未考虑输入数据对编码或解码效率的影响,因此将不同输入数据同等对待。例如,李绍俊等人[16]的算法为作者所知最为高效的编码算法,但其编码过程中需逐阶遍历,共需
2n 次查询状态视图的操作,其复杂度为O(n) 。其次,该算法未考虑输入数据的差异,每个输入数据需要等同的编码和解码时间,故其效率有待进一步提升。为此,本文对状态视图深入分析后发现并证明:无需对输入数据中前部为0的阶进行迭代编码或解码。在此基础上,提出复杂度为
O(m) 的免计前0的Hilbert编码算法(Front-Zero-Free Hilbert Encoding, FZF-HE)和免计前0的Hilbert解码算法(Front-Zero-Free Hilbert Decoding, FZF-HD),可将查询状态视图的次数降低为2m 次,其中m 为需迭代的阶数,m≤n 。本文算法考虑了不同输入数据的差别,具有输入坐标点或编码值越小,编解码效率越高的特性,因此能很好地适应数据的偏斜分布。扩展实验结果表明,本文算法在均匀分布下,效率稍高于现有最优算法,而在数据偏斜分布下,效率远优于现有算法。2. Hilbert曲线
Hilbert曲线是空间填充曲线家族的一种,可以填满整个正方形网格,其通过每个网格中心1次且仅1次。Hilbert曲线将空间中每个网格进行编码,1阶Hilbert曲线将正方形网格分割为4个小正方形,小正方形的左下、左上、右下、右上分别对应着第I象限、第II象限、第III象限、第IV象限,将小正方形的第I, II, IV, III象限中心连接起来即可构成1阶Hilbert曲线。通过在每个象限放置一个1阶曲线,并对第I象限左旋90°并垂直翻转,对第III象限右旋90°并垂直翻转可得到2阶曲线,继续迭代可进而得到任意阶曲线。1, 2, 3阶Hilbert曲线的示意图如图1所示。
给定位置点
P(X,Y) ,其中X=(x1x2···xn)2, Y=(y1y2···yn)2 ,确定P 在Hilbert曲线中的序号Z=z1z2···z2n 的过程称作编码;相应地,给定Hilbert编码Z=z1z2···z2n ,确定其对应位置点P(X,Y) 的过程称作解码。其中,X ,Y 为相应的横坐标和纵坐标,Z 为对应的编码,n 称作阶或解析度。xi ,yi ,z2i−1z2i 分别为第i 阶横坐标、纵坐标和编码值,其中1≤i≤n 。3. Hilbert曲线编码和解码
3.1 状态视图的设计
状态视图是本文Hilbert编码和解码算法的基础,故首先对Hilbert曲线的状态进行分析。
对1阶Hilbert曲线而言,有4种基本状态:状态0、状态1、状态2和状态3,分别对应曲线开口朝下、左、上和右4种情形。图2给出了4种状态下1阶Hilbert曲线对应的状态图(括号外为编码,括号内为坐标),不同状态下,每个坐标对应着不同的编码值。
对编码而言,设计状态视图的关键在于构建坐标和Hilbert编码值、坐标和下一阶状态之间的映射关系。根据图2和Hilbert曲线的旋转和翻转规则,易于构建编码状态视图如表1所示。
表 1 编码状态视图状态 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 坐标 00 01 10 11 00 01 10 11 00 01 10 11 00 01 10 11 编码 00 01 11 10 00 11 01 10 10 11 01 00 10 01 11 00 下阶状态 1 0 3 0 0 2 1 1 2 1 2 3 3 3 0 2 解码为编码的逆操作,同理,根据图2和Hilbert曲线的旋转和翻转规则,易于构建解码状态视图如表2所示。
表 2 解码状态视图状态 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 编码 00 01 10 11 00 01 10 11 00 01 10 11 00 01 10 11 坐标 00 01 11 10 00 10 11 01 11 10 00 01 11 01 00 10 下阶状态 1 0 0 3 0 1 1 2 3 2 2 1 2 3 3 0 需要注意的是,先前的工作[18]将第
i 阶横纵坐标值xi 和yi 合并为xiyi 后执行计算及状态视图查询操作,这将导致额外的位合并开销。为避免这一开销且提高状态视图的查询效率,本文设计两个3维数组Key 和Type 来表示状态视图,其中Key 实现某阶坐标值到该阶编码值的映射,而Type 实现某阶状态到下一阶状态的映射。直接通过第i 阶坐标值xi ,yi 和第i 阶状态Ti ,分别按照式(1)和式(2)查询状态视图获得对应的第i 阶编码值z2i−1z2i 和第i+1 阶状态值Ti+1 ,从而避免位合并操作的开销z2i−1z2i=Key[Ti][xi][yi] (1) Ti+1=Type[Ti][xi][yi] (2) 对应地,通过设计两个2维数组
InvKey 和InvType ,其中InvKey 实现某阶编码值到该阶坐标值的映射,而InvType 实现某阶状态到下一阶状态的映射。可根据输入的状态值Ti 和编码值z2i−1z2i ,分别按照式(3)和式(4)查询状态视图获取第i 阶坐标组合xiyi 和i+1 阶状态值Ti+1 xiyi=InvKey[Ti][z2i−1z2i] (3) Ti+1=InvType[Ti][z2i−1z2i] (4) 为便于描述,本文后续编码和解码算法中,规定第1阶的状态总为0,即Hilbert曲线第1阶总为开口向下。
3.2 编码
构建编码状态视图之后,可逐阶迭代查找状态视图,以获取每阶的编码值和下一阶状态,进而可生成最终的编码。
但通过这种方式,需要遍历输入的每一阶的横纵坐标,算法的复杂度为
O(n) ,且这种方式下所有输入数据需要相同的编码时间,并未考虑不同输入数据的差别,因此编码效率较低。故需对其进行改进,以降低复杂度。本文的编码算法依据定理1实现,因此首先给出定理1如下所示:
定理1 若横坐标
x1x2···xk 和纵坐标y1y2···yk 同时为0,则前k 阶编码为0,第k+1 阶的状态等于0(k 为偶数)或1(k 为奇数)。证明
k=1 时,因第1阶状态为0,x1 和y1 为0,则根据表1可知,第1阶的编码为00,第2阶状态为1,故定理成立。k=2 时,因第2阶状态1,x2 和y2 为0,则根据表1可知,第2阶的编码同样为00,前2阶的编码为0000,第3阶的状态为0,故定理成立。k≥3 时,因第3阶的输入(横纵坐标值及状态)与第1阶完全相同,故第3阶的输出(编码值和第4阶状态)与第1阶相同;同理,第4阶的输出与第2阶相同,进入循环状态,定理同样成立。证毕
由定理1可见,对
X ,Y 的前k 阶全为0的情形,其前k 阶编码值为0,且第k+1 阶的状态可由k 的奇偶性直接确定,故无需对前k 阶进行迭代编码,从而可提高算法的效率。如何快速检测
X 和Y 前部同时为0的最大长度k 是算法高效执行的关键。直观来看,可以从X 和Y 的第1阶开始,逐阶探测xi 和yi 是否同时为0。但如此仍需逐阶探测,效率较低。为解决这一问题,引入高效的置位检测算法[25],该算法通过二分查找的思路快速检测给定整数的第1个被置1的位置
p (最右侧位的位置为0),因此可用于快速计算k 值。为避免对X 和Y 各执行1次置位检测,本文仅对X 和Y 的最大值执行置位检测,故1次编码过程,仅需1次置位检测操作。检测到
k 值之后,对k+1 阶及之后的各阶,迭代执行状态视图查询操作,即可获取最终的Hilbert编码。本文算法可免于计算
X ,Y 前部为0的阶,故称之为FZF-HE算法,完整的FZF-HE算法如表3所示。表 3 FZF-HE算法输入:X,横坐标 Y,纵坐标 n,阶 输出:Hilbert编码Z (1) Z←0 (2) p←msb(max(X,Y))//置位检测 (3) t←(n-p-1)%2//计算第n-p阶状态 (4) for i from n-p to n (5) Z=Z<<2|Key[t][xi][yi]//查Key视图 (6) t=Type[t][xi][yi]//查Type视图 FZF-HE算法中msb为置位检测算法,
n− p−1 即为所求k 值,Key 和Type 为3.1节对应的编码状态视图。算法复杂度分析:算法的主要开销在于迭代查询状态视图的操作,共需迭代执行
m=n− k=log2(max(X,Y)) 次,每次执行2次查表操作,因此算法的时间复杂度为O(m) ,因k≥0 ,故本文算法复杂度优于O(n) 。算法空间开销主要为2个状态视图的开销,因Key 表和Type 表分别存储2位二进制表示的编码值和4种状态值,故Key 和Type 均设计为char类型,最终需要的空间为2×4×4×sizeof(char)=32 Byte,其空间开销可基本忽略不计。从上述分析来看,
max(X,Y) 越小,FZF-HE的时间复杂度越低,因此该算法适合数据偏斜分布,尤其是向原点(正方形区域的左下角)偏斜的情形。在实际应用中,数据偏斜分布是普遍存在的,故本文算法可有效提升偏斜分布数据的编码效率①。3.3 解码
Hilbert解码是将编码映射到坐标的过程,是编码的逆过程。类似地,本文实现的FZF-HD算法依据定理2实现。
定理2 若
z1z2···z2k 为0,则前k 阶对应的横坐标x1x2···xk 和纵坐标y1y2···yk 为0,第k+1 阶的状态等于0(k 为偶数)或1(k 为奇数)。证明
k=1 时,因第1阶状态为0,z1z2 为00,则根据表2的状态视图可知,第1阶的坐标x1 和y1 均为0,第2阶状态为1,故定理成立。k=2 时,因第2阶状态1,z3z4 为00,则根据状态表2可知,第2阶的坐标x2 和y2 同样为0,第3阶的状态为0,故定理成立。k≥3 时,因第3阶的输入(编码值及状态)与第1阶完全相同,故第3阶的输出(横纵坐标和第4阶状态)与第1阶相同;同理,第4阶的输出和与第2阶相同,进入循环状态,定理同样成立。证毕
基于定理2, FZF-HD无需对
Z 前部编码值为0的阶进行解码。同编码算法,本文同样采用置位检测算法[25]来检测编码Z 中第1个被置位为1的位置p ,并根据定理2直接确定编码值前部全为0的前k 阶(k=n−p/2−1 )的坐标值和第k+1 阶状态。对k+1 阶及之后的各阶,迭代查询状态视图InvKey 和InvType ,逐阶获取该阶对应的坐标值和下一阶状态。最终实现的FZF-HD算法如表4所示。表 4 FZF-HD算法输入:Z,Hilbert编码 n,阶 输出:X,横坐标 Y,纵坐标 (1) X, Y←0 (2) p←msb(Z)//置位检测 (3) t←(n-p/2-1)//计算第n-p/2阶状态 (4) for i from n-p/2 to n (5) coord=InvKey[t][z2i−1z2i] //查InvKey视图 (6) Y=Y<<1|coord&0x1 X=X<<1|coord>>1&0x1 (7) t=InvType[t][z2i−1z2i]//查InvType视图 算法复杂度分析:算法的主要开销在于迭代查询状态视图的操作,算法共需迭代
m=n− k=log2(Z) 次,每次执行2次查表操作,因此算法的时间复杂度为O(m) 。算法的空间开销主要为2个状态视图的开销,因InvKey 和InvType 分别存储2位二进制表示的坐标值和4种状态值,故均设计为char类型,最终需要的空间同样为2×4×4×sizeof(char)=32 Byte。4. 实验结果比较及分析
4.1 实验环境及数据集
实验主要硬件平台:Intel(R) Core(TM) i5-3230m CPU @ 2.60 GHz双核,内存4 GB。软件环境:windows 10 64位,Microsoft Visual Studio C++ 2015且禁用了优化(/ Od)。
本文的数据集考虑均匀分布和偏斜分布两种情况。对均匀分布,分别为4, 8, 16, 24, 32阶各随机生成106个坐标点(用于编码)及编码值(用于解码)。对偏斜分布,为考察数据偏斜程度对算法的影响,引入偏斜阶数
α 和偏斜率β 两个参数。其中α 表示前α 阶取值全为0,β 则表示生成的数据中前α 阶取值全为0的部分占全部数据的百分比。例如,α=10 和β=40% ,表示随机生成40%的数据的前10阶全为0。本文为每个特定的α 和β ,各生成1百万条32阶的随机坐标数据(用于编码)或编码数据(用于解码)作为实验数据。为考察本文算法的有效性,将本文算法与本领域代表性的算法Li[16], Burkardt[18], Moore[19], uniState[17] 4种算法进行对比。
4.2 编码
4.2.1 均匀分布
为了考察均匀分布下阶数对不同算法编码效率的影响,给定上述生成的均匀分布的数据,比较不同算法的时间开销如图3所示。
由图3可见,所有算法的编码时间均随着阶数的增加而增加。相比较而言,FZF-HE算法优于其余4种算法,且其编码时间随阶数增加的幅度最小。以32阶为例,FZF-HE编码仅需0.29 s,Li的算法需要0.52 s,其它3种算法均需1 s以上。这些结果表明FZF-HE算法通过避免对坐标前部全为0的阶进行处理能有效提升编码效率。
4.2.2 偏斜分布
为考察偏斜分布下不同算法的效率,分别从以下两个方面进行实验对比:(1)固定
α=16 ,设置β 分别为25%(对应均匀分布),50%, 75%, 100%;(2)固定β=50% ,设置α 分别为8, 12, 16, 20, 24,实验结果分别如图4和图5所示。由图4可见,在
α=16 时,FZF-HE的编码时间大幅低于其他算法,且随着β 的增大而明显降低,而其它算法则不随β 的变化而明显变化。例如,在β 为100%时,FZF-HE的编码时间为0.17 s,仅为另外4种算法中最优的Li的算法时间的47.6%。分析可知对FZF-HE而言,随着β 的增大,前α 阶为0的数据量逐步增大,需迭代编码的部分降低,故算法效率逐步提高。由图5可见,FZF-HE算法的编码时间最低,且随着
α 增大而减小,而其他算法则不随α 的变化明显变化。这是因为随着α 增大,数据中无需编码的阶增加,所以算法效率逐步提高。4.3 解码
4.3.1 均匀分布
为了考察均匀分布下阶数对不同算法解码效率的影响,给定上述生成的均匀分布的数据,比较不同算法的时间开销如图6所示。
由图6可见,解码结果与编码结果类似,总体来看,FZF-HD算法仍为最优的算法,同样表明FZF-HD算法通过避免对编码值前部全为0的阶进行解码,从而能有效提升解码效率。
4.3.2 偏斜分布
为考察偏斜分布对不同算法的影响,分别从以下两个方面进行实验对比:(1)固定
α=16 ,设置β 分别为25%(对应均匀分布),50%, 75%, 100%; (2)固定β=50% ,设置α 分别为8, 12, 16, 20, 24,实验结果分别如图7和图8示。由图7可见,解码的实验结果与编码类似,FZF-HD的解码时间最低,且随着
β 的增大而降低,其它算法则不随β 的变化而明显变化。当β 为100%时,FZF-HD解码106个32阶数据仅需0.14 s,效率同样远优于其它算法。对FZF-HD而言,随着β 的增大,前α 阶为0的数据量逐步增大,故需迭代解码的部分降低,算法效率逐步提高。由图8可见,FZF-HD算法的解码时间随着
α 增大而减小,其它算法则不随α 变化而明显变化;因为随着α 的增大,数据中无需迭代解码的阶增加,故算法效率逐步提高。5. 结束语
本文在设计高效状态视图并结合快速置位检测技术的基础上,提出高效的编码算法FZF-HE和解码算法FZF-HD,可避免对输入数据前部全为0的阶进行编码或解码,从而可降低查询状态视图的次数,提高编解码的效率。扩展实验表明,本文算法优于现有算法且更适合数据偏斜分布的场合。下一步工作拟将FZF-HE和FZF-HD算法应用于先前采用Z曲线解决的问题或领域[13,14],以扩展Hilbert编解码算法的应用范围;此外,通过并行化等手段研究进一步提高Hilbert编解码算法效率的方法。
-
Michalski K A, Mosig J R. Multilayered media Greens functions in integral equation formulations[J].IEEE Trans. on Antennas and Propagation.1997, 45(3):508-519[2]Harrington R F. Field Computation by Moment Methods. New[3]York:Macmillan,1968, Melboume, FL: Krieger, reprint 1982.[4]Chao J C, Liu Y J, Rizzo F J, Martin P A, Udpa L. Regularized integral equations for curvilinear boundary elements for[5]electromagnetic wave scattering in three dimensions. IEEE Trans[J].on Antennas and Propagation.1995, 43:(12):1416-1422[6]Hua Y, Sarkar T K. Generalized pencil-of-function method for extracting poles of an EM system from its transient response[J].IEEE Trans. Antennas Propagat.1989, 37(2):229-234[7]Aksun M I. A robust approach for the derivation of closed-form Greens functions[J].IEEE Trans. Microwave Theory Tech.1996, 44(5):651-658[8]盛新庆. 计算电磁学要论. 北京:科学出版社,2004: 21-27.[9]Rao S M, Wilton D R, Glisson A W. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Trans. Antennas Propagation., 1982, AP-30(3): 409-418. [8] Fang D G, Yang J J, Delisle G Y. Discrete image theory for horizontal electric dipoles in a multilayered medium. IEE Proc.-H, 1988, 135(5): 297-303.[10]Pasi Yla-Oijala, Matti Taskinen. Efficient formulation of closed-form Greens functions for general electric and magnetic sources in multilayered media[J].IEEE Trans. on Antennas and Propagation.2003, 51(8):2106-2115[11]Cui T J, Wiesbeck W, Herschlein A. Electromagnetic scattering by multiple three-dimensional scatterers buried under multilayered mediaPart II: Numerical implementations and results[J].IEEE Trans. Geosci. Remote Sensing.1998, 36(2):535-546 期刊类型引用(8)
1. 黄炎,李姗姗,吕明昊,范雕,谭勖立,冯进凯. 利用CPU和GPU混合并行方法快速构建海洋扰动重力梯度基准图. 武汉大学学报(信息科学版). 2025(03): 515-527 . 百度学术
2. 贾连印,范瑶,丁家满,李晓武,游进国. 高效前缀约简的三维Hilbert空间填充曲线编解码算法. 电子与信息学报. 2024(02): 633-642 . 本站查看
3. 王维晨,贾连印,王炳月,梁彬彬,卫守林. 一种自适应的Hilbert编码算法及其并行化. 计算机应用与软件. 2024(04): 236-241 . 百度学术
4. 杨飞,华一新,杨振凯,李响,赵鑫科,张晓楠. 一种基于迭代法的非均匀Hilbert空间填充曲线快速生成算法. 武汉大学学报(信息科学版). 2024(06): 1040-1050 . 百度学术
5. 龚艳茹,马立坤. 多仓库物流周期联合配送粒子群择优算法仿真. 计算机仿真. 2023(09): 92-96 . 百度学术
6. 李娜,田文文. 基于Delaunay三角形网格的多视点视频编解码算法. 科技通报. 2023(11): 19-23 . 百度学术
7. 王炳月,贾连印,范瑶,孙劭文. 一种优化的灰度图像压缩算法. 电视技术. 2022(05): 17-23 . 百度学术
8. 贾连印,孔明,王维晨,李孟娟,游进国,丁家满. 数据偏斜分布下的二维Hilbert编解码算法. 清华大学学报(自然科学版). 2022(09): 1426-1434 . 百度学术
其他类型引用(1)
-
计量
- 文章访问数: 2657
- HTML全文浏览量: 97
- PDF下载量: 719
- 被引次数: 9