硅集成电路衬底多频率点下综合耦合参数的快速提取方法

摘要

衬底多频点综合耦合参数提取方法属于IC－CAD技术领域，其特征在于，根据衬底耦合电容提取和综合参数提取的相似性，提出了任意频率下综合耦合参数的两步提取方式：第一步，先提耦合电容参数，第二步，将电容精确修正为该频率下的综合耦合参数。本发明中也同时给出了修正过程的高效实现方式。当使用本方法计算多频率点下耦合参数时，电容参数提取操作只需进行一次，从电容到综合参数的修正过程在各频点各进行一次，但本发明中的修正过程可高效完成，每次仅需很少计算即可完成。因此，本发明对多频点下的综合耦合参数提取问题具有较高效率，同时由于基于严格的数学推导，计算精度也得到了保持。

说明书

技术领域

集成电路计算机辅助设计(IC-CAD)中的三维衬底综合耦合参数提取软件设计领域。

背景技术

半导体集成电路(Integrated Circuit，IC)是当前电子工业、乃至信息产业的基石。随着电路制造技术的发展，电路的集成度不断增加，当今的很多芯片含有几千万乃至上亿个CMOS等器件。与此同时，RF(射频)电路、数模混合电路和SoC(System on a Chip，系统芯片)的广泛应用(其中一例就是手机终端的日益普及)，使得对噪声很敏感的模拟电路和高速开关操作(switch)的数字电路被制造在一块芯片上，共享同一个衬底(substrate又有人译作基片、基底等)。这种高集成度具有很多优势，比如能量消耗较少、成本较低等，但也带来了一些挑战。其中一个就是，由于共享同一块衬底，衬底将会传播一些耦合噪声，这就是所谓的“衬底耦合噪声”。一方面，数字电路的快速开关操作(switch)以及其他有源连接线会向衬底注入噪声电流。由于衬底材料(多数为硅，Silicon)并不是理想的绝缘体，会将一部分噪声传递给敏感的模拟电路，对后者的正常功能造成影响，甚至使其工作失败。同时，通过电源线、器件等和衬底之间的寄生耦合电容，来自电源线、器件等的电压变动将会引起衬底的电压浮动，这种电压浮动又会影响衬底的阈值电压的大小，进而影响衬底上其他器件的电路性能。在较高频率下，衬底的等效电阻可能带来能量损耗等，这将对电路产生重要影响。比如，RF电路中电感器件的品质因子(qualify factor)就很受衬底的电阻损耗影响。综上所述，衬底耦合的严重性使得在电路设计中各阶段都必须关注该问题。图1显示了一个示意图，其中虚线框内为衬底和其他电路的接口contact，又可称作接触、端口，下文中简写为端口。在后续的提取参数的过程中，只需要提取这些端口间的耦合参数即可。

当前，集成电路的设计流程如图2所示。首先要提出功能描述，然后经过逻辑设计、版图设计得到描述半导体工艺尺寸、结构的版图。这时需要进行“版图验证”，即通过计算机软件模拟等来验证上述设计是否能达到当初设定的要求。如果满足要求，就可进行下一步的生产制造等；否则要返回逻辑设计等进行必要的修正。重复这个迭代过程，直到版图验证表明设计确实能够满足要求为止。在版图验证中，一个重要的环节称为“寄生参数提取”，其中就包括计算衬底耦合参数等。随着电路规模的不断增大，和特征尺寸的不断缩小，以往的简单估计甚至直接忽略这些参数的做法已经很难达到足够的精度要求。只有精确提取(计算)耦合参数，才能进行正确的电路模拟与验证。

随着工艺的进步和时钟频率的提高，衬底耦合参数的提取得到了越来越多的来自学术界和工业界的关注。当前，为进行精确的提取，需对电路版图的三维模型进行数值模拟，其主要方法包括区域型解法(有限差分方法和有限元方法)、基于Green函数的方法和直接边界元方法(Boundary Element Method，BEM)等。

与区域型解法相比，直接边界元法的优点是精度高，较少的离散变量和较强处理复杂边界形状的能力。与基于Green函数的方法相比，直接边界元方法具有适应复杂结构等优点。虽然直接边界元方法有上述优点，但由于工业实际需求越来越多、三维衬底结构也越来越大，如何提高基于直接边界元方法的三维提取方法的计算效率就成了当务之急。

在工作频率比较低时，衬底耦合参数表现为纯阻性，即等效电路中端口(contact，即与其他电路器件等的接口区域，如图1虚线框所示)之间只有寄生电阻。其提取方法包含很多，其中本人在文献“Substrate Resistance Extraction with Direct Boundary Element Method，”Asiaand South Pacific Design Automation Conference(2005年亚洲及南太平洋地区设计自动化会议)，pp.208-211，Jan.2005(以下简写作ASP-DAC2005)中将直接边界元方法应用到衬底耦合电阻提取中，取得了良好的效果。

在较高频率下，衬底耦合参数(其中包含电阻、电容)随着频率改变，即等效电路中端口之间既有电阻连接，也有电容连接，且它们的数值一般都随频率变化。此类参数下文中称综合耦合参数。本人在文献“A New Boundary Element Method for Accurate Modeling of LossySubstrates with Arbitrary Doping Profiles，”Asia and South Pacific Design Automation Conference，Yokohama，Japan，pp.683-688，Jan.2006(以下简写作ASP-DAC2006)将直接边界元方法应用到综合耦合参数的提取中，得到了较好的效果。由于该参数在不同频率下具有不同的值，往往需要计算其在多个频率点的具体数值，以此来全面评估耦合性能。上述第二篇文献虽然对此做了一些改进，但仍需在每个频率点建立、并求解一个复数的线性方程组。这在变量数很大时，效率将变得很低，因为方程组求解过程往往很慢，同时也需要较多计算机内存。更糟的是，当需要计算很多个频率点下的耦合参数时，总体计算时间(等于每频点计算时间的累加和)会很高，效率很低。

发明内容

本发明给出一种有效提取综合耦合参数(包含电阻R、电容C两部分)的方法，该方法对多频率点下耦合参数的计算问题效率较高。

为便于理解，图3显示了一个典型的衬底结构的示意图。本方法称作衬底综合耦合参数提取器，简写为SubRCExtract(Substrate coupling Resistance and Capacitance Extractor)。主要创新点在于，不同于上文文献ASP-DAC2006使用直接边界元方法在每频率点下依次计算综合参数(其流程图见图4)，而是把综合耦合参数的提取其分为两个步骤：

1.ExtractCap：提取衬底端口之间耦合电容；

2.ReviseCap2RC：将该电容修正为综合参数。

其中，ExtractCap和第一篇文献中所述的电阻提取方法非常类似(详见下文)，可能运行比较慢，但只需计算一次。ReviseCap2RC需在每个频率点都进行一次，但每次费时很少。总体上讲，SubRCExtract对多频率点参数计算问题将会有更好的效率。同时，因为是基于严格的数学推导，和直接使用直接边界元提取方法(如上述文献ASP-DAC2006)相比，其中不会有精度的损失。

SubRCExtract方法的输入为描述端口和衬底介质等几何结构的信息，输出为端口之间的综合耦合参数矩阵。方法主要包括两个模块：

【模块0】读取输入文件，得到形体信息(含衬底长、宽、高、电导率、介电常数，端口的位置、大小，工作频率等)

【模块1】进行衬底耦合电容的提取，即ExtractCap。

【模块2】将耦合电容进行数学修正，得到综合参数，即ReviseCap2RC。

下面分节详细介绍各模块。为更好理解模块2中的修正方法，在模块1之后，加入使用直接边界元进行综合参数提取的原理。

1使用直接边界元提取耦合电容(ExtractCap)

以下原理和文献ASP-DAC2005中的边界元提取耦合电阻很类似，在本小节最后会给出两者的异同之处。

以下讨论基于模型衬底图3展开，并令其内所有的介质都被认为是绝缘体(即电导率σ为0)，此时，端口之间的耦合作用将是纯电容性的。其耦合电容可通过如下的直接边界元方法求解得到。

任一介质M_i内，电压u满足Lapalce方程：

²u＝0，介质M_i中    (1)

上式的混合边界条件是：

u＝u，在端口表面上(强加边界条件)    (2a)

q＝0，在其他表面上(自然边界条件)    (2b)

其中u是在端口上预设的电压(1V或者0V)。q是法向电场强度。除此之外，在介质a和b的交界面上，电压和电压移连续：

u_a＝u_b，                       (3a)

ε_aq_a＝ε_bq_b，                 (3b)

其中ε_a和ε_b分别是a和b的介电常数，q是法向电场。

求解方程式(1)，可直接得到电场分布。然后根据电场分布数值计算端口表面的电荷，其数值和耦合电容的数值紧密相关，如下文式(7)所示。下面是具体求解电场分布的方法。

使用Green公式，并且选用自由空间Green函数作为权重函数，可以将式(1)转换成在介质M_i的边界Г_i上成立的边界积分方程：

其中u_s表示某配置点s的电压(配置点的意义将在稍后解释)。将边界Г_i划分为N_i边界元Г_ij，j＝1，2，...，N_i。式(4)的离散形式为：

其中，每个边界元Г_ij上的电压u和法向电场q都是常数，所以可挪到积分号之外。其中配置点s(或称为采样点，取采样点是为了列出含有未知量的方程组，以求解出未知量)取为每一边界元的中心点。u^*是Laplace方程的基本解，等于自由空间Green函数，形式为其中r(s，v)表示配置点s和变量点v间距离的大小，变量点v也取为Г_i上每一边界元的中心点，可和点s重合。q^*是u^*沿边界外法向n的导数，表达式为 $>>->>1>>4>\pi >>r>2>>>(>s>,>v>)>>>>>>\partial ver>>>r>>(>s>,>v>)>>>‾>>>>\partial ver>>n>‾>>>>,>>>其中r(s，v)表示距离矢量。每个变量点上(即对应的边界元上)都有两个物理量，即电压u和电场q。根据边界条件(2)，如果该边界元处在端口(contact)表面，其上u已知，即预设的电压数值，而电场q未知；否则，即该边界元在非端口区，其上电场q为0，电压u未知。所以，每个边界元实际上只有一个变量，或者为u类型变量或者为q类型变量。由(5)可知，它们的系数形式分别为和$

使用交界面方程式(3a)和式(3b)，通过变量变换(具体方式见下文)，将所有介质的边界积分方程式(5)连接起来，得到统一的方程组，写作矩阵形式为：

                A_capx_cap＝b_cap，    (6’)

其中，一个配置点s将产生A_cap的一行元素，即在该配置点时由式(5)产生的u或者q的系数，显然，这些系数可表示为或乘以1或实数ε_b/ε_a，其中ε_a和ε_b分别是相互交界的介质a和介质b的介电常数。b_cap是将已知电压u(乘以系数)挪到方程右端而形成的，具体形式为其中Г_m是端口m的表面。x_cap包含各介质中电压u型和电场q型变量。当需要求解很多对端口间的耦合电容时，方程组(6’)右端将含有多个向量，简写做：

           A_capX_cap＝B，对电容提取                    (6)这里B是由多个向量组成的矩阵。待解矩阵X_cap包含不同偏压设置时的电压u和电场q数值。

求解方程组(6)，可直接得到q的数值。如果预先设定编号为m的端口电压为1V，其余端口k(k≠m)设为0V，端口m、k之间的耦合电容是

$>>>C>mk>>=>>\int >>\Gamma >k>>>ϵqd\Gamma >,>->->->>(>7>)>>>>$

其中Г_k是端口k的表面，q是Г_k上的电场，ε是紧贴着端口k的介质的介电常数。由此可见，由电场q分布可很容易计算耦合电容¹。

注意，在方程式(3)中有两对变量u_a、q_，a和u_b、q_，b，可只保留其中一对，另外一对通过变量代换得到。比如，如果保留u_b，q_，b(属于介质b)，并代入式(3)式得到u_a和q_a：u_a＝u_b，q_a＝ε_b/ε_aq_b。对介质b中的配置点而言，变量q_b的系数的形式为其中Г_bj是介质b上的某边界元；而对介质a中的配置点而言，变量q_b的系数形式将变为其中Г_aj是介质a上的某边界元。

2使用直接边界元提取综合参数

以下原理来自文献ASP-DAC2006。

实际的衬底都是有损的，构成衬底的介质的电导率并不为零。根据文献S.Ramo，J.R.Whinnery，and T.van Duzer，Fields and Waves in Communication Electronics，3^rd ed.New York：Wiley，pp.283-287，pp.572-676，1994，这样的“有损介质”可以在一定频率下可被近似看作“绝缘体”，拥有复数介电常数ε_c＝ε-jσ/ω，然后该“有损介质”就可以当作这样的“绝缘体”来处理。为方便，这里令“绝缘体”a和b的介电常数分别为ε_c，a＝ε_a-jσ_a/ω，ε_c，b＝ε_b-jσ_b/ω，其中σ_a、σ_b分别是介质a和b的电导率，ε_a和ε_b是相应的实数介电常数，ω是角频率，等于2πf，f为频率。

为计算综合参数，我们在端口上设置1V或者0V的电压。边界条件式(2a)和式(2b)保持不变。Laplace方程式(1)仍然成立，在介质a和b交界面上，电压仍是连续的，但式(3b)变为

                    ε_c，aq_a＝ε_c，bq_b.(3b’)

将式(3b’)做如下变形：

¹这里还有一个细节，根据文献ASP-DAC2005，如果分别将式(3b)中的物理量εa和εb换成σa和σb(分别为介质a和b的电导率)，其他公式保持不变，生成并求解方程组之后，可得到耦合电阻R_mk，其中 $>>>1>>R>mk>>>=>>\int >>\Gamma >k>>>\sigma >>E>n>>d\Gamma >,>>>形式上很像电容求解式(7)。$

        (σ_a+jωε_a)q_a＝(σ_b+jωε_b)q_b.           (3_freqRC)

式(3_freqRC)可以解释为全电流守恒(包含传导电流和位移电流)。使用式(3_freqRC)代替式(3b)，保持其他不变。采用类似的变量代换方式，可得到一个方程组：

             AX＝B，对综合参数提取                     (8)

其中系数矩阵A是复数矩阵。和上文A_cap类似地，A中元素的数值(即u或者q的系数)为或乘以1或复数(σ_b+jωε_b)/(σ_a+jωε_a)，其中σ_a、σ_b和ε_a、ε_b分别是相互交界的介质a和介质b的电导率和介电常数。只要所预设的偏压是相等的，式(8)中的B和式(6)中的B形成原理一致，数值也均为其中Г_m是端口m的表面。

求解式(8)，流经某个端口的电流将为这样，如果编号m和k的端口上预设电压分别为1V和0V，它们之间的综合耦合参数就是：

$>>>Z>mk>>=>>1>>>\int >>\Gamma >k>>>>(>\sigma >+>j\omega ϵ>)>>qd\Gamma >>>->->->>(>9>)>>>>$

其中Г_k是端口k的表面，q是Г_k上的电场，σ和ε分别是紧贴着端口k的介质的电导率和介电常数。

在文献ASP-DAC2006中，需要在每个频率点下建立并求解方程组(8)，得到该频率下的综合耦合参数。但是，当需要许多频率点下的综合耦合参数时，列方程、解方程的过程将重复很多次，计算代价较高。更糟糕的是，方程组(8)往往很大，而且是复数方程组，比实数方程组更难求解。

为了保证比较好的计算效率，可通过如下“修正”方式来计算多频率点下的耦合参数。

3从电容到综合参数的修正方法ReviseCap2RC

4.2节所述的综合参数和4.1节的电容提取有很多的相似性。这种相似性可被充分利用，得到一种新的综合耦合参数求解方法。首先，可按照4.1节的原理计算出耦合电容。下一步，把这些电容修正为某频点下的综合耦合参数。第一步需要建立并求解方程组(6)。它是实数方程组，其求解得到的信息和频率无关，可被存储起来，在需要时可随时读入。第二步在每个频率点均需实施一次，但是，它所需要的计算很少。在需提取多个频率点下的耦合参数时，新方法将有较大优势，因为费时的只是第一步，第二步仅需很少时间。简易流程图见图5，在后文，将给出更易操作的流程图11。

下面，我们将给出详细的从电容到综合参数的修正方法。

定理1：电容提取与综合参数提取系数矩阵相似性定理

对同一衬底几何结构和同样的边界元划分，式(8)中矩阵A和式(6)中矩阵A_cap有如下关系：A＝A_cap+UV^T，矩阵U和V^T的构造见下文。

方程组(8)中的矩阵A和方程组(6)中的矩阵A_cap有着极大相似性。首先来看矩阵A是怎样形成的。对每个衬底介质都列出式(5)，并使用交界面上变量的代换来连接成为统一的方程组。这里以M₁、M₂之间交界面Г₁₂上的u型变量(未知电压)和q型变量(未知法向电场)为例。在该交界面两侧的u、q变量满足交界面方程(3a)和(3_freqRC)，做如4.1节所述的变量代换。这里假设保留属于介质M₂的变量，分别命名为u和q₁₂，则属于介质M₁的变量为1×u和f₁₂×q₁₂(根据交界面方程式(3a)和式(3_freqRC))，其中 $>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>.>>>在此假设下，对M$₂的某配置点s₂，变量u的系数为积分值变量q₁₂的系数为其中Г_2e是介质M₂边界上的某边界元(上述两个积分值都只和几何信息有关系)；而对M₁的某配置点s₁，变量u的系数为其中Г_1e是介质M₁边界上的某边界元；变量q₁₂的系数则为其中 $>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>.>>>同样的，对其他介质交界面上的变量需做类似的代换。显然，只有交界面上某些变量的系数是比较复杂的其他的系数均可表示为或$

上述过程也可参见示意图6，重点考察交界面变量q₁₂、q₂₃的系数。S₂₁作为q₁₂对M₂中某配置点的系数，等于和M₂相关的积分由式(3_freqRC)，q₁₂对M₁中某配置点的系数为和M₁相关的积分再乘以 $>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>,>>>也就是f$₁₂S₁₂。S₃₂和f₂₃S₂₃，含义类似：S₃₂是和M₃相关的积分值； $>>>f>23>>=>>>>\sigma >3>>+>j\omega >>ϵ>3>>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>,>>>S$₂₃是和M₂相关的积分值。

为叙述方便，上文f₁₂、f₂₃称为交界面频变系数，定义式为： $>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>,>>f>23>>=>>>>\sigma >3>>+>j\omega >>ϵ>3>>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>.>>>$

同理，使用边界元进行电容提取时(详见4.1节原理)，也需进行变量代换。为简单，仅以M₁、M₂之间交界面Г₁₂上的电场变量为例。仍假设，保留属于介质M₂的电场变量，命名为q₁₂，则根据交界面方程式(3a)和式(3_b)，属于介质M₁的电场变量则为c₁₂×q₁₂，其中 $>>>c>12>>=>>>ϵ>2>>>ϵ>1>>>.>>>对介质M$₂的某配置点s₂，q₁₂的系数为其中Г_2e是介质M₂边界上的某边界元。则对M₁的某配置点s₁，q₁₂的系数为其中Г_1e是介质M₁边界上的某边界元， $>>>c>12>>=>>>ϵ>2>>>ϵ>1>>>.>>>对。该过程参见示意图7，其中$ >>>c>23>>=>>>ϵ>3>>>ϵ>2>>>.>>>$$

为叙述方便，上文c₁₂、c₂₃称为交界面电容系数，定义式为： $>>>c>12>>=>>>ϵ>2>>>ϵ>1>>>,>>c>23>>=>>>ϵ>3>>>ϵ>2>>>.>>>$

所以，如果我们在图6中直接令 $>>>f>12>>=>>c>12>>=>>>ϵ>2>>>ϵ>1>>>,>>f>23>>=>>c>23>>=>>>ϵ>3>>>ϵ>2>>>,>>>该矩阵就变成方程组(6)的系数矩阵、如图7所示的A$_cap。可以看出，两个矩阵的差别只是和交界面变量(如q₁₂、q₂₃)相关的元素。两个矩阵的差别见图8。该差别可表示为UV^T形式，其中只有V^T和频率相关。U和V可以是高度稀疏的。下面给出一种具体构造方式。

U的大多数元素均为0，仅交界面变量的系数可能非零，且非零值等于A_cap中含c₁₂、c₂₃的部分，详见图7和图8(繁琐地说，U的元素等于对交界面序号小的介质中配置点时的交界面电场变量的系数)。这样U就不需要额外计算，只需提取A_cap的相应元素即可。

V的形状很简单，基本为对角阵，但除了对应于交界面上电场变量的那些元素非零之外，其他元素都为0；非零元素计算也很简单，为交界面频变系数除以交界面电容系数再减掉1，比如(f₁₂-c₁₂)/c₁₂，(f₂₃-c₂₃)/c₂₃，数值上仅仅取决于交界面双方介质的物理性质(电导率σ、介电常数ε)和角频率ω。可见，V依赖于频率。不同频率下，U是固定的，V是变的。

另外，数学上有一个严格成立的Woodbury公式：

(A+UV^T)^-1＝A^-1-A^-1U(I+V^TA^-1U)^-1V^TA^-1，

其中A，U，V都是任意矩阵，成立的条件为A和I+V^TA^-1U都可逆。参见文献Weisstein EW.Concise Encyclopedia of Mathematics.Boca Raton，Fla.：CRC Press，1999。

根据以上定理和Woodbury公式，可得到另一定理：

定理2：修正电容参数可得综合参数定理

对同一衬底几何结构和同样的边界元划分，方程组式(8)的解A^-1B可通过式(6)的解A_cap^-1B来获得：

$>>X>=>>A>>->1>>>B>=>>(>>A>cap>>+>>UV>T>>>)>>->1>>>B>>>>$

$>>=sup>>A>cap>>->1>sup>>B>-sup>>A>cap>>->1>sup>>U>>>(>I>+>>V>T>sup>>A>cap>>->1>sup>>U>)>>>->1>>>>V>T>sup>>A>cap>>->1>sup>>B>>>$

$>>=>>X>cap>>-sup>>A>cap>>->1>sup>>U>>>(>I>+>>V>T>sup>>A>cap>>->1>sup>>U>)>>>->1>>>>V>T>>>X>cap>>->->->>(>10>)>>>>$

换句话说，综合参数可以通过修正电容参数来获得。由于X_cap和A_cap^-1B是和频率无关的，它们可只被计算一次，然后存储起来。计算任一频率下综合参数时，可重新调入X_cap和A_cap^-1B，并使用依赖于频率的(I+V^TA_cap^-1U)^-1来做必要的“修正”，进而得到该频率下综合参数的数值。

粗略看起来，式(10)虽然严格成立，但实用意义不大，因为计算(I+V^TA_cap^-1U)^-1和直接求解A^-1一样复杂、一样困难，原因是它们两个等大小，都是困难的求逆操作。但是，使用下面的定理，可高效率地计算(I+V^TA_cap^-1U)^-1：

定理3：(I+V^TA_cap^-1U)^-1的高效计算定理

求解(I+V^TA_cap^-1U)^-1等价于求一个很小矩阵M的逆。M是由(I+V^TA_cap^-1U)^-1中非零块组成，这些非零块不包含主对角线上的子单位阵。M的大小等于“不同”介质之间交界面上的边界元的总数。这里，“不同”的意思是交界面双方介质的电导率(σ)和介电常数(ε)不全相等。

具体做法如下。由于U和V^T的高度稀疏性，I+V^TA_cap^-1U、(I+V^TA_cap^-1U)^-1都和单位矩阵I很相像，参见图9。将I+V^TA_cap^-1U中非零块抽取出来(这些非零块不包含主对角线上等于单位阵的非零子块)，组成矩阵M，记下M各小块的划分方式。计算逆矩阵W＝M^-1，按照相同的方式把W划分为很多小块，然后把各块按照上文“抽取”的方式反填到和I+V^TA_cap^-1U等大小的单位阵中，这样就得到了(I+V^TA_cap^-1U)^-1。以上过程可简单表示为：

(11)的等效性很容易在理论上得到证明：只需证明图9左图和图10右图确实互为逆阵，也就是它们的乘积为单位阵；根据矩阵乘法法则，很容易验证乘积确实为单位阵。

明显地，M比I+V^TA_cap^-1U小很多很多。它的精确大小等于“不同”介质之间交界面上的变量总数，“不同”的意思是交界面双方介质的电导率(σ)和介电常数(ε)不全相等。原因举例说明如下。参看图8中矩阵V^T的组成。如果介质M₁和M₂有全部相同的σ和ε，则交界面频变系数 $>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>=>1>,>>>而电容系数$ >>>c>12>>=>>>ϵ>2>>>ϵ>1>>>=>1>.>>>所以，f$$₁₂-c₁₂必然等于0，体现在矩阵V中就是相应位置的元素为0。如果M₁和M₂的σ和ε不全相等，那么因子f₁₂-c₁₂是和频率相关的，也就是不等于0，在矩阵V相应位置处产生非零元。根据矩阵乘法规则，只有V(或者V^T)非零的元素才会在I+V^TA_cap^-1U产生(非主对角线上的)非零元。所以，M的精确大小等于“不同”介质之间交界面上的变量总数。

现在，我们可以看到，通过式(10)来求解综合参数的效率将会很高，这是因为在某频率下，只需求解小矩阵M的逆阵。另外，其他的矩阵相乘操作如A_cap^-1U、V^TX_cap等也可以因为矩阵U、V大多数元素为0而得到加速。式(10)基于数学上严格的Woodbury公式，所以是严格成立的。

式(10)内存使用量也不太多，因为只有A_cap^-1U、I+V^TA_cap^-1U和(I+V^TA_cap^-1U)^-1中的非零元，还有X_cap需要使用内存。其中的非零元所占比例非常小，所以，内存使用量不会太多。矩阵A_cap^-1U、I+V^TA_cap^-1U和(I+V^TA_cap^-1U)^-1的稀疏性可参看后文的实例。

4方法详细流程

下面给出实用化的方法流程。

1)输入几何形体信息，包括：衬底各介质M_i的长、宽、高、电导率σ、介电常数ε；端口区的位置、大小；工作频率。

2)计算衬底不同端口之间的耦合电容ExtractCap，生成A_cap，抽取A_cap中部分列的部分元素组成U；具体为：

步骤2.1.在主端口m上设1V电压，其余端口(包括k)设为0V；

步骤2.2.把介质M_i的边界划分为N_i边界元Г_ij，j＝1，2，...，N_i；

步骤2.3.分别以每个边界元中心点为配置点即采样点s，列出方程组A_capX_cap＝B；

其中，若设变量在边界元Г_ij上，矩阵A_cap的元素值为：

电压u型变量的系数；

或非交界面电场q型变量的系数；

或交界面Г_ab上q型变量，且配置点s在M_b上；

或交界面Г_ab上q型变量，且配置点s在M_a上；

其中，交界面Г_ab位于介质M_a和M_b之间，a＜b； $>>>u>*>>=>>1>>4>\pi r>>(>s>,>v>)>>>>,>>>r(s，v)为配置点s和变量点v间距离的大小，q*为u*在Г$_ij法向上n的导数，表达式为 $>>->>1>>4>\pi >>r>2>>>(>s>,>v>)>>>>>>\partial ver>>>r>>(>s>,>v>)>>>‾>>>>\partial ver>>n>‾>>>>,>>>其中r(s，v)表示距离矢量；$

步骤2.4.矩阵B的每一列对应着一个主端口m的设定；对某配置点s，B某列的某元素值为其中Г_m是主端口m的表面，q^*定义同上；

步骤2.5.求解方程组A_capX_cap＝B，得到 $>>>X>cap>>=sup>>A>cap>>->1>sup>>B>;>>>X$_cap包含电场q型变量的值；

步骤2.6.可由X_cap值计算耦合电容：

$>>>C>mk>>=>>\int >>\Gamma >k>>>>ϵ>>(>k>)>>>>q>>(>k>)>>>d\Gamma >>>$

其中Г_k是端口k的表面，q_(k)是Г_k上的电场，ε_(k)是紧贴着端口k的介质的介电常数；

步骤2.7.生成矩阵U：从A_cap中复制形式为的元素，其余元素设为0；

步骤2.8.和步骤2.4类似地求得A_cap^-1U；

3)对每个频率点，执行以下操作，将电容修正为该频率下的综合耦合参数：

步骤3.1.计算该频率点下的V矩阵：对那些对应于介质M_a和M_b之间交界面上的电场型变量的主对角元，设定它们的数值为(f_ab-c_ab)/c_ab，其中a＜b， $>>>f>ab>>=>>>>\sigma >b>>+>j\omega >>ϵ>b>>>>>\sigma >a>>+>j\omega >>ϵ>a>>>>,>>c>ab>>=>>>ϵ>b>>>ϵ>a>>>,>>>\sigma$_a、σ_b和ε_a、ε_b分别是介质a和b的电导率和介电常数；其余元素设为0；

步骤3.2.生成I+V^TA_cap^-1U，对其进行压缩得到较小矩阵M；

步骤3.3.对M求逆，得到W＝M^-1，并将W反填到单位阵，得到(I+V^TA_cap^-1U)^-1；

步骤3.4.实施矩阵乘，得[A_cap^-1U][(I+V^TA_cap^-1U)^-1][V^TX_cap]；

步骤3.5.矩阵相减，得到X_cap-{[A_cap^-1U][(I+V^TA_cap^-1U)^-1][V^TX_cap]}；

步骤3.6.由式(9)计算综合参数，并输出；

步骤3.7.所有频率点下耦合参数都已计算完毕？若否，则转向步骤3.1；否则，转向步骤4；

4)结束。

以上过程见图11。

附图说明

图1衬底耦合效应示意图。通过衬底传递的耦合噪声将影响电路的某些部分，并进而影响整体性能。

图2一个有限区域的三维互连结构示意图(需要计算其中主导体与其他导体的耦合电容)。

图3由多层介质、多个端口(和接地板)组成的衬底结构示意图。

图中标号意义：

Mi：介质i，i＝1～3

σ_i：电导率，非零

ε₁：介电常数。

Г₁₂：介质M1和介质M2的交界面

Г₂₃：介质M2和介质M3的交界面 

图4直接使用边界元方法提取综合参数的流程。

(参见文献“A New Boundary Element Method for Accurate Modeling of Lossy Substrates with ArbitraryDoping Profiles，”Asia and South Pacific Design Automation Conference，Yokohama，Japan，pp.683-688，Jan.2006)

图5本发明中较高频耦合参数计算方法的流程图。

注：灰色块中的矩阵元素非零。

图6对图3所示衬底，式(8)中矩阵A示意图。灰色块中的矩阵元素非零。

注：灰色块中的矩阵元素非零。

图7对图3所示示例，式(6)中矩阵A_cap示意图。

图8矩阵A和A_cap的差别。

注：

U的结构和A-Acap一致，其非零元数值等于Acap相应位置元素的数值。

V是非完全对角阵，由两个对角阵组成，分别为(f12-c12)/c12I1和(f23-c23)/c23I2，其中I1和I2均为单位阵。

图9左图为矩阵I+V^TA_cap^-1U的结构，右图为其压缩后的等价矩阵M。

图10把左图所示矩阵M的逆阵W＝M-1反填到一个单位阵中得到(I+V^TA_cap^-1U)^-1，如右图所示。

图11本发明中较高频耦合参数计算方法的详细流程图。

图12由三层介质组成的简单衬底算例。

图中符号及相关解释：

σ：介质的电导率，单位为1/(Ωm))；

ε_r：介质的相对介电常数，单位为1；

t：介质的高度，单位为um。

L：介质的长度，单位为um。

W：介质的宽度，单位为um。

顶层的两个平面端口，大小均为2um×2um，严格对称分布，如右俯视图所示。

具体实施方式

本发明的提取方法已用C++语言编程实现，可以在SUN SPARC系列工作站上的UNIX操作系统及PC机的Linux操作系统上运行。下面结合一个具体实例说明含SubRCExtract方法的执行过程。

图12显示了一个简单例子。为求简单，这里contact(端口)被看作平面，且放在顶部。

注意，这里为了简单说明执行过程，特意将例子的变量数降得很低。在以下的矩阵中，元素都小数点后只取一位数字。但在程序执行中，是以计算机内部的浮点数(float型或者double型)表示的。为求清晰，矩阵零元显示为空。

1)输入图示的衬底信息。

2)设置端口1电压为1V，端口2、接地板为0V。

步骤2.1.把各介质边界划分为边界元Г_ij。

步骤2.2.按照A_cap元素的积分公式，计算得其具体数值为：(矩阵A_cap为)

  6.3  -0.8  -1.4  -1.4  13.3  -1.4  -13.3  -0.8  6.3  -1.4  -1.4  13.3  -1.4  -13.3  -1.4  -1.4  6.3  -0.8  13.3  -1.4  -13.3  -1.4  -1.4  -0.8  6.3  13.3  -1.4  -13.3  -1.4  -1.4  -1.4  -1.4  35.3  -0.8  -9.3  -1.4  -1.4  -1.4  -1.4  9.3  6.3  -35.3  6.3  35.3  -1.4  -1.4  -1.4  -1.4  -0.8  -9.3  -1.4  13.3  6.3  -0.8  -1.4  -1.4  -1.4  -13.3  -1.4  13.3  -0.8  6.3  -1.4  -1.4  -1.4  -13.3  -1.4  13.3  -1.4  -1.4  6.3  -0.8  -1.4  -13.3  -1.4  13.3  -1.4  -1.4  -0.8  6.3  -1.4  -13.3  -0.8  9.3  -1.4  -1.4  -1.4  -1.4  6.3  -35.3  6.3  109.8  -1.5  -1.5  -1.5  -1.5  0.5  0.5  -1.5  45.6  6.0  -1.0  -1.7  -1.7  0.9  0.7  -1.5  45.6  -1.0  6.0  -1.7  -1.7  0.7  0.9  -1.5  45.6  -1.6  -1.6  5.9  -1.0  0.8  0.8  -1.5  45.6  -1.6  -1.6  -1.0  5.9  0.8  0.8  -0.8  27.7  -1.9  -1.0  -1.3  -1.3  7.1  1.0  -0.8  27.7  -1.0  -1.9  -1.3  -1.3  1.0  7.1

生成U大部分元素为0，其余非领元素为A_cap的一部分，即A_cap中以黑体显示的部分数字。

  -13.3  -13.3

  -13.3  -13.3  -9.3  -35.3  -9.3  -13.3  -13.3  -13.3  -13.3  -35.3

3)求解方程组A_capX_cap＝B，其中B＝

                                 1.4

                                 1.4

                                 1.4

                                 1.4

                                 -6.3

                                 0.8

                                 0

                                 0

                                 0

                                 0

                                 0

                                 0

                                 0

                          0

                          0

                          0

                          0

                          0

                          0

得 $>>>X>cap>>=sup>>A>cap>>->1>sup>>B>=>>>$

                          1.0

                          1.0

                          1.0

                          1.0

                          -0.02

                          0.9

                          -0.02

                          0.8

                          0.8

                          0.8

                          0.8

                          0.8

                          -0.02

                          0.6

                          0.6

                          0.6

                          0.6

                          0.5

                          0.5

同时也可生成 $>sup>>A>cap>>->1>sup>>U>=>>>$

 -3.3-2.8 -3.3-2.8 -3.3-2.8 -3.3-2.8 -0.9-0.8

  -6.5  -5.7  0.1  -0.8  -6.1  -8.5  -6.1  -8.5  -6.1  -8.5  -6.1  -8.5  -5.7  -11.4  0.1  0.2  -4.4  -8.9  -4.4  -8.9  -4.5  -9.0  -4.5  -9.0  -4.0  -8.0  -4.0  -8.0

4)对下一频率点(这里以计算频率为10¹⁰Hz即10GHz下耦合参数为例)：

步骤4.1生成矩阵V^T：对两个交界面计算以下系数：

$>>>f>12>>=>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>>\sigma >1>>+>j\omega >>ϵ>1>>>>=>>>2000>+>j>6.28>\times >>10>10>>\times >11.8>\times >8.85>\times >>10>>->12>>>>>10>+>j>6.28>\times >>10>10>>\times >11.8>\times >8.85>\times >>10>>->12>>>>>=>1.4>\times >>10>2>>->91.3>i>,>>>$

$>>>f>23>>=>>>>\sigma >3>>+>j\omega >>ϵ>3>>>>>\sigma >2>>+>j\omega >>ϵ>2>>>>=>>>1>\times >>10>>->12>>>+>j>6.28>\times >>10>10>>\times >3.9>\times >8.85>\times >>10>>->12>>>>>2000>+>j>6.28>\times >>10>10>>\times >11.8>\times >8.85>\times >>10>>->12>>>>>=>3.5>\times >>10>>->6>>>+>1.1>\times >>10>>->3>>>i>,>>>$

注：j和i相同，均为复数单位，也就是-1的一个平方根，下同。

$>>>c>12>>=>>>ϵ>2>>>ϵ>1>>>=>>>11.8>\times >8.85>\times >>10>>->12>>>>>11.8>\times >8.85>\times >>10>>->12>>>>>=>1>,>>c>23>>=>>>ϵ>3>>>ϵ>2>>>=>>>3.9>\times >8.85>\times >>10>>->12>>>>>11.8>\times >8.85>\times >>10>>->12>>>>>=>0.3>.>>>$

所以，根据V^T的计算公式，计算得到V^T：

  139.1-j91.

    3 -1+j3.3e- 3

步骤4.2生成I+V^TA_cap^-1U，对其进行压缩得到较小矩阵M；

$>>I>+>>V>T>sup>>A>cap>>->1>sup>>U>=>>>$

    1    1    1    1    1    1 16.5-i10.2 -108.1+i71.0    1    1    1    1    1 -0.1+i3.6e- 4 0.8+i7-3e-4 1

  1  1  1  1  1

其压缩形式矩阵M为：

16.5-i10.2              -108.1+i71.0

-0.1+i3.6e-4            0.8+i7.3e-4

步骤4.3对M求逆，得到W＝M^-1，并将W反填到单位阵，得到(I+V^TA_cap^-1U)^-1；

W＝

                    127.2-i

1.0-i5.9e-2         95.4

1.4e-1-i9.0e-3      19.4-i13.8

将W反填到足够大的单位阵I中，得到的即是(I+V^TA_cap^-1U)^-1

  1  1  1  1  1  1  1.0-i5.9e-2  127.2-i  95.4  1  1  1  1  1  1.4e-1-i9.0e-  3  19.4-i13.8  1  1

    1    1    1    1

步骤4.4实施矩阵乘，得 $>>[sup>>A>cap>>->1>sup>>U>]>[>>>(>I>+>>V>T>sup>>A>cap>>->1>sup>>U>)>>>->1>>>]>[>>V>T>>>X>cap>>]>=>>>$

(其中V^TX_cap只是将X_cap的第7行放大139.1-i91.3倍，将第13行-1+i3.3e-3倍)

             -3.3e-2+i3.0e-2

             -3.3e-2+i3.0e-2

             -3.3e-2+i3.0e-2

             -3.3e-2+i3.0e-2

             -8.8e-3+i8.1e-3

             -6.5e-2+i5.9e-2

             -1.5e-2+i6.1e-5

             -1.2e-1+i5.9e-2

             -1.2e-1+i5.9e-2

             -1.2e-1+i5.9e-2

             -1.2e-1+i5.9e-2

             -1.7e-1+i6.0e-2

             3.4e-3+i1.2e-3

             -1.4e-1+i4.6e-2

             -1.4e-1+i4.6e-2

             -1.4e-1+i4.7e-2

             -1.4e-1+i4.7e-2

             -1.2e-1+i4.2e-2

             -1.2e-1+i4.2e-2

步骤4.5矩阵相减，得到 $>>>X>cap>>->\{>[sup>>A>cap>>->1>sup>>U>]>[>>>(>I>+>>V>T>sup>>A>cap>>->1>sup>>U>)>>>->1>>>]>[>>V>T>>>X>cap>>]>\}>=>>>$

             1.0+i3.0e-2

             1.0+i3.0e-2

             1.0+i3.0e-2

             1.0+i3.0e-2

             -6.3e-3-i8.1e-3

            9.5e-1+i5.9e-2

            -5.3e-2+i6.1e-5

            9.5e-1+i5.9e-2

            9.5e-1+i5.9e-2

            9.5e-1+i5.9e-2

            9.5e-1+i5.9e-2

            9.5e-1+i6.0e-2

            -1.8e-1+i1.2e-3

            7.4e-1-4.6e-2

            7.4e-1-i4.6e-2

            7.5e-1-i4.7e-2

            7.5e-1-i4.7e-2

            6.7e-1-i4.2e-2

            6.7e-1-i4.2e-2

步骤4.6由式(9)计算综合参数，并输出：

         端口1                   端口2         GroundPlane(接地

                                               板)

端口1    -3.37e+03+j1.35e+05     -7.21e+04     1.02e+04

                                 -j6.22e+05    -j1.63e+05

4.5)所有频率点下耦合参数都已计算完毕？若否，则转向4.1)；若是，转5)。

5)结束。

使用相同的例子和相同的参数，文献ASP-DAC2006中方法的计算结果为：

         端口1                  端口2            GroundPlane(接地

                                                 板)

端口1    -3.37e+03+j1.35e+05    -7.21e+04        1.02e+04

                                -j6.22e+05       -j1.63e+05

从上两个表格可看出，文献方法的计算结果和本发明的方法结果一致，这里也验证了定理2(从电容到综合参数的修正)没有损失计算精度。

计算效率比较：

以上算例证明了方法是正确的。但是由于规模太小，耗时太少，不便于比较效率。

另一个较大的计算实例(变量数为7252)，提取任一频率下综合耦合参数，本发明中的方法只需要9.0秒，而文献ASP-DAC2006中方法则需要400秒以上，加速比在40倍以上。同时，本发明中方法需要60MB内存，而文献方法需要45MB内存，本发明中的方法需要较多内存是因为需要存储额外的矩阵，比如A_cap^-1B、M、W等，但这些额外内存的使用可以换来几十倍的加速比。