CN104679919B - 计算微波谐振电路时域响应的外推方法 - Google Patents

Abstract

本发明公开了一种计算微波谐振电路时域响应的外推方法。该方法的实现分为两个过程。第一个过程为预处理过程，对三维微波谐振电路结构建模并使用曲六面体剖分，使用时域谱元法数值分析方法计算出前一段时间的响应，并通过改进的模态叠加方法提取出能够准确描述线性系统的一系列重要特征值，这些特征值对应的特征向量用以构成特征向量矩阵；第二个过程即外推计算过程，利用预处理过程最终得到的特征向量矩阵，对时域谱元法描述的线性系统进行特征提取降阶，重新形成一个只包含重要特征的小规模线性系统，该系统的时域外推求解不仅快速而且准确，能够适用于各种具有谐振、损耗与反射结构电磁响应的外推分析。

Description

计算微波谐振电路时域响应的外推方法

方法领域

本发明属于电磁学仿真计算与信号处理的交叉领域，具体是一种电磁特性的快速仿真与外推计算方法。

背景方法

在信号处理领域，信号的外推方法可以用于实现信号参数的特征估计、雷达目标探测识别、核磁共振快速成像等方法；在数值计算仿真领域，外推方法可以极大的加快计算，有利于时间与内存资源的合理利用。在过去的几十年里，计算电磁学数值分析方法被广泛运用于各种电磁问题的分析当中，对通讯、集成电路设计、射频器件与天线设计、电磁兼容、电磁散射与目标识别等领域的发展起了极大地推动作用，不仅有利于降低实验研究的成本与风险，缩短设计与开发周期，还对理论研究有很大帮助。

国内外对于电磁响应信号参数估计和外推方法的研究已经发展多年，从1971年C.E.Baum提出的奇点展开法，到后来发展起来的Prony算法、KT算法和Matrix-Pencil算法及其它们的各种改进算法，已经可以从信噪比较低的环境中中提取信号特征，这类方法的数据主要来自于实际测量，优点是只需要知道信号的一系列采样数据，即可推测出信号特征，有利于未知源或目标的探测识别，缺点是受采样值影响较大，且外推范围有限；近年来，也有研究者将基于数值仿真的全波分析方法与信号外推方法结合来加快计算过程，采样数据来源于数值仿真计算，主要用于已知结构的响应特性快速分析，例如FDTD与改良的Matrix Pencil方法结合，来分析电路散射特性等，缺点主要是耦合度低，由于所采用的外推方法是一种纯粹的信号处理过程，脱离目标结构本身特性，大多只适用于分析衰减与非衰减的周期信号与指数型叠加信号，适用范围有一定的局限性，且受外推算法的影响，外推范围仍然有限。

发明内容

本发明的目的在于提供一种基于数值分析方法计算微波谐振电路响应的外推方法，该方法提高了数值分析方法与外推方法的耦合度，结合目标本身结构特性外推计算，使得外推范围大幅度提高，计算更加快速，外推结果更加可靠。

实现本发明目的的方法解决方案为：一种计算微波谐振电路时域响应的外推方法，步骤如下：

第一步，采用全波分析的时域谱元法对使用曲六面体剖分的目标模型进行显式格式的迭代求解。

第二步，用一定间隔时间步（采样间隔需要满足奈奎斯特采样定律）采样计算的电场构成一个正交矩阵V，用该矩阵对原问题进行降阶，求降阶系统的特征值问题，每一个采样歩进行上述过程，直到毎个采样步计算得到的重要特征值趋于稳定，预处理结束。

第三步，用重要特征值对应的特征向量构成一个Φ_r矩阵，该矩阵对预处理结束时的降阶系统进行再降阶，即特征提取，可以只保留下原系统的重要特征，构造出用于外推计算的小规模系统，通过迭代求解该系统，可以还原出原系统的解，从而实现快速而准确的外推计算。

本发明与现有方法相比，其显著优点：

(1)所采用的数值方法为时域谱元法，具有谱精度，由于产生的质量矩阵为块对角矩阵，可以快速求逆，而且可以采用显式格式的差分迭代求解，求解速度较常用的时域电磁学分析方法快，质量矩阵与刚度矩阵都满足对称性与正定半正定特性，使得广义特征值的求解收敛速度更快。

(2)采用改进型的模态叠加法进行参数估计，使得提取特征参数时只需要求解小规模的特征值问题，速度快，且降阶后的矩阵特性变好，参数提取更加准确。

(3)提取出来的特征参数用于重构一个小规模线性系统用于外推计算，这种方式提高了外推方法与结构本身特性的耦合度，使得外推的范围大大增加，结果更加可信。

附图说明

图1实例1卫星舱体模型与曲六面体剖分效果示意图。

图2实例1使用外推方法与纯数值方法的时域电场结果对比。

图3实例2含不连续结构的矩形波导模型示意图。

图4实例2观察点1使用外推方法与纯数值方法的时域电场结果对比。

图5实例2观察点2使用外推方法与纯数值方法的时域电场结果对比。

图6实例2使用外推方法与纯数值方法的S参数结果对比。

具体实施方式

下面结合附图对本发明作进一步详细描述。

结合图1，本发明为一种计算微波谐振电路时域响应的外推方法，具体实施步骤如下：

第一步，时域谱元法微分方程的建立。

对目标模型使用ANSYS软件建模，并使用20个节点的曲六面体剖分，得到模型剖分后的单元与节点的信息文件，将每个20点表示的曲六面体映射成一个三轴坐标均从-1到1的标准立方体，这里使用20个节点曲六面体的形函数为：

其中(ξ,η,ζ)表示标准立方体内任意一点坐标，(ξ_i,η_i,ζ_i)是20个剖分点映射到标准参量坐标系下的坐标。

物理坐标系下各结点的坐标(x_i,y_i,z_i)与参量坐标系下各结点的坐标(x′,y′,z′)之间的映射关系如下：

真个模型的求解域为Ω，求解域内电场满足电场时域波动方程：

采用GLL插值基函数展开单元内电场E：

其中表示物理坐标系下的展开基函数，而e_i(t)则为展开系数。

通过映射关系，给出物理坐标系与标准参量坐标系下的基函数关系：

其中，J为雅克比矩阵，为标准参量坐标系下的展开基函数，

其中分别代表参量坐标系下三个方向的基函数，r、s、t则代表标准参量坐标系下某点的坐标。

坐标映射关系雅克比矩阵的形式：

(6)式中的φ_r(ξ)、φ_s(η)、φ_t(ζ)是具有正交性质的Gauss-Lobatto-Legendre基函数，表达式为：

其中为Legendre多项式，为Legendre多项式的一阶导数，N_ξ代表基函数的插值阶数，的微分形式如下：

首先，将波动方程中的电场用(4)式展开，并对该方程使用Galerkin测试，得到：

其中i代表展开基函数索引，j代表测试基函数索引，为单元e上的参数加权积分，伽辽金测试中要求在整个求解域Ω上的残数加权积分为0，即

整理得到带求解的时域谱元法微分方程：

其中，向量为待求解向量，代表参量坐标系下的不同位置和方向的电场，为电流密度，这里作为激励源。

其他矩阵或向量形式为：

[T]矩阵为块对角矩阵，具有很快的求逆速度。

对(11)式采用中心差分格式迭代求解：

求出后，可以通过雅克比矩阵映射回物理坐标系下，得到电场：

第二步，采样与预处理过程。

预处理过程的目的是找到一个包含所有重要模态足以描述系统的特征向量矩阵，即特征的提取。

设每隔k步对(13)式求出的向量进行采样，装入V矩阵，毎采样一次向量，就将其与前一次采样的V矩阵进行正交化，得到一个正交化的V矩阵。将在该矩阵展开，即：

将(15)式代入(11)式，方程两边左乘[V]^T，得到：

其中[V^TTV]、[V^TSV]、[V^TRV]矩阵是规模为k×k的对称矩阵，k＜＜N，N为 向量中元素的个数。方程(11)的阶数被降低到了k阶，该降阶方程对应的特征多项式方程为：

{λ²[V^TTV]+λ[V^TRV]+[V^TSV]}φ＝0. (19)

求解(19)式所示的特征值方程，得到k对特征值λ与特征向量

在该预处理过程中，每次采样构成正交化的V矩阵，求解(19)式得到k个特征值及每个特征值对应的特征向量。

设置一个条件：

其中λ_i代表第i个特征值，i＝1,2...k，n代表前一次采样，n+1代表当前采样，ε₁是人为设定的一个阈值，通常取0.001或者0.01，该条件称之为稳定性条件，满足该条件的特征值为稳定的重要特征值。

再设定一个重要特征值个数M，只要满足稳定性条件的特征值个数r大于等于M，则停止预处理过程，此时r个稳定的重要特征值对应的特征向量可以构成Φ_r矩阵，使用Φ_r矩阵来表征(11)式系统的重要特征。

第三步，外推求解小规模系统。

预处理停止后，即找到了(11)式系统的所有重要模态，设重要模态的个数为r,将重要特征值对应的特征向量构成一个特征提取矩阵Φ,将(16)式的场解 在该矩阵展开：

其中k为(16)V矩阵的列数，r为重要模态的个数。

将(20)式代回(15)式，可得：

再将(21)式代回(11)式，可以得到小规模线性系统：

该系统的规模仅为k，外推求解就是对(22)式进行迭代求解，原矩阵的对称和正定特性得以保留，求解速度大大快于(11)式所示的规模为N的系统，求出 后，可以通过(21)式还原出原问题的解

下面用2个实例来验证该外推方法的准确性和快速特性。

第一个实例为一个边长为2.7m的金属卫星空舱，如图1所示，正中有一个内径0.89m外径0.9m的有耗薄介质圆筒，高度为2.7m，相对介电常数为2.0，电导率为2.0。在正中心处模拟一个点源激励，激励形式采用正弦调制高斯形式，中心频率1.0GHz，频带范围从0.5GHz-1.5GHz,计算未知量为94851。在圆筒中轴线上离中心1.05m处设定一个观察点，预处理过程V矩阵的构造使用了前20ns的时域电场值毎0.4ns采样得到，稳定性判定的阈值精度使用0.001，外推结果计算到80ns，结果如图2所示，所用的时间和加速比数据如表1所示，该结构实际上是一个复杂谐振结构的信道问题，激励源在谐振结构内部多次衰减、反射、叠加，需要保留的模态较多，这里保留了71个模态。

表1-实例1计算时间与加速比对比

第二个实例为一个矩形波导，如图3所示，波导横截面尺寸20mm×10mm，中间有一介质块作为不连续结构，尺寸为10mm×10mm×20mm，相对介电常数为6.0，波导两端加了PML边界条件，在波导一侧加面激励源，采用正弦调制高斯信号形式，中心频率为10GHz，频带范围从8.5GHz-11.5GHz。计算未知量13062个未知量，时间迭代步数为3000步，预处理消耗450步，稳定性判定的阈值精度使用0.001，最终保留43个模态。在不连续结构两侧分别设置一观察点，用来观察时域电场,分别如图4和5所示，并用来计算S参数，如图6所示，所用的时间和加速比数据如表2所示，可见计算过程更快，外推结果准确可靠。

表2-实例2计算时间与加速比对比

注1：总加速比=(SETD总时间)/(MSM_SETD时间)

迭代加速比=(SETD总时间)/(MSM_SETD时间-MSM_SETD预处理时间)

注2：SETD-Spectral Element Method of Time Domain，时域谱元法

MSM-Modal Superposition Method，模态叠加法。

Claims (4)

1.一种计算微波谐振电路时域响应的外推方法，其特征在于，步骤如下：

第一步，构造一个待求解的时域谱元法线性系统方程；

对物理坐标系下的三维微波谐振类型的结构进行建模并采用曲六面体剖分，得到模型剖分后的单元与节点信息，将物理坐标系下得到的曲六面体映射为三边坐标均为-1到1的标准立方体，映射后的坐标系称为标准参量坐标系；

整个模型为求解域Ω，求解域Ω内的电场满足电场时域波动方程：

其中ε为求解域内的介电常数，μ为磁导率，σ为电导率，为电流密度；

对于方程(1)，采用Gauss-Lobatto-Legendre插值基函数展开物理坐标系单元内电场

其中表示物理坐标系中的展开基函数，e_i(t)为展开系数，N表示单元内展开基函数的个数；

然后，对方程(1)使用Galerkin测试，得到：

其中i代表展开基函数索引，j代表测试基函数索引，V^e表示单元e内，为单元e上的残数加权积分，Galerkin测试要求在整个求解域Ω上的残数加权积分为0，即必须满足

最后，整理(3)式得到待求解的时域谱元法线性系统方程：

其中，向量为待求解向量，代表标准参量坐标系下的不同位置和方向的电场在某时刻的值，为电流密度，这里作为激励源；

(4)式中其他矩阵或向量形式为：

采用中心差分格式对方程(4)迭代求解：

求出标准参量坐标系下的电场后，通过雅克比矩阵[J]映射回物理坐标系，得到物理坐标系下电场：

至此，完成了(4)式所示的时域谱元法线性系统方程的构造；

第二步，进行预处理得到重要特征构成的特征向量矩阵Φ_r；

首先，创建一个空的V矩阵，行为向量的元素个数，初始列数为0；

对第一步构造的系统方程进行(5)式所示的时间步差分迭代求解，每一定时间间隔对进行采样，将采样的向量追加为V矩阵的一列向量，对V矩阵正交化,V矩阵的列数为k；

将用V矩阵展开，将其代入(4)式，方程两边同时左乘[V]^T，得到：

其中[V^TTV]、[V^TSV]、[V^TRV]矩阵是规模为k×k的对称矩阵，k＜＜N，N为向量中元素的个数；方程(6)的阶数被降低到了k阶，该降阶方程对应的特征多项式方程为：

求解(7)式所示的特征值方程，得到k对特征值λ与特征向量根据重要特征值的判定准则来判断是否停止预处理，预处理停止后，取出能够表征重要特征的r个特征向量构成Φ_r矩阵；

第三步，外推计算小规模系统；

用第二步最终构造的Φ_r矩阵，展开方程(6)的解

其中k为(6)式中V矩阵的列数,r为重要特征值的个数，x_i为每个特征向量的展开系数；

将(8)式代入(6)式，方程两边同时左乘Φ_r ^T，得到：

该(9)式所示的系统规模仅为r×r，r＜k＜＜N；

迭代求解(9)式所示的小规模系统，得到再根据

还原出系统(4)的解从而实现快速准确的外推计算。

2.根据权利要求1所述的计算微波谐振电路时域响应的外推方法，其特征在于：所述步骤一中，物理坐标系下得到的曲六面体映射为三边坐标均为-1到1的标准立方体的步骤如下：

使用20个节点的曲六面体形函数P_i表示一个曲六面体：

其中(ξ,η,ζ)表示标准参量坐标系下标准立方体内任意一点的坐标，(ξ_i,η_i,ζ_i)是20个曲六面体节点映射到标准参量坐标系下的坐标；

标准参量坐标系下任意一点的坐标(ξ,η,ζ)对应的物理坐标系下的坐标为(x,y,z)，映射关系如下：

采用雅克比矩阵表示点(ξ,η,ζ)与(x,y,z)之间的映射关系：

。

3.根据权利要求1所述的计算微波谐振电路时域响应的外推方法，其特征在于：所述步骤一中，采用Gauss-Lobatto-Legendre插值基函数的步骤如下：

设为标准参量坐标系下的展开基函数，

通过映射关系，物理坐标系下基函数与标准参量坐标系下的基函数之间的关系表示为：

其中分别代表标准参量坐标系下三个方向ξ、η、的基函数，r、s、t则代表标准参量坐标系下某点的坐标，|J|为雅克比矩阵；φ_r(ξ)、φ_s(η)、φ_t(ζ)是具有正交性质的Gauss-Lobatto-Legendre基函数，表达式为：

其中为Legendre多项式，为Legendre多项式的一阶导数，N_ξ代表基函数的插值阶数，的微分形式如下：

。

4.根据权利要求1所述的计算微波谐振电路时域响应的外推方法，其特征在于：所述步骤二中的判定准则如下：

对于权利要求1中第二步的预处理过程，每次采样构成正交化的V矩阵，求解(7)式得到k个特征值及每个特征值对应的特征向量；

设置一个条件：

其中λ_i代表第i个特征值，i＝1,2...k，n代表前一次采样，n+1代表当前采样，ε₁是设定的一个阈值，取0.001或者0.01，该条件称之为稳定性条件，满足该条件的特征值为稳定的重要特征值；

再设定一个重要特征值个数M，只要满足稳定性条件的特征值个数r大于等于M，则停止预处理过程，此时r个特征值对应的特征向量构成矩阵Φ_r。