振动对于航空、航天、船舶、机械、电子、大型桥梁等许多工程结构常常是造成恶性破坏的直接原因。因此，在大型结构设计中，振动设计与控制是至关重要的。对工程结构进行定量、准确的动力学分析，解决工程结构中普遍存在的振动问题，首先必须建立结构的动力学模型。有限元方法具有分析速度快、设计周期短、与结构试验相比成本低等优点，已成为结构动力系统建模最常用、最有效的方法之一。由有限元方法可建立无阻尼结构系统的运动微分方程

$ {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{\ddot x}}\left( t \right) + {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{x}}\left( t \right) = \mathit{\boldsymbol{f}}\left( t \right) $

式中：x(t)∈Rⁿ为位移向量；n为有限元模型的自由度；M_a, K_a∈R^n×n分别为系统有限元模型的质量和刚度矩阵; f(t)∈Rⁿ为载荷向量。通常，有限元模型的质量和刚度矩阵都是对称半正定矩阵，并且具有稀疏结构。结构系统的固有频率ω或特征值λ(λ=ω²)以及相应的振型或特征向量x由如下特征方程确定

$ {\mathit{\boldsymbol{K}}_a}\mathit{\boldsymbol{x}} = \lambda {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{x}} $

(1)

应用特征值问题的数值方法^[1]，可求得广义特征值问题(1)的特征值λ_i^(a)和相应的特征向量x_i^(a)。

随着试验模态分析技术、信号处理和动态测试技术的发展，由振动试验可获得实际结构系统的一些低阶频率ω_i^(e)(i=1, 2, …, m≤n)，即特征值λ_i^(e)=(ω_i^(e))²以及相应的振型x_i^(e)。由于实际结构的复杂性，有限元方法必须对结构的边界条件、连接条件和约束条件进行简化处理，还必须对结构的几何特性做一些假设，因而有限元方法所建立的结构有限元模型往往不能准确地反映实际结构的动态特性，有限元模型的计算结果λ_i^(a)和x_i^(a)与实测结果λ_i^(e)和x_i^(e)之间存在着误差。另一方面，受试验条件和成本的限制，实测信息的不完整导致试验建模无法获得完备的数学模型。有限元模型修正就是将有限元建模和试验建模结合起来，先用有限元方法建立结构的有限元模型，然后依据系统动态特性的试验测量数据修正有限元模型，使得修正后结构模型的低阶动态特性参数与试验值趋于一致，并且使有限元模型在保持拓扑结构的同时改变量尽可能小。数学上，有限元模型修正问题归结为求矩阵，使其满足低阶动态特性、保持物理和结构性质的同时，并使矩阵的改变量达到最小。修正的有限元动力学模型可用于结构响应分析和控制。

四十多年来，结构有限元模型修正的理论和方法取得了许多成果^[2-3]，并且在工程实际中得到了重要应用^[4-5]。因为有限元模型质量和刚度矩阵是实际结构的较好近似，所以最早发展并被广泛使用的有限元模型修正方法是以系统矩阵或矩阵元素为修正对象的矩阵型和元素型修正方法。Berman等^[6-7]利用测量的部分振型修正质量矩阵M，要求修正的质量矩阵M满足对称性M^T=M和正交性条件X_e^TMX_e=I_m，并使质量矩阵的改变量最小，应用Lagrange乘子法，导出了修正质量矩阵M的表达式。Zhang和Lee等^[8-9]分别利用矩阵变换和Moore-Penrose广义逆给出了修正质量矩阵M的显式表示。但这些方法修正的质量矩阵M不仅改变了解析质量矩阵M_a的带状或稀疏结构，而且未必是半正定的。为保持修正质量矩阵的稀疏结构，Wei和Zhang^[10], Cha^[11]以质量矩阵的待修正元素为对象, 发展了修正质量矩阵的元素型修正方法。Yuan和Dai^[12]考虑了解析质量矩阵的一个主子矩阵准确、修正其余元素的问题，利用广义奇异值分解给出了修正质量矩阵的表达式。这类方法具有灵活、方便等特点，但难以保证修正质量矩阵M的半正定性。

为了保证修正的质量矩阵既满足对称半正定性，又具有稀疏结构，本文以质量矩阵改变量的F-范数为目标函数，并以待修正质量矩阵应具有的性质，如满足正交关系或特征方程，对称性，半正定性和稀疏性作为约束条件，数学上将问题归结为带约束的矩阵最佳逼近问题。应用矩阵的QR分解，给出了该问题有解的条件。基于循环投影方法^[13]，提出了求解这类问题的数值方法，提供了两个数值例子以说明所给方法的有效性。

为了表述方便，对使用的记号作如下说明:rank(A)和A⁺分别表示矩阵A的秩和Moore-Penrose广义逆。对A, B∈R^m×n, (A, B)=tr(B^TA)和A×B分别表示矩阵A与B的内积和Hadamard积，基于此内积，R^m×n构成Hilbert空间，并且由此内积导出的长度即为矩阵的F-范数‖*‖_F。SR^n×n表示n阶实对称矩阵的全体，SR₀^n×n表示n阶对称半正定矩阵的集合，OR^n×n表示n阶正交矩阵的全体。对A∈SR^n×n，A≥0表示A为实对称半正定矩阵。sparse(M)=sparse(M_a)表示矩阵M和M_a具有相同的稀疏结构。I_n表示n阶单位矩阵。记Λ_e=diag(λ₁^(e), …, λ_m^(e))，X_e=[x₁^(e), …, x_m^(e)]，并设Λ_e>0，rank(X_e)=m，M_a=(m_aij), K_a=(k_aij)∈SR₀^n×n，并且M_a正定。

1 利用正交关系修正质量矩阵

为了使修正质量矩阵M既满足如下正交性条件

$ \mathit{\boldsymbol{X}}_e^T\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m} $

(2)

又是对称半正定的，并保持稀疏结构，本文将修正质量矩阵M的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\mathit{\boldsymbol{X}}_e^T\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0}\\ {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $

(3)

设测量的振型矩阵X_e的QR分解为

$ {\mathit{\boldsymbol{X}}_e} = \mathit{\boldsymbol{Q}}\left( \begin{array}{l} \mathit{\boldsymbol{R}}\\ {\bf{0}} \end{array} \right) $

(4)

式中：Q=[Q₁, Q₂]∈OR^n×n, Q₁∈R^n×m，R为m阶可逆上三角矩阵。将式(4)代入式(2)，可导出满足正交条件(2)的对称矩阵为

$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{22}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

(5)

式中：M₁₂∈R^m×(n－m)和M₂₂∈SR^{(n－m)×(n－m)}是任意的。

为了求得矩阵方程(2)的对称半正定解，需要如下引理。

引理1^[14] 设E∈SR^n×n, F∈R^n×k, G∈SR^k×k，$ \boldsymbol{H }= \left( {\begin{array}{*{20}{c}} \boldsymbol{E}&\boldsymbol{F}\\ {{\boldsymbol{F}^{\rm{T}}}}&\boldsymbol{G} \end{array}} \right)$，则H≥0当且仅当E≥0, G－F^TE⁺F≥0，且rank[E, F]=rank(E)。

由引理1，即得矩阵方程(2)的对称半正定解M为

$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{\mathit{\boldsymbol{M}}_{12}^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{12}} + \mathit{\boldsymbol{G}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

(6)

式中：M₁₂∈R^m×(n－m)和G∈SR₀^{(n－m)×(n－m)}是任意的。

为了给出问题(3)解的存在性和唯一性，需要如下引理。

引理2(最佳逼近定理)^[15] 设H是一个Hilbert空间，W是H中的一个非空闭凸集，则对g∈H，存在唯一的向量g^*∈W，使得

$ \left\| {\mathit{\boldsymbol{g}} - {\mathit{\boldsymbol{g}}^ * }} \right\| = \mathop {\min }\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{W}}} \left\| {\mathit{\boldsymbol{g}} - \mathit{\boldsymbol{x}}} \right\| $

(7)

式中：x表示Hilbert空间中向量x的长度或由长度导出的范数。

问题(7)的解g^*通常称为g在W中的最佳逼近或称g^*为g在W中的投影，记为g^*=P_W(g)。

如果引理2中Hilbert空间H的非空闭凸集W= $\bigcap\limits_{i = 1}^k {{\boldsymbol{W}_i}} $，其中W_i(i=1, 2, …, k)都是Hilbert空间H中的闭凸集，并且对g∈H，仅能计算g在每个W_i上的投影P_Wi(g)。Boyle和Dykstra^[13]给出了计算g在W中投影的一个迭代方法，称为循环投影方法，并证明了循环投影方法的收敛性。

由引理2可得如下结论。

定理1  设测量的振型矩阵X_e的QR分解为式(4)，如果存在M₁₂∈R^m×(n－m)和G∈SR₀^{(n－m)×(n－m)}，使得式(6)中的质量矩阵M满足sparse(M)=sparse(M_a)，则问题(3)的可行域是非空闭凸集，从而问题(3)存在唯一解。

记

$ {S_{{M_1}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m}} \right.,{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}}} \right\} $

$ {S_{{M_2}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \right.} \right\} $

$ {S_{{M_3}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $

则S_M3=SR₀^n×n，问题(3)的可行域为S_M1∩S_M2∩S_M3，并且是R^n×n中的闭凸集。因此问题(3)等价于求M_a在S_M1∩S_M2∩S_M3中的最佳逼近。下面依次考虑M_a在S_M1, S_M2, S_M3上的投影。

对解析质量矩阵M_a∈SR^n×n，记

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{Q}} = \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_{a11}}}&{{\mathit{\boldsymbol{M}}_{a12}}}\\ {\mathit{\boldsymbol{M}}_{a12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{a22}}} \end{array}} \right)}\\ {{\mathit{\boldsymbol{M}}_{aij}} = \mathit{\boldsymbol{Q}}_i^{\rm{T}}{\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{Q}}_j}\;\;\;\;i,j = 1,2} \end{array} $

(8)

对M∈S_M1，由式(5，8)和F-范数的正交不变性，可得

$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2 = \left\| {\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{22}}} \end{array}} \right) - {\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{Q}}} \right\|_{\rm{F}}^2 = }\\ {\left\| {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}} - {\mathit{\boldsymbol{M}}_{a11}}} \right\|_{\rm{F}}^2 + 2\left\| {{\mathit{\boldsymbol{M}}_{12}} - {\mathit{\boldsymbol{M}}_{a12}}} \right\|_{\rm{F}}^2 + }\\ {\left\| {{\mathit{\boldsymbol{M}}_{22}} - {\mathit{\boldsymbol{M}}_{a22}}} \right\|_{\rm{F}}^2} \end{array} $

当且仅当M₁₂=M_a12, M₂₂=M_a22，$ \mathop {\min }\limits_{\boldsymbol{M} \in \boldsymbol{S}{\boldsymbol{R}^{n \times n}}} \boldsymbol{M} - \boldsymbol{M}_{a{\rm{F}}}^2$取得最小值。因此, 解析质量矩阵M_a在S_M1中的投影为

$ {P_{{S_{{M_1}}}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{a12}}}\\ {\mathit{\boldsymbol{M}}_{a12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{a22}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

由式(4)可得X_e⁺=R^－1Q₁^T。利用Moore-Penrose广义逆的性质和式(8)，上式可化为

$ {P_{{S_{{M_1}}}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = {\mathit{\boldsymbol{M}}_a} + {\left( {\mathit{\boldsymbol{X}}_e^ + } \right)^{\rm{T}}}\mathit{\boldsymbol{X}}_e^ + - {\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{X}}_e^ + {\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{X}}_e^ + $

(9)

对给定的解析质量矩阵M_a∈SR^n×n，使用如下指标集S_Ma刻画矩阵M_a的稀疏结构

$ \begin{array}{*{20}{c}} {{S_{{M_a}}} = \left\{ {\left( {i,j} \right)\left| {{m_{aij}} = 0,} \right.} \right.}\\ {\left. {{\mathit{\boldsymbol{M}}_a} = \left( {{m_{aij}}} \right) \in \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{R}}^{n \times n}}} \right\}} \end{array} $

(10)

由指标集S_Ma的定义，容易证明如下结论。

引理3  设T=(t_ij)∈SR^n×n，则问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - \mathit{\boldsymbol{T}}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $

(11)

有唯一解M=(m_ij)∈SR^n×n，且m_ij= $\left\{ \begin{array}{l} 0,\left( {i,j} \right) \in {S_{{M_a}}}\\ {t_{ij}},\left( {i,j} \right) \notin {S_{{M_a}}}。\end{array} \right. $

由引理3可知，问题(11)的M就是T在S_M2中的投影，即

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}} = {P_{{S_{{M_2}}}}}\left( \mathit{\boldsymbol{T}} \right) = \left( {{m_{ij}}} \right)}\\ {{m_{ij}} = \left\{ \begin{array}{l} 0\;\;\;\;\left( {i,j} \right) \in {S_{{M_a}}}\\ {t_{ij}}\;\;\;\left( {i,j} \right) \notin {S_{{M_a}}} \end{array} \right.} \end{array} $

(12)

为了计算一个矩阵在S_M3中的投影，需要如下引理。

引理4^[16]  设A∈R^n×n，$ \boldsymbol{\hat A} = \frac{1}{2}$(A+A^T)的谱分解为$\boldsymbol{\hat A} $=Udiag(θ₁, …, θ_n)U^T，其中θ₁, …, θ_n为矩阵$\boldsymbol{\hat A} $的特征值, U∈OR^n×n，则存在唯一的对称半正定矩阵A₊，使得

$ {\left\| {{\mathit{\boldsymbol{A}}_ + } - \mathit{\boldsymbol{A}}} \right\|_{\rm{F}}} = \mathop {\min }\limits_{\mathit{\boldsymbol{X}} \in \mathit{\boldsymbol{SR}}_0^{n \times n}} {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{A}}} \right\|_{\rm{F}}} $

(13)

并且A₊=Udiag(β₁, …, β_n)U^T，其中β_i=max{θ_i, 0}(i=1, 2, …, n)。

对矩阵A∈R^n×n，问题(13)的解称为矩阵A的正逼近，记为A₊。实际上，A₊是A∈R^n×n在S_M3中的投影，即A₊=P_SM3(A)=P_SR0^n×n(A)。

于是，可给出求解问题(3)的数值算法如下：

算法1  修正质量矩阵的循环投影方法

输入：M_a∈SR^n×n，X_e∈R^n×m和收敛准则ε，并置M_{1, 0}=M_a，I_{0, j}=0(j=1, 2, 3)。

输出：修正的质量矩阵M。

对i=1, 2, …，执行以下运算直到收敛

对j=1, 2, 3，计算

$ {\mathit{\boldsymbol{M}}_{i,j}} = {P_{{S_{{M_j}}}}}\left( {{\mathit{\boldsymbol{M}}_{i,j - 1}},{\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

$ {\mathit{\boldsymbol{I}}_{i,j}} = {\mathit{\boldsymbol{M}}_{i,j}} - \left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

置M_{i+1, 0}=M_{i, 3}, 如果‖M_{i+1, 0}－M_{i, 0}‖_F≤ε，则M=M_{i+1, 0}，停止；否则继续迭代。

定理2^[13]  如果问题(3)的可行域非空，则对任意j(j=1, 2, 3)，算法1所产生的矩阵序列{M_{n, j}}_n=1^∞收敛于问题(3)的唯一解。

由算法1修正的质量矩阵不仅满足正交性条件(2)，而且是对称半正定的，与解析质量矩阵M_a具有相同的稀疏结构。

2 利用特征方程修正质量矩阵

假定解析刚度矩阵K_a准确，并设测量的特征值Λ_e和相应的振型X_e满足

$ \mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} $

(14)

基于待修正质量矩阵M应满足如下特征方程

$ \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e} $

(15)

和对称半正定性，戴华^[17]将修正质量矩阵的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \end{array} $

(16)

记

$ {D_M} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $

则D_M是问题(16)的可行域。

为了给出D_M中元素的表达式，需要如下引理。

引理5^[18]  设X, B∈R^n×m，如果rank(X)=m，并且X有QR分解X=Q$\left( {\begin{array}{*{20}{c}} \boldsymbol{R}\\ \boldsymbol{0} \end{array}} \right) $，其中Q=[Q₁, Q₂]∈OR^n×n, Q₁∈R^n×m，R为m阶可逆上三角矩阵，则AX=B有解A∈SR₀^n×n当且仅当

$ {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{X}} \ge 0,{\rm{rank}}\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right) = {\rm{rank}}\left( \mathit{\boldsymbol{B}} \right)。$

在有解的情况下，AX=B的对称半正定解可表示为A=BX⁺+[(BX⁺)^T+(I_n－XX⁺)B(X^TB)⁺B^T](I_n－XX⁺)+Q₂GQ₂^T，其中G∈SR₀^{(n－m)×(n－m)}是任意的。

定理3  设测量的特征值Λ_e和相应的振型X_e满足式(14)，且X_e的QR分解为式(4)，则问题(16)的可行域非空，并且问题(16)存在唯一解

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}{\left[ {\mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right]_ + }\mathit{\boldsymbol{Q}}_2^{\rm{T}} $

(17)

其中M₀=K_aX_eΛ_e^－2X_e^TK_a。

证明  由引理5知，矩阵方程MX_eΛ_e=K_aX_e存在对称半正定解，并且其通解可表示为

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e^{ - 2}\mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{K}}_a} + {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{GQ}}_2^{\rm{T}} $

(18)

式中：G∈SR₀^{(n－m)×(n－m)}是任意的，从而问题(16)的可行域是非空闭凸集。由引理2可知，问题(16)存在唯一解。

对任意M∈D_M，M可表示为

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{GQ}}_2^{\rm{T}} = {\mathit{\boldsymbol{M}}_0} + \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

由F-范数的正交不变性，可得

$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2 = \left\| {\mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} - \left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right)} \right\|_{\rm{F}}^2}=\\ {\left\| {\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right) - {\mathit{\boldsymbol{Q}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right)\mathit{\boldsymbol{Q}}} \right\|_{\rm{F}}^2 = }\\ {\left\| {\mathit{\boldsymbol{Q}}_1^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_1}} \right\|_{\rm{F}}^2 + 2\left\| {\mathit{\boldsymbol{Q}}_1^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right\|_{\rm{F}}^2 + }\\ {\left\| {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right\|_{\rm{F}}^2} \end{array} $

从而║M－M_a║_F²=min当且仅当║G－Q₂^T(M_a－M₀)Q₂║_F²=min。由引理4可知，当G=[Q₂^T(M_a－M₀)Q₂]₊时，║M－M_a║_F²取最小值。将G=[Q₂^T(M_a－M₀)Q₂]₊代入式(18)即得式(17)。

定理3给出的修正质量矩阵M不仅满足特征方程(15)，而且是对称半正定矩阵。袁永新和戴华^[19]利用奇异值分解给出了问题(16)解的表达式。但这些表达式不能保持解析质量矩阵的稀疏结构。

为了使修正质量矩阵M既满足特征方程(15)，又是对称半正定的，并保持稀疏结构，本文将修正质量矩阵M的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0}\\ {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $

(19)

由定理3和引理2，可得如下结论。

定理4  设测量的特征值Λ_e和振型X_e满足式(14)，且X_e的QR分解为式(4)。如果存在矩阵G∈SR₀^{(n－m)×(n－m)}，使得M=K_aX_eΛ_e^－2X_e^TK_a+Q₂GQ₂^T满足sparse(M)=sparse(M_a)，则问题(19)的可行域非空，从而问题(19)存在唯一解。

记

$ {S_1} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $

$ {S_2} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \right.} \right\} $

则问题(19)的可行域D_M=S₁∩S₂。

对给定的矩阵M_a∈SR^n×n，由定理3可得M_a在S₁上的投影为

$ {P_{{S_1}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}{\left[ {\mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right]_ + }\mathit{\boldsymbol{Q}}_2^{\rm{T}} $

(20)

因为S₂=S_M2，对给定的矩阵T=(t_ij)∈SR^n×n，由式(12)容易计算T在S₂上的投影P_S2(T)=P_SM2(T)。

于是，给出求解问题(19)的数值算法如下：

算法2  修正质量矩阵的循环投影方法

输入：K_a, M_a∈SR^n×n, X_e∈R^n×m, Λ_e∈R^m×m和收敛准则ε，置M_{1, 0}=M_a，I_{0, j}=0(j=1, 2)。

输出：修正的质量矩阵M。

对i=1, 2, …，执行以下运算直到收敛。

对j=1, 2，计算

$ {\mathit{\boldsymbol{M}}_{i,j}} = {P_{{S_j}}}\left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

$ {\mathit{\boldsymbol{I}}_{i,j}} = {\mathit{\boldsymbol{M}}_{i,j}} - \left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

置M_{i+1, 0}=M_{i, 2}, 如果‖M_{i+1, 0}－M_{i, 0}‖_F≤ε，则M=M_{i+1, 0}，停止；否则继续迭代。

由文献[13]中定理2，可得算法2的收敛性结果。

定理5  如果问题(19)的可行域非空，则对任意j(j=1, 2)，算法2所产生的矩阵序列{M_{n, j}}_n=1^∞收敛于问题(19)的唯一解。

3 数值结果

本节给出两个数值例子，以说明算法1和算法2的有效性。算法1和算法2用MATLAB R2015b编程, 程序在主频为2.90 GHz的计算机上实现。收敛准则ε=10^－7。用M和$\boldsymbol{\tilde M} $分别表示精确的质量矩阵和计算的修正质量矩阵，Iter和CPU分别表示执行算法的迭代步数和计算时间(单位: s)，记

$ {\mathit{\boldsymbol{E}}_1} = {\left\| {\mathit{\boldsymbol{X}}_e^{\rm{T}}\mathit{\boldsymbol{\tilde M}}{\mathit{\boldsymbol{X}}_e} - \mathit{\boldsymbol{I}}} \right\|_{\rm{F}}} $

$ {\mathit{\boldsymbol{E}}_2} = {\left\| {\mathit{\boldsymbol{\tilde M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} - {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e}} \right\|_{\rm{F}}},{\rm{Err}} = \frac{{{{\left\| {\mathit{\boldsymbol{\tilde M}} - \mathit{\boldsymbol{M}}} \right\|}_{\rm{F}}}}}{{{\mathit{\boldsymbol{M}}_{\rm{F}}}}} $

例1  考虑图 1所示的无阻尼弹簧-质点系统，其中m_i=1 kg(i=1, 2, …, 5)，k_i=0.5 N/m(i=1, 2, …, 6)，则系统精确的质量矩阵与刚度矩阵分别为

$ \mathit{\boldsymbol{M}} = {\rm{diag}}\left( {1,1,1,1,1} \right) $

$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{K}}_a} = \left( {\begin{array}{*{20}{c}} 1&{ - 0.5}&{}&{}&{}\\ { - 0.5}&1&{ - 0.5}&{}&{}\\ {}&{ - 0.5}&1&{ - 0.5}&{}\\ {}&{}&{ - 0.5}&1&{ - 0.5}\\ {}&{}&{}&{ - 0.5}&1 \end{array}} \right) $

取广义特征值问题Kx=λMx的二个最小特征值及其相应的特征向量作为测量特征对, 分别记为Λ_e, X_e, 即

$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\rm{diag}}\left( {0.134\;0,0.500\;0} \right) $

$ {\mathit{\boldsymbol{X}}_e} = \left( {\begin{array}{*{20}{c}} {0.288\;7}&{ - 0.500\;0}\\ {0.500\;0}&{ - 0.500\;0}\\ {0.577\;4}&0\\ {0.500\;0}&{0.500\;0}\\ {0.288\;7}&{0.500\;0} \end{array}} \right) $

容易验证：Λ_e和X_e满足式(14)。M的估计为M_a=M+μE_r×M，其中E_r为随机生成的五阶对称矩阵, 其元素随机分布在区间[-1.0, 1.0], 参数μ表示对M相对扰动的大小。利用算法1和算法2计算的修正质量矩阵$\boldsymbol{\tilde M} $与精确的质量矩阵M完全相同。数值结果见表 1。

例2  考虑一端固定、另一端自由的离散变截面杆^[20]，基于线性形函数的解析质量矩阵与刚度矩阵为

$ \mathit{\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} {0.4}&{0.1}&{}&{}&{}&{}&{}&{}&{}&{}\\ {0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}&{}&{}\\ {}&{0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}&{}\\ {}&{}&{0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}\\ {}&{}&{}&{0.1}&{0.6}&{0.2}&{}&{}&{}&{}\\ {}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}&{}&{}\\ {}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}&{}\\ {}&{}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}\\ {}&{}&{}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}\\ {}&{}&{}&{}&{}&{}&{}&{}&{0.2}&{0.4} \end{array}} \right) $

$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{K}}_a} = \left( {\begin{array}{*{20}{c}} 2&{ - 1}&{}&{}&{}&{}&{}&{}&{}&{}\\ { - 1}&2&{ - 1}&{}&{}&{}&{}&{}&{}&{}\\ {}&{ - 1}&2&{ - 1}&{}&{}&{}&{}&{}&{}\\ {}&{}&{ - 1}&2&{ - 1}&{}&{}&{}&{}&{}\\ {}&{}&{}&{ - 1}&4&{ - 3}&{}&{}&{}&{}\\ {}&{}&{}&{}&{ - 3}&6&{ - 3}&{}&{}&{}\\ {}&{}&{}&{}&{}&{ - 3}&6&{ - 3}&{}&{}\\ {}&{}&{}&{}&{}&{}&{ - 3}&6&{ - 3}&{}\\ {}&{}&{}&{}&{}&{}&{}&{ - 3}&6&{ - 3}\\ {}&{}&{}&{}&{}&{}&{}&{}&{ - 3}&3 \end{array}} \right) $

这里刚度矩阵和质量矩阵的元素缩放了一定比例。取广义特征值问题Kx=λMx的4个最小特征值及其相应的特征向量作为测量特征对, 分别记为Λ_e, X_e。M的估计为M_a=M+μE_r×M，其中E_r为随机生成的10阶对称矩阵, 其元素随机分布在区间[-1.0, 1.0]。利用算法1和算法2计算的修正质量矩阵$\boldsymbol{\tilde M} $与精确的质量矩阵M完全相同。数值结果见表 2。

例1, 2的数值结果说明由算法1和算法2修正的质量矩阵不仅满足正交性条件或特征方程, 而且是对称半正定矩阵, 还具有要求的稀疏结构。因此，算法1, 2是有效的，并且由表 1, 2可见，如果刚度矩阵已知，由特征方程修正质量矩阵的计算效率更高。

4 结束语

本文研究了有限元模型质量矩阵的修正问题。建立了以待修正质量改变量的F-范数为目标函数并以质量矩阵应具有的性质如满足正交关系或特征方程，对称性，半正定性和稀疏性为约束条件的矩阵最佳逼近问题。给出了问题有解的条件。提出了求解带约束矩阵最佳逼近问题的循环投影方法。所获得的修正质量矩阵不仅满足正交性条件或特征方程，而且是对称半正定的，并与解析质量矩阵具有相同的稀疏结构。本文所发展的修正质量矩阵的理论和方法可应用于修正有限元模型的刚度矩阵或同时修正有限元模型的质量和刚度矩阵。