4.5 广义最小二乘法（GLS ） GLS----Generalized Least Squares 1. 基本原理广义最小二乘法的基本思想在于引入一个所谓成形滤波器（白化滤波器），把相关噪声)(k ξ转化成白噪声)(k ε。

由方程（4-4）、（4-5），系统的差分方程可以表示为)()()()()(11k k u z b k y z a ξ+=-- （4-114）式中n n z a z a z a z a ----++++=ΛΛ221111)(nn z b z b z b b z b ----++++=ΛΛ221101)(如果知道有色噪声序列)(k ξ的相关性，则可以把)(k ξ看成白噪声通过线性系统后所得的结果。

这种线性系统通常称为成形滤波器，其差分方程为)()()()(11_k z d k zc εξ---= （4-115）式中)(k ε是均值为零的白噪声序列，)()(11_---z d 、z c 是1-z 的多项式。

令 _111212_1()()1()m m c z f z f z f z f z d z ------==+++L L （4-116）有 )()(1)()()()(11k z f k k k z f εξεξ--==或 （4-117）即1212(1)()()m m f z f z f z k k ξε---++++=L L （4-118）或)()()2()1()(21k m k f k f k f k m εξξξξ+-------=ΛΛ ()1,,n k n N =++L L（4-119）这一噪声模型（自回归模型）的阶m ，一般事先是不知道的，实际经验表明，若指定m为2或3，就可以获得令人满意的描述)(k ξ的模型。

把方程（4-119）看作输入为零的差分方程，并由此式来写出N 个方程。

⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧++-+---+--+-=+++-+---+-=+++-+-----=+)()()2()1()()2()2()()1()2()1()1()1()()1(212121N n m N n f N n f N n f N n n m n f n f n f n n m n f n f n f n m m m εξξξξεξξξξεξξξξΛΛM ΛΛΛΛ写成向量矩阵形式为εξ+Ω=f （4-120）其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++=)()1(N n n ξξξM ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=m f f f M 1，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++=)()1(N n n εεεM ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-+--+--+--+--+--+----=Ω)()2()1()2()()1()1()1()(m N n N n N n m n n n m n n n ξξξξξξξξξM Λ(4-120)式所示的线性组合关系是辨识问题的基本表达形式,称作最小二乘格式。

应用我们熟知的最小二乘法，可求出f 的估值。

1ˆ[]T T fξ-=ΩΩΩ （4-121） 考虑在噪声模型（4-117）下，系统的差分方程。

将（4-117）式得到的)()(1)(1k z f k εξ-=，代入（4-114）式有)()(1)()()()(111k z f k u z b k y z a ε---+= （4-122）变换成 )()()()()()()(1111k k u z f z b k y z f z a ε+=---- （4-123）令 )()()(1k y k y z f =- （4-124）)()()(1k u k u z f =- （4-125）于是可得 )()()()()(11k k u z b k y z a ε+=-- （4-126）即1201()(1)(2)()()(1)()()n n y k a y k a y k a y k n b u k b u k b u k n k ε=-------++-++-+L L （4-127）在（4-126）或（4-127）式中，)(k ε为不相关的随机序列（白噪声），故可以用最小二乘法得到θ的无偏估计（即n o n b b b a a ΛΛ11,,,）。

由此可见，广义最小二乘法（GLS ）是建立在最小二乘法（LS ）的基础之上的。

【3】基本最小二乘法只是广义最小二乘法在1)(1=-z f 时的特例。

2. 计算步骤：广义最小二乘法的关键问题是如何用比较简便的方法找到成形滤波器的系数。

其计算是逐次逼近法。

下面以（4-121）和（4-123）式为依据来讨论。

第一步：应用输入、输出数据)2,1)(()(N n k k y k u +=Λ和，按最初的模型a ()()()()()k k ub k y ξ+Z =Z--11求出θ的最小二乘估计(1)(1)(1)T (1)(1)(1)(1)(1)(1)121ˆˆˆ[]ˆˆˆˆˆˆTn o n a b a a a b b b θ=⎡⎤=⋯⋯⋯⋯⎣⎦这个估值是不精确的，它只是被估参数的一次近似。

第二步：计算残差)()1(κe ，并拟合成形滤波器的模型：)()(ˆ)()(ˆ)(1)1(1)1()1(k u b k y ae --Z -Z =κ （4-128） 或 (1)(1)(1)(1)(1)(1)11ˆˆˆˆˆ()()(1)()()(1)()n o ne y a y k a y k n b u k b u k b u k n κκ=+-+⋯⋯+-----⋯⋯-- ()N n n k +⋯⋯+=,1 （4-129）)(k e 称为广义残差，用广义残差e )1(k 代替ζ（k ）注意：残差)(k e 是模型噪声，而()k ε是系统噪声，可用最小乘法拟合出一个成型滤波器的模型。

)()()1()()1()1(1)1(k m k e f k e f k e m ε=--⋯⋯--=N n n k ++=ΛΛ,1 （4-130）利用（4-121）式得到f 的最小二乘估值)1()1(1)1()1()1(][ˆe fT T ΩΩΩ=- （4-131） 其中 (1)(1)(1)(1)12ˆˆˆˆ[]T mf f f f =L LT N n e n e n e e ]2)1([)1()1()1()1(）（）（+++=ΛΛ(1)(1)(1)(1)(1)()(1)(1)()e n e n m e n N e n m N ⎡⎤--+⎢⎥Ω=⎢⎥⎢⎥-+--+⎣⎦L LMML L第三步：应用而所得的成形滤波器，对输入输出数据滤波由（4-124）（4-125）式有：⎪⎭⎪⎬⎫Z =Z =--)()(ˆ)()()(ˆ)(1)1()1(1)1()1(k u f k u k y f k y （4-132） 或 ⎪⎭⎪⎬⎫-++-+=-++-+=)()1(ˆ)()()(ˆ)1(ˆ)()()1()1(1)1()1()1(1)1(m k u f k u f k u k u m k y f k y f k y k y m mΛΛΛΛ （4-133） 第四步：求出参数θ的第二次估值θˆ（2）按模型a(z -1))()()()()1(1)1(k k u z b k yε+=-重新估计θ，得其最小=乘估值)2(ˆθ，并转入第二步。

以下将重复以上步骤：由)2(ˆθ，按步骤2计算残差()(2)k e，估值)2(ˆf ；按步骤3计算(2)(2)()();yk u k 、按步骤求θ的第3次估值)3(ˆθ。

进行多次循环，直到θ的第i 次估值)(ˆi θ收敛为止。

循环程序的收敛性判断：lim()1ˆ()1i i fz -→∞= （4-134） 这意味着残差e )(k 已经白噪声化，数据不需要继续滤波)(ˆi θ是参数θ的一个良好估计。

广义最小二乘辩识算法的程序框图：图4.12 广义最小二乘辩识算法的程序框图3．算例：将GLS 算法与LS 算法进行比较 有一单输入一单输出系统)()1()2()1()(121k k u b k y a k y a k y ξ+-+---=θ的真值为 []121[-0.5 0.5 1.0]T a a b θ==输入u （k ）是具有零均值和单位方差的独立高斯随机变量序列；)(k ξ 的成形滤波器模型为：11()()(1.00.85)()()f z k z k k ξξε--=+=即)(85.0,0.1110相关较大是为了使残差强烈f f f ==)(k ε 是具有方差20.64δ=的零均值白噪声。

辩识结果：利用了N=300的输入输出数据，用广义最小二乘法（GLS ）进行迭代计算。

每次迭代计算，都计算出残差的均方误差：Ne e T =2ˆσ，计算结果绘于图X 中。

图4.13 辨识结果图4. 优缺点：优：能够克服当存在有色噪声干扰时，基本最小二乘估计的有偏性，估计效果较好，在实际中得到较好的应用。

缺：——计算量大，每个循环要调用两次最小二乘法及一次数据滤波，——求差分方程（4-123）的参数估值，是一个非线性最优化问题，不一定总能保证算法对最优解的收敛性。

广义最小二乘法本质上是一种逐次逼近法。

对于循环程序的收敛性还没有给出证明。

——GLS 算法的最小二乘指标函数J 中可能存在一个以上局部极小值，（特别在信噪比不大时，J 可能是多举的）。

GLS 方法的估计结果往往取决于所选用参数的初始估值。

参数估计初值应选得尽量接近优参数。

在没有验前信息的情况下，最小二乘估值被认为是最好的初始条件。

——广义最小二乘法的收敛速度不是很高。

（图4.11即是一个说明）。

4.6 递推广义最小二乘法（RGLS ）RGLS---Recursive Generalized Least Squares 首先回顾整批GLS 辨识方法中的有关算式：11()()()()()a z y k b z u k k ξ--=+ （4-114）)()()(1k k z f εξ=- （4-117）)()()()()()()(1111k k u z f z b k y z f z a ε+=---- （4-123）)()()(1k y z f k y -= （4-124） )()()(1k u z f k u -= （4-125）)()()()()(11k k u z b k y z a ε+=-- （4-126）（4-126）式的形式和（4-4）、（4-5）式完全一样，相应有ξφθεθφ+=↔+=Y Y （4-33）（4-135）为建立递推计算公式，当观测数据长度为N 时，把（4-135）式写成：N N N N Y εθφ+= （4-136）广义最小二乘法的递推计算过程可分成二部分：（1）按递推最小二乘法（RLS ），随着N 的增大，不断计算N θˆ（逐步接近于无偏）和Nf ˆ（逐步使噪声白化）； （2）在递推过程中，N θˆ和Nf ˆ是时变的，则过滤信号)().(k y k u 及残差)(k e 是由时变系统产生，要不断计算)().(k y k u 及)(k e 。