CN105929253A - 采用偶极矩模型对pcb电路电磁场进行反演的方法 - Google Patents

Abstract

一种采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于包括有下列步骤：一：利用电磁场探头测试近场电磁场强度数据；二：设置偶极子阵列；三：计算采样点和偶极子阵列之间的映射矩阵T；四：计算偶极矩矩阵X(α)；五：计算偶极矩矩阵X(α)在高于PCB电路任意高度观察平面上的切向电磁场分量；六：把步骤五求得的切向电磁场分量代入公式(32)和(33)进行计算，求解步骤五中平面上的法向电磁场分量。优点在于：该矩阵能够很好的反映PCB电路的电压和电流分布，同时可以对电磁场的辐射进行准确的计算，降低了电磁场分布的测试成本；为系统级电磁兼容设计指标量化验证环节提供了有力依据；实现了PCB电路电磁场的反演，得到PCB电路的近场耦合和远场辐射特性。

Description

采用偶极矩模型对PCB电路电磁场进行反演的方法

【技术领域】

本发明涉及一种对电磁场分布进行反演的方法，具体来说是对未知悉物理特征的PCB电路利用偶极矩模型进行等效重建，反演PCB电路电磁场近场耦合和远场辐射的一种方法。

【背景技术】

随着电气工程、电子科学技术、计算机技术、控制理论与控制工程等技术的飞速发展，PCB电路在各行各业的应用占据不可或缺的地位。PCB电路工作频率频段不断扩展，集成度和复杂度越来越高，PCB电路中的任何一段走线都有可能成为辐射电磁场的源，不同的PCB电路辐射之间的耦合，是造成系统级电磁兼容问题的主要原因之一。

国家技术监督局发布的文件明确指出，任何不符合国家标准规定，无线电干扰严重的产品禁止生产和使用。对不同PCB电路的辐射进行准确的测量是电磁兼容设计中指标量化环节有效的考量手段。但是在对PCB电路辐射的电磁场强度进行测量时，并不能获得空间中每一点处的辐射强度，这就有可能在系统组装完毕后形成预料外的电磁兼容性问题。同时另一个造成电磁兼容设计存在缺陷主要原因是PCB电路的具体物理特征受到各种原因有不能被使用者知悉，因此不能通过前期的仿真工作对其辐射进行准确的计算，在系统整合过程中形成了潜在的电磁干扰。

【发明内容】

为了解决PCB电路板产生的电磁场分布难以获取的问题，本发明提出一种采用偶极矩模型对PCB电路电磁场进行反演的方法。该方法通过测量空间中有限点处的近场电场强度和磁场强度形成辐射强度数值矩阵F，设置偶极矩阵列与辐射强度数值矩阵之间的映射矩阵T，求解得到偶极矩矩阵X，继而对PCB电路的电磁场分布进行反演。

本发明中由电磁场理论可知，任意的电小源都可以由6个偶极矩分量等效代替：三个电偶极矩Px，Py，Pz和三个磁偶极矩分量Mx，My，Mz。这些偶极子组合起来产生的场与原来的电小源产生的场相同，因此利用偶极矩模型反演PCB电路的电磁场是可行的。为了减小电路中的回路电感，PCB电路的设计中走线紧贴电源平面或者参考地平面，这样用来等效PCB电路的偶极矩在满足理想导体的边界条件前提下减少为3个：Mx，My和Pz。垂直电偶极子Pz可以描述PCB走线与参考地平面之间的电压分布，水平磁偶极子Mx和My可以描述PCB走线上的电流分布。用N×N个偶极子的偶极子阵列来等效PCB电路，每个偶极子包含三个偶极矩分量，分别为Pz，Mx，My。在M×M个采样数据点的近场采样平面上，每个采样数据点的水平场强分别为Ex，Ey，Hx，Hy。设定偶极子的位置坐标为(x’,y’,d),采样点的位置坐标为(x,y,z)，d表示的是偶极子距离参考地平面的高度，z表示采样点距离参考地平面的高度，其中z>d>0，采样点(x,y,z)处的切向电场强度和磁场强度可以由以下公式计算出来

$E_x={\mathrm{\τ}}_E\{\[\frac{(z-d)(x-x^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(z+d)(x-x^{\′})}{r_2^2}{q}_1\left(r_2\right)\]P_z+\[\frac{(z-d)}{r_1}{q}_3\left(r_1\right)+\frac{(z+d)}{r_2}{q}_3\left(r_2\right)k_0{M}_y\]\}---\left(1\right)$

$E_y={\mathrm{\τ}}_E\{\[\frac{(y-y^{\′})(z-d)}{r_1^2}{q}_1\left(r_1\right)+\frac{(y-y^{\′})(z+d)}{r_2^2}{q}_1\left(r_2\right)\]P_z+\[-\frac{(z-d)}{r_1}{q}_3\left(r_1\right)-\frac{(z+d)}{r_2}{q}_3\left(r_2\right)\]k_0{M}_x\}---\left(2\right)$

$\begin{array}{c}{H}_x={\mathrm{\τ}}_H\{-\[\frac{(y-y^{\′})}{r_1}{q}_3\left(r_1\right)+\frac{(y-y^{\′})}{r_2}{q}_3\left(r_2\right)\]P_z+\[-\frac{{(y-y^{\′})}^2+{(z-d)}^2}{r_1^2}{q}_1\left(r_1\right)+q_2\left(r_1\right)\\ -\frac{{(y-y^{\′})}^2+{(z+d)}^2}{r_2^2}{q}_1\left(r_2\right)+q_2\left(r_2\right)\]k_0{M}_x+\[\frac{(x-x^{\′})(y-y^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(x-x^{\′})(y-y^{\′})}{r_2^2}{q}_1\left(r_2\right)\]k_0{M}_y\}\end{array}---\left(3\right)$

$\begin{array}{c}{H}_y={\mathrm{\τ}}_H\{\[\frac{(x-x^{\′})}{r_1}{q}_3\left(r_1\right)+\frac{(x-x^{\′})}{r_2}{q}_3\left(r_2\right)\]P_z+\[\frac{(x-x^{\′})(y-y^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(x-x^{\′})(y-y^{\′})}{r_2^2}{q}_1\left(r_2\right)\]k_0{M}_x\\ +\[-\frac{{(x-x^{\′})}^2+{(z-d)}^2}{r_1^2}{q}_1\left(r_1\right)+q_2\left(r_1\right)-\frac{{(x-x^{\′})}^2+{(z+d)}^2}{r_2^2}{q}_1\left(r_2\right)+q_2\left(r_2\right)\]k_0{M}_y\end{array}---\left(4\right)$

在以上公式中，j是虚数单位，k₀和η₀分别表示自由空间中的波数和波阻抗，f为频率，c为真空中的光速，η₀＝120π，Pz表示垂直电偶极矩，它是一个复数，单位为A·m，Mx和My表示水平磁偶极矩，单位为A·m²，r₁表示近场采样数据点与偶极子阵列中每个偶极子的距离，r₂表示近场采样数据点与偶极子阵列中每个偶极子的镜像的距离。

根据镜像原理，面对较难求解的给定系统，可以构造和希望求解系统的电位场解在原求解区域相等的一组适当配置的电荷系统来等效替代，从而简化求解过程。因此，可以构造虚拟的镜像偶极子来求解靠近参考地平面的偶极子的辐射场。

r₁和r₂可由下式表示：

r₁＝[(x-x′)²+(y-y′)²+(z-d)²]^1/2 (5)

r₂＝[(x-x′)²+(y-y′)²+(z+d)²]^1/2 (6)

q₁(r)、q₂(r)、q₃(r)可以表示为:q₁(r)、q₂(r)、q₃(r)没有物理意义，这是为了简化公式(1)-(4)的，否则公式显得太复杂；

$q_1\left(r\right)=\[\frac{3}{{\left(k_{0}r\right)}^2}+j\frac{3}{k_{0}r}-1\]f\left(r\right)---\left(7\right)$

$q_2\left(r\right)=\[\frac{2}{{\left(k_{0}r\right)}^2}+j\frac{2}{k_{0}r}\]f\left(r\right)---\left(8\right)$

$q_3\left(r\right)=\[\frac{1}{k_{0}r}+j\]f\left(r\right)---\left(9\right)$

其中，r为空间点的矢径的模值，e为自然常数；

因此用来等效PCB电路的偶极子的偶极矩(Pz,Mx,My)和近场采样数据场强(Ex,Ey,Hx,Hy)存在下述的映射关系：

$F=\left(\begin{array}{c}{\[Ex\]}_{M^2\×1}\\ {\[Ey\]}_{M^2\×1}\\ {\[Hx\]}_{M^2\×1}\\ {\[Hy\]}_{M^2\×1}\end{array}\right)=TX=T\left(\begin{array}{c}{\[Pz\]}_{N^2\×1}\\ {\[Mx\]}_{N^2\×1}\\ {\[My\]}_{N^2\×1}\end{array}\right)---\left(10\right)$

其中[Ex]，[Ey]，[Hx]，[Hy]分别表示x方向和y方向上的电场强度和磁场强度矩阵，近场采样点的个数为M×M。是待求的偶极矩矩阵，[Pz]是偶极子阵列的垂直电偶极矩矩阵，[Mx]、[My]是偶极子阵列的水平磁偶极矩，每种偶极子个数为N×N。映射矩阵描述的是空间中任意一点电磁场与偶极子矩阵之间关系，映射矩阵可由下式表示：

$T=\left(\begin{array}{ccc}{T}_{ExPz}& T_{ExMx}& T_{ExMy}\\ T_{EyPz}& T_{EyMx}& T_{EyMy}\\ T_{HxPz}& T_{HxMx}& T_{HxMy}\\ T_{HyPz}& T_{HyMx}& T_{HyMy}\end{array}\right)---\left(11\right)$

上式每一个子矩阵的维数为M²×N²，表示的是由某一种偶极矩产生的某一方向上的场强。以T_ExMy为例，它表示所有y方向上的磁偶极矩My在x方向上产生的电场强度Ex。T矩阵总共有4M²×3N²个元素，每一个子矩阵可由下式表示：

$T_{EyPz}(a,b)={\mathrm{\τ}}_E\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0-d)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0+d)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]---\left(12\right)$

T_ExMx(a,b)＝0 (13)

$T_{ExMy}(a,b)={\mathrm{\τ}}_E\[\frac{(z_0-d)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{(z_0+d)}{r_{+}}{q}_3\left(r_{+}\right)\]k_0---\left(14\right)$

$T_{EyPz}(a,b)={\mathrm{\τ}}_E\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0-d)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0+d)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]---\left(15\right)$

$T_{EyMx}(a,b)={\mathrm{\τ}}_E\[-\frac{(z_0-d)}{r_{-}}{q}_3\left(r_{-}\right)-\frac{(z_0+d)}{r_{+}}{q}_3\left(r_{+}\right)\]k_0---\left(16\right)$

T_EyMy(a,b)＝0 (17)

$T_{HxPz}(a,b)=-{\mathrm{\τ}}_H\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)}{r_{+}}{q}_3\left(r_{+}\right)\]---\left(18\right)$

$\begin{array}{c}{T}_{HxMx}(a,b)={\mathrm{\τ}}_H\[-\frac{{\left(y\right(a)-y^{\′}(b\left)\right)}^2+{(z_0-d)}^2}{{r_{-}}^2}{q}_1\left(r_{-}\right)+q_2\left(r_{-}\right)\\ -\frac{{\left(y\right(a)-y^{\′}(b\left)\right)}^2+{(z_0+d)}^2}{{r_{+}}^2}{q}_1\left(r_{+}\right)+q_2\left(r_{+}\right)\]k_0\end{array}---\left(19\right)$

$T_{HxMy}(a,b)={\mathrm{\τ}}_H\left[\begin{array}{c}\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\\ \frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\end{array}\right]k_0---\left(20\right)$

$T_{HyPz}(a,b)={\mathrm{\τ}}_H\[\frac{\left(x\right(a)-x^{\′}(b\left)\right)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{\left(x\right(a)-x^{\′}(b\left)\right)}{r_{+}}{q}_3\left(r_{+}\right)\]---\left(21\right)$

$T_{HyMx}(a,b)={\mathrm{\τ}}_H\[\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]k_0---\left(22\right)$

$\begin{array}{c}{T}_{HyMy}(a,b)={\mathrm{\τ}}_H\[-\frac{(z_0-d)+{\left(x\right(a)-x^{\′}(b\left)\right)}^2}{{r_{-}}^2}{q}_1\left(r_{-}\right)+q_2\left(r_{-}\right)\\ -\frac{{(z_0+d)}^2+{\left(x\right(a)-x^{\′}(b\left)\right)}^2}{{r_{+}}^2}{q}_1\left(r_{+}\right)+q_2\left(r_{+}\right)\]k_0\end{array}---\left(23\right)$

上述公式中(x(a),y(a)，z₀)为控件中任意一点的坐标，(x′(b),y′(b),d)表示偶极子的坐标，其中a和b分别表示空间中任意一点和偶极子的序号，并且满足以下约束条件a＝1,2,3...M²，b＝1,2,3...N²。r_-表示空间中任意一点与偶极子的距离，r₊表示任意控件一点与偶极子镜像的距离。

根据公式(5)和(6)可以得到采样点与偶极子之间响应的位置关系，将每一个采样点的坐标(x,y,z)用(x(a),y(a),z₀)表示，每一个偶极子的坐标(x',y',d)用(x'(b),y'(b),d)表示，因此可以得到近场采样数据点与偶极子的距离为：

r₁(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀-d)²]^1/2 (24)

近场采样数据点与偶极子的镜像的距离为：

r₂(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀+d)²]^1/2 (25)

将采样点和偶极子及其镜像之间相应位置关系r₁(a,b)和r₂(a,b)分别用r_-和r₊表示，代入公式(12)-(23)计算得到采样点与偶极子之间的映射矩阵T。

通过电磁场探头测得PCB电路的电磁场分布，通过数值计算求解偶极子阵列作为等效的辐射源是典型的反问题模型。反问题的重要特征就是病态，即使采样数据极小的误差都会引起近似解与真实解之间极大的误差。因此为了能够减弱反问题病态特性的影响，采用正则化方法来求解偶极子的偶极矩。

argmin{||F-TX||²+α||X||²} (26)

X(α)＝[T′T+α²I]^-1T′F (27)

式中“argmin”表示使得(26)取得最小值，α被称作正则化系数。在本发明中采用通过对数坐标尺度描述||X_α||和||F-TX_α||的L曲线准则，α对应L曲线上曲率最大的点。将求解得到的正则化参数α代入(27)可以求得与干扰不敏感的偶极矩矩阵X(α)。|| ||表示一个矩阵的范数，F表示近场采样数据场强矩阵。

在任意观察平面上的设置M×M观察点，x方向上等间隔取M个观察点，y方向上等间隔取M个观察点，每一个观察点的坐标为(x”,y”,z”)，z”是任意观察平面与参考地平面之间的高度距离，且满足z”＞d＞0。将每一个观察的坐标(x”,y”,z”)用(x(a),y(a),z₀)表示，每一个偶极子的坐标(x',y',d)用(x'(b),y'(b),d)表示，因此可以得到观察点与偶极子的距离为：

r₁'(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀-d)²]^1/2 (28)

观察点与偶极子的镜像的距离为：

r₂'(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀+d)²]^1/2 (29)

将观察点和偶极子及其镜像之间相应位置关系r₁'(a,b)和r₂'(a,b)分别用r_-和r₊表示，代入公式(12)-(23)计算得到采样点与偶极子之间的映射矩阵T_v。

通过等效偶极矩矩阵X(α)和映射矩阵T_v可以根据公式(30)求解任意观察平面的电磁场切向分量。

$T_{v}X\left(\mathrm{\α}\right)=F_v=\left[\begin{array}{c}E{x}_v\\ {\mathrm{Ey}}_v\\ {\mathrm{Hx}}_v\\ {\mathrm{Hy}}_v\end{array}\right]---\left(30\right)$

其中F_v是在观察平面上电磁场强度矩阵，该矩阵有4M²×1个元素。其中Ex_v，Ey_v，Hx_v，Hy_v分别表示x方向和y方向上的电场强度和磁场强度，每个矩阵有M²×1个元素。

根据简单无耗媒质中的麦克斯韦方程组可以得到

▽×E(r)＝-jωμH(r) (31)

$\▿\×\tilde{H}\left(r\right)=j\mathrm{\ω}\mathrm{\ϵ}\tilde{E}\left(r\right)---\left(32\right)$

式中E(r),H(r)分别表示电场和磁场的复数形式；▽×E(r)表示E(r)的旋度，▽×H(r)表示H(r)的旋度，ω为角频率，μ磁导率，ε为介电常数。把(31)、(32)展开，由z方向上的分量相等可以得到

其中Hz是观察平面上法向磁场分量，Ez是观察平面上法向电场分量。利用中心差分商代替偏导数，得到(33)、(34)的离散形式：

$Hz(m,n)=\frac{j}{2\mathrm{\ω}\mathrm{\μ}}\{\frac{\[{\mathrm{Ey}}_v(m+1,n)-{\mathrm{Ey}}_v(m-1,n)\]}{\mathrm{\Δ}x}-\frac{\[{\mathrm{Ex}}_v(m,n+1)-{\mathrm{Ex}}_v(m,n-1)\]}{\mathrm{\Δ}y}\}---\left(35\right)$

$Ez(m,n)=-\frac{j}{2\mathrm{\ω}\mathrm{\ϵ}}\{\frac{\[{\mathrm{Hy}}_v(m+1,m)-{\mathrm{Hy}}_v(m-1,n)\]}{\mathrm{\Δ}x}-\frac{\[{\mathrm{Hx}}_v(m,n+1)-{\mathrm{Hx}}_v(m,n-1)\]}{\mathrm{\Δ}y}\}---\left(36\right)$

其中m,n表示每一个观察点在x方向和y方向上的序号，且满足2≤m,n≤M-1，Δx、Δy分别表示x方向和y方向上近场采样点的间距，实现了求解任意观察平面上的电磁场的分布。

本发明一种采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于包括有下列步骤：

步骤一：利用电磁场探头测试近场电磁场强度数据；

步骤二：设置偶极子阵列；

步骤三：计算采样点和偶极子阵列之间的映射矩阵T；

步骤四：计算偶极矩矩阵X(α)；

步骤五：计算偶极矩矩阵X(α)在高于PCB电路任意高度观察平面上的切向电磁场分量；

步骤六：把步骤五求得的切向电磁场分量代入公式(32)和(33)进行计算，求解步骤五中平面上的法向电磁场分量。

在步骤一中，对PCB电路进行辐射强度测量。首先确定PCB电路的相对位置坐标，如图2所示。在PCB电路上方高度为d的采样层平面，均匀测量M×M个点的场强度，沿x方向等间距取M个采样点，采样间隔为ΔX，沿y方向等间距取M采样点，采样间隔为ΔY，其中ΔX＝ΔY。记录每一个采样点坐标为(x，y，z)，将采集到的x方向的电场强度记为E_x，y方向的电场强度记为E_y，x方向的磁场强度记为H_x，y方向的磁场强度记为H_y，每一个矩阵有M²×1个元素，场强归一化结果参见图3，平面场强分布结果参见图3A。将矩阵E_x、E_y、H_x和H_y组合成一个矩阵F，其中

该矩阵有4M²×1个元素。

在步骤二中，参见图2，在PCB上方高度为d的平面上均匀放置N×N个偶极子，沿x方向等间距放置N个偶极子，沿y方向等间距放置N个偶极子，记录每一个偶极子坐标(x’,y’,d)，将步骤一中的采样点的坐标(x，y，z)和步骤二中的偶极子坐标(x’,y’,d)代入公式(24)和(25)中计算得到采样点和偶极子及其镜像的相对位置关系；

在步骤三中，通过步骤二求解采样点和偶极子及其镜像的相对位置关系，利用公式(12)-(25)可以求解出采样点和偶极子之间的映射矩阵T，该矩阵具有4M²×3N²个元素。

在步骤四中，利用步骤一和步骤三求解得到的电磁场强度矩阵F和映射矩阵T，通过公式(10)可知X(α)为待求的偶极矩矩阵，该矩阵有3N²×1个元素。利用正则化方法中的L曲线准则计算正则化系数α，将α代入公式(27)求解偶极矩矩阵X(α)。

在步骤五中，任意观察平面上的设置M×M观察点，x方向上等间隔取M个观察点，y方向上等间隔取M个观察点，每一个观察点的坐标为(x”,y”,z”)，z”是任意观察平面与参考地平面之间的高度距离，且满足z”＞d＞0。利用公式(12)-(23)，(28)和(29)求解得到观察点和偶极子之间的映射矩阵T_v。利用MATLAB软件计算公式(30)可以得到观察平面的电磁场切向分量。

在步骤六中，在任意观察平面上利用步骤五中求解得到切向分量的电场磁场强度Ex_v,Ey_v,Hx_v,Hy_v代入公式(35)和(36)进行计算，求解观察平面上的法向电磁场分量。

本发明采用偶极矩模型对PCB电路电磁场进行反演的方法优点在于：

(1)对于未知物理特性的PCB电路，通过有限个近场电磁场强度的采样数据与偶极子阵列之间的关系，构造出一个等效的偶极矩矩阵，该矩阵能够很好的反映PCB电路的电压和电流分布，同时可以对电磁场的辐射进行准确的计算，大大降低了电磁场分布的测试成本。

(2)相对于传统的近场扫描而言，此发明不仅仅得到PCB电路在近场采样平面的场分布情况，同时得到PCB电路在其他观察平面上的电磁场分布情况，为系统级电磁兼容设计指标量化验证环节提供了有力依据。

(3)利用等效偶极矩矩阵求解观察平面的电磁场切向分量，根据麦克斯韦方程组中心差分算法可以求得该观察平面上电磁场的法向分量，实现了PCB电路电磁场的反演，通过设置不同的观察平面可以得到PCB电路的近场耦合和远场辐射特性。

【附图说明】

图1是本发明的系统配置图。

图2是偶极子阵列构成的等效源模型。

图3是采样层归一化场强值曲线图。

图3A是采样层归一化场强值平面分布图。

图4是采样层与观察层重合时传统计算和本发明计算的归一化场强值曲线图。

图5是采样层与观察层不重合时传统计算和本发明计算的归一化场强值曲线图。

图5A是采样层与观察层不重合时本发明计算归一化场强值平面分布图。

图5B是采样层与观察层不重合时传统计算归一化场强值平面分布图。

图6是本发明操作流程图。

【具体实施方式】

下面将结合附图对本发明做法进一步的详细说明。

参见图1所示的近场扫描平台，该平台包括有测试夹具、电场探头、磁场探头、频谱分析仪和计算机。测试夹具夹持电场探头或磁场探头，电场探头和磁场探头通过导线与频谱分析仪连接。三者的连接是保证近场扫描平台的正常工作。

测试夹具控制探头的测量位置和精度。

电场探头和磁场探头用于测量PCB电路辐射的电磁场分布信息。

频谱分析仪用于将探头获取的场强度信息进行数值显示和存储。

计算机中的MATLAB软件(R2015b)对频谱仪获取的数据进行处理。

参见图6，本发明采用偶极矩模型对PCB电路辐射的电磁场进行反演的方法包括有下列步骤。

第一步：利用电磁场探头测试近场电磁场强度数据。

对PCB电路进行辐射强度测量。首先确定PCB电路的相对位置坐标，如图2所示。在PCB电路上方高度为d的采样层平面，均匀测量M×M个点的场强度，沿x方向等间距取M个采样点，采样间隔为ΔX，沿y方向等间距取M采样点，采样间隔为ΔY，其中ΔX＝ΔY。记录每一个采样点坐标为(x，y，z)，将采集到的x方向的电场强度记为E_x，y方向的电场强度记为E_y，x方向的磁场强度记为H_x，y方向的磁场强度记为H_y，每一个矩阵有M²×1个元素，场强归一化结果参见图3，平面场强分布结果参见图3A。将矩阵E_x、E_y、H_x和H_y组合成一个矩阵F，其中

该矩阵有4M²×1个元素。

第二步：设置偶极子阵列；

在步骤二中，参见图2，在PCB上方高度为d的平面上均匀放置N×N个偶极子，沿x方向等间距放置N个偶极子，沿y方向等间距放置N个偶极子，记录每一个偶极子坐标(x’,y’,d)，将步骤一中的采样点的坐标(x，y，z)和步骤二中的偶极子坐标(x’,y’,d)代入公式(24)和(25)中计算得到采样点和偶极子及其镜像的相对位置关系；

第三步：计算映射矩阵T；

在步骤三中，通过步骤二求解采样点和偶极子及其镜像的相对位置关系，利用公式(12)-(25)可以求解出采样点和偶极子之间的映射矩阵T，该矩阵具有4M²×3N²个元素。

第四步：计算偶极矩矩阵X；

在步骤四中，利用步骤一和步骤三求解得到的电磁场强度矩阵F和映射矩阵T，通过公式(10)可知X(α)为待求的偶极矩矩阵，该矩阵有3N²×1个元素。利用正则化方法中的L曲线准则计算正则化系数α，将α代入公式(27)求解偶极矩矩阵X(α)。

第五步：计算偶极矩矩阵X(α)在高于PCB电路任意高度观察平面上的切向电磁场分量；

在步骤五中，任意观察平面上的设置M×M观察点，x方向上等间隔取M个观察点，y方向上等间隔取M个观察点，每一个观察点的坐标为(x”,y”,z”)，z”是任意观察平面与参考地平面之间的高度距离，且满足z”＞d＞0。利用公式(12)-(23)，(28)和(29)求解得到观察点和偶极子之间的映射矩阵T_v。利用MATLAB软件计算公式(30)可以得到观察平面的电磁场切向分量。

第六步：把步骤五求得的切向电磁场分量代入公式(32)和(33)进行计算，求解步骤五中平面上的法向电磁场分量。

在步骤六中，在任意观察平面上利用步骤五中求解得到切向分量的电场磁场强度Ex_v,Ey_v,Hx_v,Hy_v代入公式(35)和(36)进行计算，求解观察平面上的法向电磁场分量。

实施例

利用电场和磁场探头在采样层采集PCB电路近场电场和磁场的强度，通过设置偶极子阵列求解映射矩阵T，利用正则化算法对偶极矩模型矩阵进行求解。通过改变不同的观察层得到不同的映射矩阵，进而计算得到观察层上的场强分布。通过计算数值和仿真数值计算相对误差，传统计算方法和本发明计算方法在不同观察层上产生的相对误差如表1所示。

表1 实施例中传统计算方法和本发明计算方法在不同观察层上产生的相对误差

根据表1中的数据可以看出采用本发明计算方法和采用传统计算方法在观察层与采样层重叠时计算结果是一致的，参见图4。但是传统计算方法当观察层面与采样层重叠时并不能求解准确的电磁场分布，归一化场强值(观察层3)参见图5，本发明计算平面(观察层3)场强分布结果参见图5A和传统计算平面(观察层3)场强分布结果参见图5B。本发明对PCB电路电磁场的反演是准确的。

Claims (7)

1.一种采用偶极矩模型对PCB电路电磁场进行反演的方法，其中，任意的电小源都由6个偶极矩分量等效代替：三个电偶极矩Px，Py，Pz和三个磁偶极矩分量Mx，My，Mz；

PCB电路的设计中走线紧贴电源平面或者参考地平面，用来等效PCB电路的偶极矩在满足理想导体的边界条件前提下减少为Mx，My和Pz；垂直电偶极子Pz描述PCB走线与参考地平面之间的电压分布，水平磁偶极子Mx和My描述PCB走线上的电流分布；用N×N个偶极子的偶极子阵列来等效PCB电路，每个偶极子包含三个偶极矩分量，分别为Pz，Mx，My；在M×M个采样数据点的近场采样平面上，每个采样数据点的水平场强分别为Ex，Ey，Hx，Hy；设定偶极子的位置坐标为(x’,y’,d),采样点的位置坐标为(x,y,z)，d表示的是偶极子距离参考地平面的高度，z表示采样点距离参考地平面的高度，其中z>d>0，采样点(x,y,z)处的切向电场强度和磁场强度由以下公式计算出来

$E_x={\mathrm{\τ}}_E\{\[\frac{(z-d)(x-x^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(z+d)(x-x^{\′})}{r_2^2}{q}_1\left(r_2\right)\]P_z+\[\frac{(z-d)}{r_1}{q}_3\left(r_1\right)+\frac{(z+d)}{r_2}{q}_3\left(r_2\right)\]k_0{M}_y\}---\left(1\right)$

$E_y={\mathrm{\τ}}_E\{\[\frac{(y-y^{\′})(z-d)}{r_1^2}{q}_1\left(r_1\right)+\frac{(y-y^{\′})(z+d)}{r_2^2}{q}_1\left(r_2\right)\]P_z+\[-\frac{(z-d)}{r_1}{q}_3\left(r_1\right)-\frac{(z+d)}{r_2}{q}_3\left(r_2\right)\]k_0{M}_x\}---\left(2\right)$

$\begin{array}{c}{H}_x={\mathrm{\τ}}_H\{-\[\frac{(y-y^{\′})}{r_1}{q}_3\left(r_1\right)+\frac{(y-y^{\′})}{r_2}{q}_3\left(r_2\right)\]P_z+\[-\frac{{(y-y^{\′})}^2+{(z-d)}^2}{r_1^2}{q}_1\left(r_1\right)+q_2\left(r_1\right)\\ -\frac{{(y-y^{\′})}^2+{(z+d)}^2}{r_2^2}{q}_1\left(r_2\right)+q_2\left(r_2\right)\]k_0{M}_x+\[\frac{(x-x^{\′})(y-y^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(x-x^{\′})(y-y^{\′})}{r_2^2}{q}_1\left(r_2\right)\]k_0{M}_y\}\end{array}---\left(3\right)$

$\begin{array}{c}{H}_y={\mathrm{\τ}}_H\{\[\frac{(x-x^{\′})}{r_1}{q}_3\left(r_1\right)+\frac{(x-x^{\′})}{r_2}{q}_3\left(r_2\right)\]P_z+\[\frac{(x-x^{\′})(y-y^{\′})}{r_1^2}{q}_1\left(r_1\right)+\frac{(x-x^{\′})(y-y^{\′})}{r_2^2}{q}_1\left(r_2\right)\]k_0{M}_x\\ +\[-\frac{{(x-x^{\′})}^2+{(z-d)}^2}{r_1^2}{q}_1\left(r_1\right)+q_2\left(r_1\right)-\frac{{(x-x^{\′})}^2+{(z+d)}^2}{r_2^2}{q}_1\left(r_2\right)+q_2\left(r_2\right)\]k_0{M}_y\end{array}---\left(4\right)$

其中，j是虚数单位，k₀和η₀分别表示自由空间中的波数和波阻抗，f为频率，c为真空中的光速，η₀＝120π，Pz表示垂直电偶极矩，它是一个复数，单位为A·m，Mx和My表示水平磁偶极矩，单位为A·m²，r₁表示近场采样数据点与偶极子阵列中每个偶极子的距离，r₂表示近场采样数据点与偶极子阵列中每个偶极子的镜像的距离；

r₁和r₂由下式表示：

r₁＝[(x-x′)²+(y-y′)²+(z-d)²]^1/2 (5)

r₂＝[(x-x′)²+(y-y′)²+(z+d)²]^1/2 (6)

q₁(r)、q₂(r)、q₃(r)表示为:q₁(r)、q₂(r)、q₃(r)没有物理意义，这是为了简化公式(1)-(4)的，否则公式显得太复杂；

$q_1\left(r\right)=\[\frac{3}{{\left(k_{0}r\right)}^2}+j\frac{3}{k_{0}r}-1\]f\left(r\right)---\left(7\right)$

$q_2\left(r\right)=\[\frac{2}{{\left(k_{0}r\right)}^2}+j\frac{2}{k_{0}r}\]f\left(r\right)---\left(8\right)$

$q_3\left(r\right)=\[\frac{1}{k_{0}r}+j\]f\left(r\right)---\left(9\right)$

其中，r为空间点的矢径的模值，e为自然常数；

因此用来等效PCB电路的偶极子的偶极矩(Pz,Mx,My)和近场采样数据场强(Ex,Ey,Hx,Hy)存在下述的映射关系：

$F=\left(\begin{array}{c}{\[Ex\]}_{M^2\×1}\\ {\[Ey\]}_{M^2\×1}\\ {\[Hx\]}_{M^2\×1}\\ {\[Hy\]}_{M^2\×1}\end{array}\right)=TX=T\left(\begin{array}{c}{\[Pz\]}_{N^2\×1}\\ {\[Mx\]}_{N^2\×1}\\ {\[My\]}_{N^2\×1}\end{array}\right)---\left(10\right)$

其中[Ex]，[Ey]，[Hx]，[Hy]分别表示x方向和y方向上的电场强度和磁场强度矩阵，近场采样点的个数为M×M；是待求的偶极矩矩阵，[Pz]是偶极子阵列的垂直电偶极矩矩阵，[Mx]、[My]是偶极子阵列的水平磁偶极矩，每种偶极子个数为N×N；映射矩阵描述的是空间中任意一点电磁场与偶极子矩阵之间关系，映射矩阵由下式表示：

$T=\left(\begin{array}{ccc}{T}_{ExPz}& T_{ExMx}& T_{ExMy}\\ T_{EyPz}& T_{EyMx}& T_{EyMy}\\ T_{HxPz}& T_{HxMx}& T_{HxMy}\\ T_{HyPz}& T_{HyMx}& T_{HyMy}\end{array}\right)---\left(11\right)$

每一个子矩阵的维数为M²×N²，表示的是由某一种偶极矩产生的某一方向上的场强；以T_ExMy为例，它表示所有y方向上的磁偶极矩My在x方向上产生的电场强度Ex；T矩阵总共有4M²×3N²个元素，每一个子矩阵由下式表示：

$T_{EyPz}(a,b)={\mathrm{\τ}}_E\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0-d)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0+d)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]---\left(12\right)$

T_ExMx(a,b)＝0 (13)

$T_{ExMy}(a,b)={\mathrm{\τ}}_E\[\frac{(z_0-d)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{(z_0+d)}{r_{+}}{q}_3\left(r_{+}\right)\]k_0---\left(14\right)$

$T_{EyPz}(a,b)={\mathrm{\τ}}_E\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0-d)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)(z_0+d)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]---\left(15\right)$

$T_{EyMx}(a,b)={\mathrm{\τ}}_E\[-\frac{(z_0-d)}{r_{-}}{q}_3\left(r_{-}\right)-\frac{(z_0+d)}{r_{+}}{q}_3\left(r_{+}\right)\]k_0---\left(16\right)$

T_EyMy(a,b)＝0 (17)

$T_{HxPz}(a,b)=-{\mathrm{\τ}}_H\[\frac{\left(y\right(a)-y^{\′}(b\left)\right)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{\left(y\right(a)-y^{\′}(b\left)\right)}{r_{+}}{q}_3\left(r_{+}\right)\]---\left(18\right)$

$\begin{array}{c}{T}_{HxMx}(a,b)={\mathrm{\τ}}_H\[-\frac{{\left(y\right(a)-y^{\′}(b\left)\right)}^2+{(z_0-d)}^2}{{r_{-}}^2}{q}_1\left(r_{-}\right)+q_2\left(r_{-}\right)\\ -\frac{{\left(y\right(a)-y^{\′}(b\left)\right)}^2+{(z_0+d)}^2}{{r_{+}}^2}{q}_1\left(r_{+}\right)+q_2\left(r_{+}\right)\]k_0\end{array}---\left(19\right)$

$T_{HxMy}(a,b)={\mathrm{\τ}}_H\left[\begin{array}{c}\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\\ \frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{+}}^2}{q}_1\left(r_{+}\right)+\end{array}\right]k_0---\left(20\right)$

$T_{HyPz}(a,b)={\mathrm{\τ}}_H\[\frac{\left(x\right(a)-x^{\′}(b\left)\right)}{r_{-}}{q}_3\left(r_{-}\right)+\frac{\left(x\right(a)-x^{\′}(b\left)\right)}{r_{+}}{q}_3\left(r_{+}\right)\]---\left(21\right)$

$T_{HyMx}(a,b)={\mathrm{\τ}}_H\[\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{-}}^2}{q}_1\left(r_{-}\right)+\frac{\left(x\right(a)-x^{\′}(b\left)\right)\left(y\right(a)-y^{\′}(b\left)\right)}{{r_{+}}^2}{q}_1\left(r_{+}\right)\]k_0---\left(22\right)$

$\begin{array}{c}{T}_{HyMy}(a,b)={\mathrm{\τ}}_H\[-\frac{(z_0-d)+{\left(x\right(a)-x^{\′}(b\left)\right)}^2}{{r_{-}}^2}{q}_1\left(r_{-}\right)+q_2\left(r_{-}\right)\\ -\frac{{(z_0+d)}^2+{\left(x\right(a)-x^{\′}(b\left)\right)}^2}{{r_{+}}^2}{q}_1\left(r_{+}\right)+q_2\left(r_{+}\right)\]k_0\end{array}---\left(23\right)$

上述公式中(x(a),y(a),z₀)为控件中任意一点的坐标，(x′(b),y′(b),d)表示偶极子的坐标，其中a和b分别表示空间中任意一点和偶极子的序号，并且满足以下约束条件a＝1,2,3...M²，b＝1,2,3...N²；r_-表示空间中任意一点与偶极子的距离，r₊表示任意控件一点与偶极子镜像的距离；

根据公式(5)和(6)得到采样点与偶极子之间响应的位置关系，将每一个采样点的坐标(x,y,z)用(x(a),y(a),z₀)表示，每一个偶极子的坐标(x',y',d)用(x'(b),y'(b),d)表示，因此得到近场采样数据点与偶极子的距离为：

r₁(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀-d)²]^1/2 (24)

近场采样数据点与偶极子的镜像的距离为：

r₂(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀+d)²]^1/2 (25)

将采样点和偶极子及其镜像之间相应位置关系r₁(a,b)和r₂(a,b)分别用r_-和r₊表示，代入公式(12)-(23)计算得到采样点与偶极子之间的映射矩阵T；

采用正则化方法来求解偶极子的偶极矩；

argmin{||F-TX||²+α||X||²} (26)

X(α)＝[T′T+α²I]^-1T′F (27)

式中“argmin”表示使得(26)取得最小值，α被称作正则化系数；采用通过对数坐标尺度描述||X_α||和||F-TX_α||的L曲线准则，α对应L曲线上曲率最大的点；将求解得到的正则化参数α代入(27)求得与干扰不敏感的偶极矩矩阵X(α)；|| ||表示一个矩阵的范数，F表示近场采样数据场强矩阵；

在任意观察平面上的设置M×M观察点，x方向上等间隔取M个观察点，y方向上等间隔取M个观察点，每一个观察点的坐标为(x”,y”,z”)，z”是任意观察平面与参考地平面之间的高度距离，且满足z”＞d＞0；将每一个观察的坐标(x”,y”,z”)用(x(a),y(a),z₀)表示，每一个偶极子的坐标(x',y',d)用(x'(b),y'(b),d)表示，因此得到观察点与偶极子的距离为：

r₁'(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀-d)²]^1/2 (28)

观察点与偶极子的镜像的距离为：

r₂'(a,b)＝[(x(a)-x′(b))²+(y(a)-y′(b))²+(z₀+d)²]^1/2 (29)

将观察点和偶极子及其镜像之间相应位置关系r₁'(a,b)和r₂'(a,b)分别用r_-和r₊表示，代入公式(12)-(23)计算得到采样点与偶极子之间的映射矩阵T_v；

通过等效偶极矩矩阵X(α)和映射矩阵T_v根据公式(30)求解任意观察平面的电磁场切向分量；

$T_{v}X\left(\mathrm{\α}\right)=F_v=\left[\begin{array}{c}{\mathrm{Ex}}_v\\ {\mathrm{Ey}}_v\\ {\mathrm{Hx}}_v\\ {\mathrm{Hy}}_v\end{array}\right]---\left(30\right)$

其中F_v是在观察平面上电磁场强度矩阵，该矩阵有4M²×1个元素；其中Ex_v，Ey_v，Hx_v，Hy_v分别表示x方向和y方向上的电场强度和磁场强度，每个矩阵有M²×1个元素；

根据简单无耗媒质中的麦克斯韦方程组得到

$\▿\×E\left(r\right)=-j\mathrm{\ω}\mathrm{\μ}H\left(r\right)---\left(31\right)$

$\▿\×\tilde{H}\left(r\right)=j\mathrm{\ω}\mathrm{\ϵ}\tilde{E}\left(r\right)---\left(32\right)$

式中E(r),H(r)分别表示电场和磁场的复数形式；表示E(r)的旋度，表示H(r)的旋度，ω为角频率，μ磁导率，ε为介电常数；把(31)、(32)展开，由z方向上的分量相等得到

其中Hz是观察平面上法向磁场分量，Ez是观察平面上法向电场分量；利用中心差分商代替偏导数，得到(33)、(34)的离散形式：

$Hz(m,n)=\frac{j}{2\mathrm{\ω}\mathrm{\μ}}\{\frac{\[{\mathrm{Ey}}_v(m+1,n)-{\mathrm{Ey}}_v(m-1,n)\]}{\mathrm{\Δ}x}-\frac{\[{\mathrm{Ex}}_v(m,n+1)-{\mathrm{Ex}}_v(m,n-1)\]}{\mathrm{\Δ}y}\}---\left(35\right)$

$Ez(m,n)=-\frac{j}{2\mathrm{\ω}\mathrm{\ϵ}}\{\frac{\[{\mathrm{Hy}}_v(m+1,m)-{\mathrm{Hy}}_v(m-1,n)\]}{\mathrm{\Δ}x}-\frac{\[{\mathrm{Hx}}_v(m,n+1)-{\mathrm{Hx}}_v(m,n-1)\]}{\mathrm{\Δ}y}\}---\left(36\right)$

其中m,n表示每一个观察点在x方向和y方向上的序号，且满足2≤m,n≤M-1，Δx、Δy分别表示x方向和y方向上近场采样点的间距，实现了求解任意观察平面上的电磁场的分布；

其特征在于，包括有下列步骤：

步骤一：利用电磁场探头测试近场电磁场强度数据；

步骤二：设置偶极子阵列；

步骤三：计算采样点和偶极子阵列之间的映射矩阵T；

步骤四：计算偶极矩矩阵X(α)；

步骤五：计算偶极矩矩阵X(α)在高于PCB电路任意高度观察平面上的切向电磁场分量；

步骤六：把步骤五求得的切向电磁场分量代入公式(32)和(33)进行计算，求解步骤五中平面上的法向电磁场分量。

2.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤一中，对PCB电路进行辐射强度测量；首先确定PCB电路的相对位置坐标，在PCB电路上方高度为d的采样层平面，均匀测量M×M个点的场强度，沿x方向等间距取M个采样点，采样间隔为ΔX，沿y方向等间距取M采样点，采样间隔为ΔY，其中ΔX＝ΔY；记录每一个采样点坐标为(x，y，z)，将采集到的x方向的电场强度记为E_x，y方向的电场强度记为E_y，x方向的磁场强度记为H_x，y方向的磁场强度记为H_y，每一个矩阵有M²×1个元素，将矩阵E_x、E_y、H_x和H_y组合成一个矩阵F，其中

该矩阵有4M²×1个元素。

3.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤二中，在PCB上方高度为d的平面上均匀放置N×N个偶极子，沿x方向等间距放置N个偶极子，沿y方向等间距放置N个偶极子，记录每一个偶极子坐标(x’,y’,d)，将步骤一中的采样点的坐标(x，y，z)和步骤二中的偶极子坐标(x’,y’,d)代入公式(24)和(25)中计算得到采样点和偶极子及其镜像的相对位置关系。

4.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤三中，通过步骤二求解采样点和偶极子及其镜像的相对位置关系，利用公式(12)-(25)求解出采样点和偶极子之间的映射矩阵T，该矩阵具有4M²×3N²个元素。

5.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤四中，利用步骤一和步骤三求解得到的电磁场强度矩阵F和映射矩阵T，通过公式(10)可知X(α)为待求的偶极矩矩阵，该矩阵有3N²×1个元素；利用正则化方法中的L曲线准则计算正则化系数α，将α代入公式(27)求解偶极矩矩阵X(α)。

6.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤五中，任意观察平面上的设置M×M观察点，x方向上等间隔取M个观察点，y方向上等间隔取M个观察点，每一个观察点的坐标为(x”,y”,z”)，z”是任意观察平面与参考地平面之间的高度距离，且满足z”＞d＞0；利用公式(12)—(23)，(28)和(29)求解得到观察点和偶极子之间的映射矩阵T_v；利用MATLAB软件计算公式(30)得到观察平面的电磁场切向分量。

7.根据权利要求1所述的采用偶极矩模型对PCB电路电磁场进行反演的方法，其特征在于：在步骤六中，在任意观察平面上利用步骤五中求解得到切向分量的电场磁场强度Ex_v,Ey_v,Hx_v,Hy_v代入公式(35)和(36)进行计算，求解观察平面上的法向电磁场分量。