CN116449324A - 一种应用于电磁散射的h-自适应矩量方法 - Google Patents
一种应用于电磁散射的h-自适应矩量方法 Download PDFInfo
- Publication number
- CN116449324A CN116449324A CN202310357007.8A CN202310357007A CN116449324A CN 116449324 A CN116449324 A CN 116449324A CN 202310357007 A CN202310357007 A CN 202310357007A CN 116449324 A CN116449324 A CN 116449324A
- Authority
- CN
- China
- Prior art keywords
- grid
- charge density
- grid cell
- super
- error
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 84
- 230000006870 function Effects 0.000 claims abstract description 53
- 238000011084 recovery Methods 0.000 claims abstract description 36
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 19
- 230000005284 excitation Effects 0.000 claims abstract description 12
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 2
- 238000003780 insertion Methods 0.000 abstract description 4
- 230000037431 insertion Effects 0.000 abstract description 4
- 238000004088 simulation Methods 0.000 description 17
- 230000003044 adaptive effect Effects 0.000 description 12
- 230000008569 process Effects 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 241001274197 Scatophagus argus Species 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 239000004020 conductor Substances 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 239000000243 solution Substances 0.000 description 3
- 230000005672 electromagnetic field Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012966 insertion method Methods 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000012088 reference solution Substances 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- Software Systems (AREA)
- Operations Research (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种应用于电磁散射的h‑自适应矩量方法,包括:创建用于计算待测目标散射特性的电磁模型,进行平面波激励下的初始网格划分;利用基于RWG基函数的矩量法求解当前各网格单元的电流系数,继而计算对应网格单元中超收敛点处的面电荷密度;基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法计算对应网格单元的局部误差;基于各网格单元的局部误差计算全局误差;判断是否满足迭代停止条件;若否选择部分网格单元基于B‑W插点算法进行网格细化,得到重新划分的网格并返回电流系数求解步骤;若是停止迭代,输出当前求解出的电流系数。本发明能够保证计算精度、鲁棒性、减少误差估计、减少网格细化数量并减少计算资源消耗。
Description
技术领域
本发明属于电磁场领域,具体涉及一种应用于电磁散射的h-自适应矩量方法。
背景技术
矩量法(Method of Moments,MoM)作为一种经典的数值分析方法,凭借计算精度高、普适性强等优点,已广泛应用于电磁计算中,如隐身与反隐身技术需要对目标的电磁散射特性进行精确分析等。而针对电磁领域的高精度计算,精准的网格剖分是一个必要的前提,原因在于,如果剖分尺寸过大,无法准确描述目标的细小结构或场源附近电磁流的复杂变化,会产生数值误差;如果剖分尺寸过小,会使网格数量呈非线性增长,导致所需的计算资源急剧增加。针对网格剖分,剖分正确性的判断和重新剖分都有可能造成大量的时间浪费。
目前,一种通过控制局部网格尺寸以减小数值解误差的方案,即h-自适应细化,可有效解决上述问题。该方案包括误差估计器和控制器;误差估计器用于计算网格单元误差,而控制器决定对哪些网格进行细化并设置迭代终止条件。
2018年,S K Kim提出了一种基于电流密度不连续性误差估计的h-自适应矩量法,该方法利用网格单元之间的电流密度不连续性计算网格单元的误差值,这种误差估计方法消耗时间极少,但鲁棒性较低。
与之相比,其余常用的误差估计方法通常为残差误差估计,即通过改变选配过程中的检验函数或边界条件进行误差估计。这种方法虽然鲁棒性高,但需要进行与矩阵填充相当的计算并广泛使用数值积分,导致误差估计时间很长。而其余计算时间较短的误差估计方法如不连续性误差估计在很多情况下又会失效。可见,在保证计算精度、鲁棒性要求的同时,减少误差估计时间是目前的一个技术难点。并且,同时尽可能地减少网格细化数量以减少计算资源消耗也是本领域内一个亟待解决的技术问题。
发明内容
为了解决现有技术中存在的上述问题,本发明提供了一种应用于电磁散射的h-自适应矩量方法。本发明要解决的技术问题通过以下技术方案实现:
创建用于计算待测目标散射特性的电磁模型,进行平面波激励下所述电磁模型的初始网格划分,得到待测目标表面划分出的三角形网格;
利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度;
基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差;
基于各网格单元的局部误差计算全局误差;
基于所述全局误差判断是否满足迭代停止条件;
若否,选择部分网格单元基于B-W插点算法进行网格细化,得到重新划分的网格;并返回所述利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度的步骤;
若是,停止迭代,输出当前求解出的各网格单元的电流系数。
在本发明的一个实施例中,所述进行平面波激励下所述电磁模型的初始网格划分,包括:
将所述电磁模型以λ4的尺寸进行初始网格划分;其中,λ表示所述平面波的波长。
在本发明的一个实施例中,所述利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度,包括:
利用所述基于RWG基函数的矩量法和当前划分的网格对待测目标表面的等效电磁流进行离散,将利用网格离散后产生的基函数组成线性方程组,求解所述线性方程组得到各网格单元的电流系数;
利用各网格单元的电流系数和面电荷密度计算公式,计算出对应网格单元中超收敛点处的面电荷密度。
在本发明的一个实施例中,所述基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差,包括:
针对属于各网格单元的每一个顶点,利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值;
针对每个网格单元,将该网格单元的三个顶点分别得到的恢复电荷密度值利用最小二乘法进行线性拟合,构造出该网格单元的线性恢复电荷密度函数;
为保证RWG基函数关于三角形对中正负三角形上电荷量的特性,通过将该网格单元的线性恢复电荷密度函数与常数相加的形式,得到该网格单元的修正恢复电荷密度函数;
基于该网格单元的修正恢复电荷密度函数和原超收敛点处的面电荷密度,计算该网格单元的局部误差。
在本发明的一个实施例中,所述利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值,包括:
将包含该顶点的所有网格单元中超收敛点处的面电荷密度进行平均,将得到的平均值确定为该顶点的恢复电荷密度值。
在本发明的一个实施例中,计算该网格单元的局部误差所采用的公式,包括:
其中,表示第i个网格单元的局部误差;ρi',rec表示第i个网格单元的修正恢复电荷密度函数;ρi表示第i个网格单元中超收敛点处的面电荷密度;Si表示第i个网格单元。
在本发明的一个实施例中,所述基于各网格单元的局部误差计算全局误差,包括:
针对每个网格单元,将其三角形面积作为自身权重对该网格单元的局部误差进行加权,得到该网格单元的加权局部误差;并求取所有网格单元的加权局部误差之和,得到加权局部误差和;
将所述加权局部误差和利用所有网格单元的三角形面积和进行归一化;
将归一化结果求取平方根,得到全局误差。
在本发明的一个实施例中,所述基于所述全局误差判断是否满足迭代停止条件,包括:
当所述全局误差降低为预设的初始全局误差的1/Q以下,或者迭代次数达到预设的最大迭代次数,则确定满足迭代停止条件;否则,确定不满足迭代停止条件;其中Q为大于1的自然数。
在本发明的一个实施例中,所述选择部分网格单元基于B-W插点算法进行网格细化,包括:
计算所有网格单元中最大的加权局部误差与预设百分比的乘积得到加权局部误差阈值;
确定加权局部误差大于所述加权局部误差阈值的网格单元作为选定网格单元;
利用基于B-W插点算法对各选定网格单元进行网格细化。
在本发明的一个实施例中,利用基于B-W插点算法对任一选定网格单元进行网格细化,包括:
在当前Delaunay三角化Mi中搜索到覆盖待插入点pi的基单元;其中,当前Delaunay三角化Mi表示当前划分出的三角形网格;
利用网格单元的相邻关系在基单元的相邻网格单元处标记所有不符合外接圆准则的网格单元,并将这类网格单元的集合记为Bi;
在当前Delaunay三角化Mi中删除Bi所覆盖的网格单元形成空腔Ci,连接pi和Ci边界上的顶点,在空腔内形成新的三角化Bi+1,得到新的Delaunay三角化Mi+1。
本发明的有益效果:
本发明实施例所提供的应用于电磁散射的h-自适应矩量方法中,首先创建用于计算待测目标散射特性的电磁模型,进行平面波激励下所述电磁模型的初始网格划分,得到待测目标表面划分出的三角形网格;其次,利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度;接下来,基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差;之后,基于各网格单元的局部误差计算全局误差;然后基于所述全局误差判断是否满足迭代停止条件;若否,选择部分网格单元基于B-W插点算法进行网格细化,得到重新划分的网格,并返回利用基于RWG基函数的矩量法求解电流系数的步骤;若是,停止迭代,输出当前求解出的各网格单元的电流系数。
本发明实施例以控制网格数量的同时得到精准的数值结果为目标,从实现准确而高效的误差估计器和控制器入手,提出一种SPR误差估计器,该误差估计器将超收敛点恢复(SPR)法与RWG基函数特性相结合,能够高效且准确地求得网格单元的误差,且鲁棒性强、计算时间短;在此基础上实现的基于该误差估计器的h-自适应矩量法计算精度高,能有效控制计算所需资源,可以做到无人工干预的情况下对网格自适应细化,并且能够在设置好迭代停止条件后自动停止,可以有效地控制网格数量并得到准确的仿真结果。
附图说明
图1为本发明实施例所提供的一种应用于电磁散射的h-自适应矩量方法的流程示意图;
图2为本发明实施例中三角形网格及网格单元的超收敛点示意图;
图3(a)~图3(d)为本发明实施例中利用基于B-W插点算法对选定网格单元进行网格细化的过程示意图;
图4为本发明实施例仿真实验中的圆锥体的电磁模型和平面波激励示意图;
图5为本发明实施例仿真实验中圆锥体网格数量和全局误差随自适应迭代变化图;
图6为本发明实施例仿真实验中圆锥体初始网格分布图;
图7为本发明实施例仿真实验中圆锥体自适应细化后网格分布图;
图8为本发明实施例仿真实验中圆锥体自适应细化与初始网格双站RCS对比图;
图9为本发明实施例仿真实验中圆锥体常规剖分网格分布图;
图10为本发明实施例仿真实验中圆锥体自适应细化与常规剖分双站RCS对比图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
为了便于理解本发明实施例的方案,在此先对相关技术和本发明实施例的发明构思进行简介。
针对目标隐身技术、反隐身探测技术等应用场景,获取目标的雷达散射截面具有重要的工程意义。但是,通过实验获取数据需要同时考虑场地、资金、技术等问题,实施起来较为困难;而电磁仿真作为一种高效、方便、实用的手段,在目标雷达散射截面获取中越来越受到青睐。本发明实施例首先对现有矩量法等算法进行研究。其中,基于RWG基函数的矩量法是通过三角形网格对目标表面的等效电磁流进行离散,不精准的网格剖分会在求解时带来误差,如一些曲面在剖分时会被平面网格取代,进而造成模型形变;一些细小结构处由于电磁流变化剧烈,对网格的要求也同样很高。本发明实施例为了量化这个误差并依据这个误差对网格尺寸进行控制,并在控制网格数量的同时得到精准的数值结果,针对目标的电磁散射特性计算,从如何实现准确而高效的误差估计器和控制器入手,提出了一种基于超收敛点恢复法(Superconvergent Patch Recovery,SPR)误差估计的h-自适应矩量法。具体的,本发明实施例是将超收敛点恢复法引入矩量法的误差估计中,并结合RWG基函数的特性对误差估计方法进行优化,配合快速B-W插点算法对网格进行细化。
如图1所示,本发明实施例所提供的一种应用于电磁散射的h-自适应矩量方法,可以包括如下步骤:
S1,创建用于计算待测目标散射特性的电磁模型,进行平面波激励下所述电磁模型的初始网格划分,得到待测目标表面划分出的三角形网格;
本发明实施例中的待测目标可以是材料为理想电导体(PEC)的任意几何结构的电磁散射体(electromagnetic scattering),如飞机、轮船、导弹等。
为了计算待测目标的散射特性,需要根据所述待测目标的结构、尺寸以及边界面的连续性,构建与设计要求相对应的几何模型;并将该几何模型的外表面设置为理想电导体(Perfect Electric Conductor,PEC)材料,创建出待测目标的电磁模型。之后,为所述电磁模型设置平面波激励,再进行初始网格划分。
可选的一种实施方式中,所述进行平面波激励下所述电磁模型的初始网格划分,包括:
将所述电磁模型以λ/4的尺寸进行初始网格划分;
其中,λ表示所述平面波的波长。
经过初始网格划分,待测目标表面生成三角形网格,每个三角形对应一个网格单元。
S2,利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度;
可选的一种实施方式中,S2可以包括:
S21,利用所述基于RWG基函数的矩量法和当前划分的网格对待测目标表面的等效电磁流进行离散,将利用网格离散后产生的基函数组成线性方程组,求解所述线性方程组得到各网格单元的电流系数;
该步骤的具体过程可以参见相关技术理解:
[In]=[zmn]-1·[gm]
其中,In为电流系数,zmn为阻抗矩阵,gm为激励矩阵。
最终,通过RWG基函数表征金属表面的电流密度为:
其中,In为电流系数,Ns为基函数个数,为基函数。
可以理解的是,通过求解得到的电流系数可以获取到表面电流密度,依据现有公式可以通过表面电流密度计算得到目标的雷达散射截面、天线的增益、方向系数、微波器件的S参数等。
S22,利用各网格单元的电流系数和面电荷密度计算公式,计算出对应网格单元中超收敛点处的面电荷密度。
超收敛点恢复法指出三角形的质心上具有超收敛性,因为质心距离存在不连续性的边最远,因此超收敛点即质心上的电磁变量计算精度最高,也就是说,三角形质心处的电磁变量可用于误差分析,因此,本发明实施例将超收敛点恢复法引入矩量法的误差估计中,针对每个网格单元,利用面电荷密度计算公式和该网格单元的电流系数,计算出该网格单元的质心处的面电荷密度。
其中,面电荷密度计算公式为:
其中,为表面电流密度,ω为角频率。
S3,基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差;
该步骤由本发明实施例的SPR误差估计器实现。可选的一种实施方式中,S3可以包括以下步骤:
S31,针对属于各网格单元的每一个顶点,利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值;
其中,利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值,包括:
将包含该顶点的所有网格单元中超收敛点处的面电荷密度进行平均,将得到的平均值确定为该顶点的恢复电荷密度值。
具体请参见图2所示的三角形网格及网格单元的超收敛点示意图;其中,示例性地给出了五个三角形,每个三角形的质心分别以ρ1,ρ2,ρ3,ρ4,ρ5表示,五个三角形的共有顶点表示为K。该步骤对网格单元上某一顶点分析,如图2所示,通过对包含该顶点的所有网格单元的超收敛点处的面电荷密度进行平均,所得的平均值即为该顶点的恢复电荷密度值:
其中,ρk表示该顶点的恢复电荷密度值;ρi表示包含该顶点的所有网格单元的超收敛点处的面电荷密度;h表示包含该顶点的所有网格单元的数量。
可以理解的是,通过SPR方法,每个三角形对应的网格单元的顶点都得到了一个恢复电荷密度值。
S32,针对每个网格单元,将该网格单元的三个顶点分别得到的恢复电荷密度值利用最小二乘法进行线性拟合,构造出该网格单元的线性恢复电荷密度函数;
该步骤中的最小二乘法拟合求解过程可以表示为:
1)设ρrec=Pa;
其中,ρrec表示该网格单元的线性恢复电荷密度函数;P表示坐标向量;a表示系数向量;
P=[1,x,y,z];a=[a1,a2,a3,a4]T
2)由最小二乘法拟合公式
其中,(xi,yi,zi)表示对应点的三个坐标值。
可知,当F(a)的最小值时:
解得:
a=A-1b
其中:
然后将a带入ρrec=pa中即可求得线性恢复电荷密度函数ρrec。
S33,为保证RWG基函数关于三角形对中正负三角形上电荷量的特性,通过将该网格单元的线性恢复电荷密度函数与常数相加的形式,得到该网格单元的修正恢复电荷密度函数;
根据RWG基函数特性:对于任意一组三角形对,其正负三角形上的电荷量为常数,且总和为零。为保证此特性,本发明实施例针对每个网格单元,需要将构造出的线性恢复电荷密度函数加上一个常数,使得该线性恢复电荷密度函数在三角形中心点的恢复电荷密度与该点原始的面电荷密度相同。基于此目标,可以确定一个常数,将该网格单元的线性恢复电荷密度函数加上这个常数的和确定为该网格单元的修正恢复电荷密度函数。
S34,基于该网格单元的修正恢复电荷密度函数和原超收敛点处的面电荷密度,计算该网格单元的局部误差。
具体的,计算该网格单元的局部误差所采用的公式,包括:
其中,表示第i个网格单元的局部误差;ρi',rec表示第i个网格单元的修正恢复电荷密度函数;ρi表示第i个网格单元中超收敛点处的面电荷密度;Si表示第i个网格单元。
针对每个网格单元,按照上述步骤均可以计算出一个局部误差。
本发明实施例对误差估计器的处理过程进行了改进,具体利用基于RWG基函数矩量法是通过三角形网格对目标进行离散的特征,将超收敛点恢复法引入矩量法的误差计算中。由于恢复方法依据邻接网格的面电荷密度得到的是更高阶、更准确的数值结果,以恢复解作为参考解进行误差估计相比直接对比网格单元面电荷密度的不连续性更精准,因此鲁棒性强;并且,计算每个网格单元的误差只需要对邻接网格的面电荷密度进行简单的数乘运算且只涉及到一次数值积分,因此计算时间短。因此,本发明实施例的该误差估计器将超收敛点恢复(SPR)法与RWG基函数特性相结合,能够高效且准确地求得网格单元的误差,且鲁棒性强、计算时间短,能够在保证准确性的同时兼顾高效性,减少网格细化数量并节省误差估计消耗的时间。
S4,基于各网格单元的局部误差计算全局误差;
该步骤由控制器实现。可选的一种实施方式中,所述基于各网格单元的局部误差计算全局误差,包括:
S41,针对每个网格单元,将其三角形面积作为自身权重对该网格单元的局部误差进行加权,得到该网格单元的加权局部误差;并求取所有网格单元的加权局部误差之和,得到加权局部误差和;
S42,将所述加权局部误差和利用所有网格单元的三角形面积和进行归一化;
S43,将归一化结果求取平方根,得到全局误差。
上述步骤以公式表示为:
其中,eglobal表示全局误差;表示第i个网格单元的局部误差;Ai表示第i个网格单元的三角形面积;N表示网格单元的数量。
S5,基于所述全局误差判断是否满足迭代停止条件;
其中,该步骤可以包括:
当所述全局误差降低为预设的初始全局误差的1/Q以下,或者迭代次数达到预设的最大迭代次数,则确定满足迭代停止条件;否则,确定不满足迭代停止条件;即上述两个条件均不满足时,确定不满足迭代停止条件。
其中Q为大于1的自然数。具体数值可以根据需要设置,比如,可选的一种实施方式中,Q为4等。
同样的,预设的最大迭代次数也可以根据需要设置,比如可以为10次等等。
若否,执行S6,选择部分网格单元基于B-W插点算法进行网格细化,得到重新划分的网格;并返回所述利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度的步骤;即返回S2步骤。
具体的,针对S6步骤对应的选择部分网格单元基于B-W插点算法进行网格细化,得到重新划分的网格,是指依据各网格单元的局部误差,对大误差区域的网格进行细化,由控制器实现。
具体的,S6可以包括:
S61,计算所有网格单元中最大的加权局部误差与预设百分比的乘积得到加权局部误差阈值;
该步骤中得到的加权局部误差阈值可以表示为:
其中,u%表示预设百分比,具体数值可以根据需要合理设置,比如可以为70%等,u为大于50的自然数;表示所有网格单元中最大的加权局部误差。这里采用的加权局部误差是防止细化过程中网格面积差异过大,导致网格整体质量变差,同时还能保证未被包含在本次细化的误差较大的网格单元满足下一次迭代的大误差区域标准。
S62,确定加权局部误差大于所述加权局部误差阈值的网格单元作为选定网格单元;
本发明实施例中,大误差区域,即选定网格单元,定义为加权局部误差大于所述加权局部误差阈值的网格单元。以预设百分比为例表示为:
通过上述处理可以确定至少一个选定网格单元。
S63,利用基于B-W插点算法对各选定网格单元进行网格细化。
B-W插点算法基于Delaunay三角化的外接圆准则,属于增量插点法。假设Tn是点集Pn={p1,p2,…,pn}的Delaunay三角化,对Tn中任意一个单纯形体t,求外接球半径和球心分别为rt和qt。在Tn中插入新点pn+1,定义所有外接球包含pn+1的单纯形体组成集合为B,
B={t|t∈Tn,d(pn+1,qt)<rt}
若保证pn+1在Pn的凸包内,则B非空,删除B内所有单纯形体,在Tn内形成空腔C,理论上可保证空腔C为星形(star-shaped),且其上所有点对pn+1可视。连接pn+1和空腔C上的点,可得到点集Pn+1={p1,p2,…,pn,pn+1}的Delaunay三角化Tn+1。
集合上述理论,请参见图3(a)~图3(d)理解单次B-W增量插点法的处理过程,具体的,利用基于B-W插点算法对任一选定网格单元进行网格细化,包括以下步骤:
①在当前Delaunay三角化Mi中搜索到覆盖待插入点pi的基单元;
其中,当前Delaunay三角化Mi表示当前划分出的三角形网格;请参见图3(a),其中所有的三角形组成当前Delaunay三角化Mi。待插入点pi为图3(b)中绿色小方框表示的点。基单元是图3(b)中包含pi的那一个绿色边的三角形。
②利用网格单元的相邻关系在基单元的相邻网格单元处标记所有不符合外接圆准则的网格单元,并将这类网格单元的集合记为Bi;
其中,Bi为图3(b)中黑色边内的三个三角形。
③在当前Delaunay三角化Mi中删除Bi所覆盖的网格单元形成空腔Ci,连接pi和Ci边界上的顶点,在空腔内形成新的三角化Bi+1,得到新的Delaunay三角化Mi+1。
其中,Ci为图3(c)中黑边组成的空腔;Bi+1为图3(d)中黑色边内的五个三角形;Mi+1由图3(d)中所有三角形组成。
通过下式可以形象地描述整个过程:
Mi+1=Mi-Bi+Bi+1
具体过程可以参见相关技术理解,在此不做详细说明。
本发明实施例基于B-W插点算法进行网格细化,能够满足自适应过程中因为迭代不断改变需要细化的网格单元以满足数值模拟精度的要求,能够实现无人工干预对网格进行自适应细化,得到高精度的仿真结果。
若是,执行S7,停止迭代,输出当前求解出的各网格单元的电流系数。
可以理解的是,利用各网格单元的电流系数可以进一步实现相关应用,比如,得到目标的雷达散射截面、天线的增益、方向系数、微波器件的S参数等,在此不作详细说明。
本发明实施例所提供的方案中,以控制网格数量的同时得到精准的数值结果为目标,从实现准确而高效的误差估计器和控制器入手,提出一种SPR误差估计器,该误差估计器将超收敛点恢复(SPR)方法与RWG基函数特性相结合,能够高效且准确地求得网格单元的误差,且鲁棒性强、计算时间短,在此基础上实现的基于该误差估计器的h-自适应矩量法计算精度高,能有效控制计算所需资源,可以做到无人工干预的情况下对网格自适应细化,并且能够在设置好迭代停止条件后自动停止,可以有效地控制网格数量并得到准确的仿真结果。
为了验证本发明实施例方法的有效性,以下以仿真实验进行说明。
仿真用的圆锥体半径为100mm,高为500mm,材料为PEC,频率为15GHz垂直极化平面波向圆锥体尖端方向入射,相应的电磁模型和平面波激励如图4所示。
大误差区域,即选定网格单元的定义为:
采用依据SPR误差估计器结果的控制器对圆锥体进行自适应细化。自适应细化的初始网格设定为最大尺寸1/4λ的常规剖分,初始网格数量为6480。全局误差在第六次迭代后降低到初始全局误差的1/4以下,满足迭代停止条件,自适应细化结束。请参见图5所示的圆锥体网格数量和全局误差随自适应迭代变化图。每次迭代对应的网格数量和全局误差请见表1所示。
表1圆锥体自适应细化每次迭代后的网格数量与全局误差
一方面,为了体现误差估计器的准确性,从全局误差的变化角度进行分析。图5是网格数量和全局误差随自适应迭代的变化,表1是自适应细化每次迭代后的网格数量和全局误差。图6和图7是初始和自适应细化后的网格分布图,并给出了局部特写,其中,图6的网格数量为6480;图7的网格数量为15120。由图5可得:在SPR误差估计器地标记下,每次自适应迭代在增加较少网格数量的同时,全局误差都出现了明显地下降。这反映了SPR误差估计器在本例中的准确性。
另一方面,从该模型不同网格剖分下的仿真结果角度进行分析。以该模型在商业软件FEKO中剖分88544网格数量计算的双站RCS(雷达散射截面)作为参考结果。由图8可得,模型经过自适应细化后的双站RCS相比初始网格更接近参考结果,尤其是在θ=-180°~0°区间有明显提升。其中,-180°≤θscat≤180°,f=15GHz。/>表示观测方向与+Z轴的夹角;θscat表示观测方向与+X轴的夹角;在这里/>-180°≤θscat≤180°表示的就是xoz平面。
最后,为了体现自适应细化相比常规剖分的优势。如图9,对圆锥体进行0.1λ的常规剖分,网格数量为32836。图10是自适应细化和常规剖分双站RCS的对比,由图10可得:三条曲线结果重合,自适应细化在网格数量远少于常规剖分的情况下,仿真结果与常规剖分基本相同。由表2,自适应细化相比常规剖分减少了约55%的网格数量,节省了约80%的内存消耗,在满足高计算精度的同时,减少计算资源的效果十分明显。
表2圆锥体自适应细化和常规剖分计算资源对比
剖分类型 | 自适应细化 | 常规剖分 | 减少比例(%) |
网格数量 | 15120 | 32836 | 54.0 |
内存占用(MB) | 7848 | 37017 | 78.8 |
综上,本发明实施例提出的SPR误差估计器在计算时间上相比普遍采用的残差误差估计器存在明显优势,略长于不连续误差估计器;在准确性方面与残差误差估计器性能相当,同时明显高于不连续误差估计器。本发明实施例能够在保证计算精度的同时,实现尽可能短的误差估计时间和尽可能少的网格细化数量,能够大大降低仿真时间,减少计算资源消耗。
以上所述仅为本发明的较佳实施例而已,并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等,均包含在本发明的保护范围内。
Claims (10)
1.一种应用于电磁散射的h-自适应矩量方法,其特征在于,包括:
创建用于计算待测目标散射特性的电磁模型,进行平面波激励下所述电磁模型的初始网格划分,得到待测目标表面划分出的三角形网格;
利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度;
基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差;
基于各网格单元的局部误差计算全局误差;
基于所述全局误差判断是否满足迭代停止条件;
若否,选择部分网格单元基于B-W插点算法进行网格细化,得到重新划分的网格;并返回所述利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度的步骤;
若是,停止迭代,输出当前求解出的各网格单元的电流系数。
2.根据权利要求1所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述进行平面波激励下所述电磁模型的初始网格划分,包括:
将所述电磁模型以λ/4的尺寸进行初始网格划分;其中,λ表示所述平面波的波长。
3.根据权利要求2所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述利用基于RWG基函数的矩量法,求解当前划分出的各网格单元的电流系数,利用求解出电流系数计算对应网格单元中超收敛点处的面电荷密度,包括:
利用所述基于RWG基函数的矩量法和当前划分的网格对待测目标表面的等效电磁流进行离散,将利用网格离散后产生的基函数组成线性方程组,求解所述线性方程组得到各网格单元的电流系数;
利用各网格单元的电流系数和面电荷密度计算公式,计算出对应网格单元中超收敛点处的面电荷密度。
4.根据权利要求3所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述基于各网格单元中超收敛点处的面电荷密度和超收敛点恢复法,计算对应网格单元的局部误差,包括:
针对属于各网格单元的每一个顶点,利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值;
针对每个网格单元,将该网格单元的三个顶点分别得到的恢复电荷密度值利用最小二乘法进行线性拟合,构造出该网格单元的线性恢复电荷密度函数;
为保证RWG基函数关于三角形对中正负三角形上电荷量的特性,通过将该网格单元的线性恢复电荷密度函数与常数相加的形式,得到该网格单元的修正恢复电荷密度函数;
基于该网格单元的修正恢复电荷密度函数和原超收敛点处的面电荷密度,计算该网格单元的局部误差。
5.根据权利要求4所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述利用包含该顶点的所有网格单元中超收敛点处的面电荷密度,得到该顶点的恢复电荷密度值,包括:
将包含该顶点的所有网格单元中超收敛点处的面电荷密度进行平均,将得到的平均值确定为该顶点的恢复电荷密度值。
6.根据权利要求5所述的应用于电磁散射的h-自适应矩量方法,其特征在于,计算该网格单元的局部误差所采用的公式,包括:
其中,表示第i个网格单元的局部误差;ρi',rec表示第i个网格单元的修正恢复电荷密度函数;ρi表示第i个网格单元中超收敛点处的面电荷密度;Si表示第i个网格单元。
7.根据权利要求6所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述基于各网格单元的局部误差计算全局误差,包括:
针对每个网格单元,将其三角形面积作为自身权重对该网格单元的局部误差进行加权,得到该网格单元的加权局部误差;并求取所有网格单元的加权局部误差之和,得到加权局部误差和;
将所述加权局部误差和利用所有网格单元的三角形面积和进行归一化;
将归一化结果求取平方根,得到全局误差。
8.根据权利要求7所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述基于所述全局误差判断是否满足迭代停止条件,包括:
当所述全局误差降低为预设的初始全局误差的1/Q以下,或者迭代次数达到预设的最大迭代次数,则确定满足迭代停止条件;否则,确定不满足迭代停止条件;其中Q为大于1的自然数。
9.根据权利要求8所述的应用于电磁散射的h-自适应矩量方法,其特征在于,所述选择部分网格单元基于B-W插点算法进行网格细化,包括:
计算所有网格单元中最大的加权局部误差与预设百分比的乘积得到加权局部误差阈值;
确定加权局部误差大于所述加权局部误差阈值的网格单元作为选定网格单元;
利用基于B-W插点算法对各选定网格单元进行网格细化。
10.根据权利要求9所述的应用于电磁散射的h-自适应矩量方法,其特征在于,利用基于B-W插点算法对任一选定网格单元进行网格细化,包括:
在当前Delaunay三角化Mi中搜索到覆盖待插入点pi的基单元;其中,当前Delaunay三角化Mi表示当前划分出的三角形网格;
利用网格单元的相邻关系在基单元的相邻网格单元处标记所有不符合外接圆准则的网格单元,并将这类网格单元的集合记为Bi;
在当前Delaunay三角化Mi中删除Bi所覆盖的网格单元形成空腔Ci,连接pi和Ci边界上的顶点,在空腔内形成新的三角化Bi+1,得到新的Delaunay三角化Mi+1。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310357007.8A CN116449324A (zh) | 2023-04-04 | 2023-04-04 | 一种应用于电磁散射的h-自适应矩量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310357007.8A CN116449324A (zh) | 2023-04-04 | 2023-04-04 | 一种应用于电磁散射的h-自适应矩量方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116449324A true CN116449324A (zh) | 2023-07-18 |
Family
ID=87135022
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310357007.8A Pending CN116449324A (zh) | 2023-04-04 | 2023-04-04 | 一种应用于电磁散射的h-自适应矩量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116449324A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118278227A (zh) * | 2024-06-04 | 2024-07-02 | 西安中电科西电科大雷达技术协同创新研究院有限公司 | 适用于通用误差指示器的高效网格细化自适应矩量方法 |
-
2023
- 2023-04-04 CN CN202310357007.8A patent/CN116449324A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118278227A (zh) * | 2024-06-04 | 2024-07-02 | 西安中电科西电科大雷达技术协同创新研究院有限公司 | 适用于通用误差指示器的高效网格细化自适应矩量方法 |
CN118278227B (zh) * | 2024-06-04 | 2024-08-09 | 西安中电科西电科大雷达技术协同创新研究院有限公司 | 适用于通用误差指示器的高效网格细化自适应矩量方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Tao et al. | KD-tree based fast ray tracing for RCS prediction | |
CN106654566B (zh) | 一种飞行器天线罩的快速厚度设计方法 | |
CN116449324A (zh) | 一种应用于电磁散射的h-自适应矩量方法 | |
CN116071519B (zh) | 基于调和映射生成网格模型的图像处理方法和装置 | |
CN112381862B (zh) | 一种cad模型与三角网格全自动配准方法和装置 | |
CN111276822A (zh) | 一种天线方向图主瓣可控的天线阵列设计方法 | |
CN109887024A (zh) | 一种点云法线估算新方法 | |
CN109241556B (zh) | 一种随机粗糙体目标建模方法及存储介质 | |
CN112733364A (zh) | 一种基于阻抗矩阵分块的箔条云散射快速计算方法 | |
CN114755652A (zh) | 基于aca与cat的获取电大尺寸目标宽带rcs方法 | |
CN113255194A (zh) | 一种基于参数曲面的自主智能直六面体剖分方法 | |
CN111914364B (zh) | 基于高阶矩量法与投影的频选天线罩建模方法 | |
CN113658077A (zh) | 一种基于曲率的分区域双边海量点云降噪方法 | |
CN105869210A (zh) | 三维地质表面模型中的插值数据处理方法 | |
CN116401498A (zh) | 一种基于mb-crwg基函数的电磁散射分析方法、设备及介质 | |
CN115795936A (zh) | 一种用于电磁器件设计的参数-拓扑混合优化方法 | |
Sun et al. | Automatic quadrilateral mesh generation and quality improvement techniques for an improved combination method | |
CN115937460A (zh) | 基于最优传输的保特征表面重建方法 | |
CN115754961A (zh) | 一种基于属性散射中心近场修正模型的回波生成方法 | |
Tian et al. | Accelerated hybrid method for electromagnetic scattering from multiple complex targets above a rough surface | |
CN109636914B (zh) | 一种基于自适应滚球算法的曲面重建方法和系统 | |
Xiang et al. | The application of Barycentric subdivision method for numerical integration in method of moments | |
CN105869209A (zh) | 三维地质表面模型中的畸形三角形数据处理方法 | |
Altinozen et al. | Green coordinates for generation of conformal antenna geometries | |
CN118278227B (zh) | 适用于通用误差指示器的高效网格细化自适应矩量方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |