对线性回归分析、坐标转换和自回归分析等问题进行建模时，往往需要根据观测数据的实际情况对模型参数进行选择，选取的模型参数过少会导致所建模型带有系统误差，而参数过多会在一定程度上降低解算精度^[1]。由这几类问题的数学模型特点可以发现，对于平差模型优选，均可归结为最佳线性回归模型的选择问题。对于经典线性回归模型，线性假设法^[2-3]是目前国内外公认的模型优选有效方法，如利用其确定多项式回归模型阶数^[4]、判断坐标转换模型的相似变换条件是否成立^[5]、检验三维坐标转换中尺度参数的显著性^[6]以及确定一维自回归模型和多维自回归模型阶数^[7-8]等。研究表明，利用线性假设法需满足一定假设条件：每次被检验的两个模型中至少存在一个符合观测数据实际，否则得到的结果不可靠^[9-10]。但在实际建模时，每次检验的两个模型是否符合观测数据实际往往未知，因此，采用该方法解决含有多个待选模型的最佳模型选择问题时，结果的可靠性和有效性有时无法保证^[11]。

此外，经典线性回归模型认为仅因变量包含测量误差，而忽略自变量的测量误差。但在实际应用中，回归模型的自变量和因变量均源自观测数据(或测量平差值)，不可避免地含有误差。因此，传统回归分析基于经典线性回归模型并采用最小二乘法求解回归系数，其解具有有偏性^[12]。对此，同时顾及自变量和因变量测量误差的回归模型得到发展，采用整体最小二乘法求解回归系数^[13-15]，其解具有渐进无偏性^[16]。这些研究成果均以给定的平差模型为基础，而在实际应用中，利用这些方法往往应先考虑模型的适用性问题，只有选用合理的平差模型才能够准确表达观测量之间的物理或几何关系。

目前，除线性假设检验法外，还有白噪声检验准则^[10]、残差平方和准则^{[1, 17]}、赤池弘治信息准则^[18]等模型优选方法，但这些方法同样仅适用于经典线性回归模型，无法直接用于自变量和因变量均含误差的回归模型。由此可见，研究有效的适用于同时考虑自变量和因变量测量误差的最佳线性回归模型的选择理论和方法，是亟待解决的问题。为此，文本基于同时考虑自变量和因变量测量误差的线性回归模型，结合附有参数约束的测量平差理论和含有多个备选假设的假设检验理论，研究线性回归模型的优选方法。

1 附有约束条件的线性回归模型 1.1 顾及自变量和因变量测量误差的线性回归模型

线性回归分析、坐标转换和自回归分析等问题的平差模型均可看作线性回归模型。在传统线性回归分析模型基础上，为自变量引入误差向量，可得：

$ {\boldsymbol{y}} = {\boldsymbol{\bar A\xi }} + {{\boldsymbol{e}}_y} = {\boldsymbol{A\xi }} + {{\boldsymbol{E}}_A}{\boldsymbol{\xi }} + {{\boldsymbol{e}}_y} $

(1)

式中，y为n×1观测值向量；A为n×t系数矩阵真值矩阵；ξ为t×1待估模型参数；e_y为n×1随机误差向量；A为n×t系数矩阵；E_A为A对应的n×t改正数矩阵，其中，固定值对应部分为0，相同随机元素对应的改正数相同。

式(1)中各项误差的随机模型为：

$ \left[ \begin{array}{l} {{\boldsymbol{e}}_{As}}\\ {{\boldsymbol{e}}_y} \end{array} \right]\sim \left( {\left[ \begin{array}{l} {\boldsymbol{0}}\\ {\boldsymbol{0}} \end{array} \right], \;{\sigma _0}\left[ \begin{array}{l} {{\boldsymbol{Q}}_{As}}\;\;\;{\boldsymbol{0}}\\ {\boldsymbol{0}}\;\;\;\;\;{{\boldsymbol{Q}}_y} \end{array} \right]} \right) $

(2)

式中，σ₀为单位权中误差；Q_y为观测量误差的协因数矩阵；设A中含有u个不同的随机量，a为将A中不同的随机量按列依次提取组成的u×1向量，e_As为a对应的误差向量，Q_As为a的协因数阵。

式(1)可视为非线性高斯-赫尔默特模型。将E_A中随机量作为待估参数，式(1)可表示为：

$ \left\{ \begin{array}{l} {\boldsymbol{a}} = {\boldsymbol{\bar a}} + {{\boldsymbol{e}}_{As}}\\ {\boldsymbol{y}} = {\boldsymbol{A\xi }} + {{\boldsymbol{E}}_A}{\boldsymbol{\xi }} + {{\boldsymbol{e}}_y} \end{array} \right.$

(3)

式中，a为a的真值向量，E_A为系数矩阵A对应的改正数矩阵，观测值为(a^T  y^T)^T，待估参数为(a^T  ξ^T)^T。

显然，式(3)为非线性模型，将其转化为线性模型并采用高斯-牛顿法迭代求解。设第i次迭代后，参数ξ的估值为ξ_(i)，e_As的估值为e_As(i)，将式(3)右端在(ξ_(i)，e_As(i))处采用泰勒级数展开并取至一阶项：

$ \left\{ \begin{array}{l} {\boldsymbol{a}} = {\boldsymbol{\bar a}} + {{\boldsymbol{e}}_{As}}\\ {\boldsymbol{y}} = {\boldsymbol{A}}{{\boldsymbol{\xi }}_{(i)}} + {{\boldsymbol{E}}_{A(i)}}{{\boldsymbol{\xi }}_{(i)}} + {\boldsymbol{A}}δ{\boldsymbol{\xi }} + {{\boldsymbol{E}}_{A(i)}}δ{\boldsymbol{\xi }} + {\boldsymbol{R}}{{\boldsymbol{e}}_{As}} + {{\boldsymbol{e}}_y} \end{array} \right. $

(4)

式中，δξ为ξ_(i)的微小改正值；R为与ξ_(i)及E_A有关的n×u矩阵，满足Re_As=E_Aξ_(i)。

构造目标函数：

$ \varphi \left( {{{\boldsymbol{e}}_y}, \;{{\boldsymbol{e}}_{As}}} \right) = {\boldsymbol{e}}_y^{\mathop{\rm T}\nolimits} {\boldsymbol{Q}}_y^{ - 1}{{\boldsymbol{e}}_y} + {\boldsymbol{e}}_{As}^{\mathop{\rm T}\nolimits} {\boldsymbol{Q}}_{As}^{ - 1}{{\boldsymbol{e}}_{As}} $

(5)

对各变量求偏导并令导数为0，可得：

$ {{\boldsymbol{\hat \beta }}_{(i + 1)}} = {\left( {{\boldsymbol{B}}_{(i)}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{B}}_{(i)}}} \right)^{ - 1}}{\boldsymbol{B}}_{(i)}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{l}}_{(i)}} $

(6)

$ {{\boldsymbol{v}}_{(i + 1)}} = {{\boldsymbol{B}}_{(i)}}{\left( {{\boldsymbol{B}}_{_{(i)}}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{B}}_{(i)}}} \right)^{ - 1}}{\boldsymbol{B}}_{_{(i)}}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{l}}_{(i)}} - {{\boldsymbol{l}}_{(i)}} $

(7)

式中，$ {{\boldsymbol{\hat \beta }}_{(i + 1)}} = \left( \begin{array}{l} {{\boldsymbol{e}}_{As(i + 1)}}\\ \delta {{\boldsymbol{\xi }}_{(i + 1)}} \end{array} \right)$, $ {{\boldsymbol{B}}_{(i)}} = \left( \begin{array}{l} {\boldsymbol{E}}\;\;\;\;{\boldsymbol{0}}\\ {\boldsymbol{R}}\;\;{{\boldsymbol{A}}_{(i)}} \end{array} \right)$, ${{\boldsymbol{l}}_{(i)}} = \left( \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{0}}\\ {\boldsymbol{y}} - \left( {{\boldsymbol{A}}{{\boldsymbol{\xi }}_{(i)}} + {{\boldsymbol{E}}_{A(i)}}{{\boldsymbol{\xi }}_{(i)}}} \right) \end{array} \right) $, $ {\boldsymbol{Q}} = \left( \begin{array}{l} {{\boldsymbol{Q}}_{As}}\\ \;\;\;\;\;\;{{\boldsymbol{Q}}_y} \end{array} \right)$, ${{\boldsymbol{v}}_{(i + 1)}} = \left( \begin{array}{l} {{\boldsymbol{e}}_{As(i + 1)}}\\ {{\boldsymbol{e}}_{y(i + 1)}} \end{array} \right) $。

以参数的最小二乘解为初值，根据式(6)和式(7)进行迭代计算，直至‖δξ_(i+1)‖₂＜ε(ε为很小的常数)时迭代结束。

1.2 附有参数约束的线性回归模型

在线性回归分析、坐标转换和自回归分析等问题对应的回归模型中，参数个数少的均可通过对参数个数多的附加一定的参数约束得到。因此，建模时若给定模型参数选择的最大范围，并称待选模型中参数个数最多的为一般回归模型(general regression model, GRM)，则其他待选模型均可通过对GRM依次添加相应参数约束得到。

假定对某实际问题，共有(f+1)个待选回归模型，根据模型中未知参数个数由多至少依次排序，设其分别为F₁= 0、F₂= 0 …F_f= 0、F_(f+1)= 0。选择参数个数最多的F₁= 0为GRM，则第j个待选模型F_j= 0可通过在式(3)基础上增加一个约束方程G_jξ=b_j得到，即

$ \left\{ \begin{array}{l} {\boldsymbol{a}} = {\boldsymbol{\bar a}} + {{\boldsymbol{e}}_{As}}\\ {\boldsymbol{y}} = {\boldsymbol{A\xi }} + {{\boldsymbol{E}}_A}{\boldsymbol{\xi }} + {{\boldsymbol{e}}_y}\\ {{\boldsymbol{G}}_j}\left( {\boldsymbol{\xi }} \right) = {{\boldsymbol{b}}_j} \end{array} \right. $

(8)

式中，j=2, 3, …, (f+1)；G_jξ=b_j为第j个待选模型对应的参数约束条件。

由式(4)、式(6)和式(7)迭代求解无约束线性回归模型的系数矩阵误差和未知参数值，设在第m次迭代后参数向量满足收敛条件。

结合式(4)，将式(8)在(ξ_(m)，e_As(m))处采用泰勒级数展开并取至一阶项：

$ \left\{ \begin{array}{l} {\boldsymbol{a}} = {\boldsymbol{\bar a}} + {{\boldsymbol{e}}_{As}}\\ {\boldsymbol{y}} = {\boldsymbol{A}}{{\boldsymbol{\xi }}_{(m)}} + {{\boldsymbol{E}}_{A(m)}}{{\boldsymbol{\xi }}_{(m)}} + {\boldsymbol{A}}δ{\boldsymbol{\xi }} + \\ \;\;\;\;\;\;\;\;{{\boldsymbol{E}}_{A(m)}}δ{\boldsymbol{\xi }} + {\boldsymbol{R}}{{\boldsymbol{e}}_{As}} + {{\boldsymbol{e}}_y}\\ {{\boldsymbol{G}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right) + {{{\boldsymbol{G'}}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right)\delta {\boldsymbol{\xi }} = {{\boldsymbol{b}}_j} \end{array} \right. $

(9)

式中，G′ _j(·)为G_j(·)的一阶导数。

按求条件极值的方法构造目标函数：

$ \varphi \left( {{{\boldsymbol{e}}_y}, \;{{\boldsymbol{e}}_{As}}} \right) = {\boldsymbol{e}}_y^{\mathop{\rm T}\nolimits} {\boldsymbol{Q}}_y^{ - 1}{{\boldsymbol{e}}_y} + {\boldsymbol{e}}_{As}^{\mathop{\rm T}\nolimits} {\boldsymbol{Q}}_{As}^{ - 1}{{\boldsymbol{e}}_{As}} +\\ \;\;\;\;\;\; 2{{\boldsymbol{K}}^{\rm{T}}}\left( {{{{\boldsymbol{G'}}}_j}\left( {{\xi _{(m)}}} \right)\delta {\boldsymbol{\xi }} + {{\boldsymbol{G}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right) - {{\boldsymbol{b}}_j}} \right) $

(10)

式中，K为对应约束条件的联系数向量(拉格朗日算子)。

对各变量求偏导并令导数为0，可得：

$ \begin{array}{l} \;\;{{{\boldsymbol{\hat \beta }}}_{(m + 1)}} = \left( {{\boldsymbol{N}}_{BB(m)}^{ - 1} - {\boldsymbol{N}}_{BB(m)}^{ - 1}{{{\boldsymbol{G'}}}_j}{{\left( {{{\boldsymbol{\xi }}_{(m)}}} \right)}^{\rm{T}}}} \right.\\ \;\;\;\;\left. {{\boldsymbol{N}}_{GG(m)}^{ - 1}{{{\boldsymbol{G'}}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right){\boldsymbol{N}}_{BB(m)}^{ - 1}} \right){{\boldsymbol{W}}_{(m)}} - \\ {\boldsymbol{N}}_{BB(m)}^{ - 1}{{{\boldsymbol{G'}}}_j}{\left( {{{\boldsymbol{\xi }}_{(m)}}} \right)^{\rm{T}}}{\boldsymbol{N}}_{GG(m)}^{ - 1}\left( {{{\boldsymbol{G}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right) - {{\boldsymbol{b}}_j}} \right) \end{array} $

(11)

$ {{\boldsymbol{v}}_{(m + 1)}} = {{\boldsymbol{B}}_{(m)}}{\hat \beta _{(m + 1)}} - {{\boldsymbol{l}}_{(m)}} $

(12)

式中，$ {{\boldsymbol{N}}_{GG(m)}} = {{\boldsymbol{G'}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right){\left( {{\boldsymbol{B}}_{(m)}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{B}}_{(m)}}} \right)^{ - 1}}{{\boldsymbol{G'}}_j}{\left( {{{\boldsymbol{\xi }}_{(m)}}} \right)^{\rm{T}}}$，${{\boldsymbol{N}}_{BB(m)}} = {\boldsymbol{B}}_{(m)}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{B}}_{(m)}} $，${{\boldsymbol{W}}_{(m)}} = {\boldsymbol{B}}_{(m)}^{\mathop{\rm T}\nolimits} {{\boldsymbol{Q}}^{ - 1}}{{\boldsymbol{l}}_{(m)}} $，B_(m)、Q、l_(m)含义同式(6)和式(7)。

拉格朗日算子及其对应的方差协方差阵为：

$ {\boldsymbol{K}} = {\boldsymbol{N}}_{GG(m)}^{ - 1}\left( {{{{\boldsymbol{G'}}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right){\boldsymbol{N}}_{BB(m)}^{ - 1}{{\boldsymbol{W}}_{(m)}} + {{\boldsymbol{G}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right) - {{\boldsymbol{b}}_j}} \right) $

(13)

$ {{\boldsymbol{D}}_{KK}} = \hat \sigma _0^2{{\boldsymbol{Q}}_{KK}} = \hat \sigma _0^2{\boldsymbol{N}}_{GG(m)}^{ - 1} $

(14)

式中，$ {\hat \sigma _0}$为单位权中误差。

需要说明的是，在实际坐标转换问题中，仿射基准变换模型允许不同坐标轴采用不同的尺度和旋转参数，这样更具一般性，可将其作为GRM；对于线性回归分析和自回归分析问题，可首先选择N/3~2N/3(N为样本观测值个数)之间的整数作为预设拟合最高阶数^[7]，若在后续分析中接近预设最高阶数仍未得到最佳模型，可再调整预设拟合最高阶数。

2 最佳线性回归模型选择

以众多待选模型中参数最多的模型作为GRM，则除GRM外的待选模型可表示为附有参数约束的线性回归模型，进而可以基于附有参数约束的测量平差理论和含有多个备选假设的假设检验理论，从中选出既符合观测数据实际，又具有最少模型参数的线性回归模型，称为初选回归模型(primary regression model, PRM)。通过分析不同模型平差后的参数精度，进一步筛选出最佳回归模型(optimum regression model, ORM)。

2.1 模型初选

要从众多待选模型中选出最佳模型，首先要对各模型中相应的附加参数约束进行多个备选假设检验和显著性分析，以定位出不兼容的参数约束。考虑到最佳模型应是参数解算精度高且模型参数尽可能少的平差模型^[17]，提出原假设和备选假设分别为：

$ \left\{ \begin{array}{l} {H_0}:{F_{(f + 1)}} = 0\;{\rm{is\;the\;primary\;regression\;model}}\\ {H_1}:{F_1} = 0\;{\rm{is\;the\;primary\;regression\;\;model}}\\ {H_2}:{F_2} = 0\;{\rm{is\;the\;primary\;regression\;\;model}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ {H_f}:{F_f} = 0\;{\rm{is\;the\;primary\;regression\;\;model}} \end{array} \right. $

(15)

式中，H₀为原假设；H₁，H₂，…，H_f为f个备选假设。当F_(f+1)=0中参数约束G_(f+1)ξ=b_(f+1)对回归模型影响不显著时接受H₀，则F_(f+1)=0为PRM；否则拒绝H₀，备选模型中PRM需通过进一步假设检验分析来确定。

式(10)中拉格朗日算子是为求解条件极值而引入的中间未知量，已有研究表明^{[11, 19]}，其数值大小反映的是参数约束对回归模型影响强度的大小，当约束条件对无约束线性回归模型影响显著时，拉格朗日算子往往较大，反之较小。基于此，构造检验统计量：

$ T = \mathop {\max }\limits_{s = 1, 2...c'} \left| {\frac{{{k_s}}}{{{{\hat \sigma }_{{k_s}{k_s}}}}}} \right| $

(16)

式中，c′为原假设H₀中参数约束条件的个数；k_s为采用参数最多的模型平差后得到的第s个约束条件对应的拉格朗日算子，$ {\hat \sigma _{{k_s}{k_s}}}$为其标准差。H₀为真时，T~t(n－t+c′－1)；反之，T~t′(n－t+c′－1, δ)，δ为非中心参数。

当参数约束与观测值不兼容时，求得的单位权中误差往往偏大。由式(16)可知，这可能会使检验统计量T普遍偏小，进而导致“纳伪”错误增加^[20]。为此，对于式(14)中的单位权中误差$ {\hat \sigma _0}$，可采用无约束线性回归模型求得：

$ {\hat \sigma _0} = \sqrt {{\boldsymbol{v}}_{{\rm{(}}m + {\rm{1)}}}^{\rm{T}}{\boldsymbol{Q}}{{\boldsymbol{v}}_{{\rm{(}}m + {\rm{1)}}}}/\left( {n - t} \right)} $

(17)

式中，v_(m+1)为无约束线性回归模型在满足收敛条件时的观测值残差向量，Q为其对应的协因数阵。

给定显著水平α，其表示H₀为真时拒绝H₀的概率，只有c′个约束条件全部成立时才接受H₀。设c′个约束条件之间相互独立，且每个约束条件成立的概率为1－α₁，则c′个约束条件全部成立的概率P₁为：

$ {P_1} \approx \prod\limits_{i = 1}^{c'} {\left( {1 - {\alpha _1}} \right)} = {\left( {1 - {\alpha _1}} \right)^{c'}} $

(18)

于是有：

$ \begin{array}{l} \;\;\;\;1 - \alpha = {\left( {1 - {\alpha _1}} \right)^{c'}} = 1 + c'\left( { - {\alpha _1}} \right) + \\ \left( {c'\left( {c' - 1} \right)/2} \right){\left( { - {\alpha _1}} \right)^2} + \left( {c'\left( {c' - 1} \right)\left( {c' - 2} \right)} \right./\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\left( {3 \times 2} \right)} \right){\left( { - {\alpha _1}} \right)^3} + \cdots \end{array} $

(19)

由于α₁为远小于1的常数，因此式(19)可忽略高次项，从而得到：

$ 1 - \alpha \approx 1 - c'{\alpha _1} $

(20)

$ {\alpha _1} \approx \alpha /c' $

(21)

因此对于统计检验量T，当T < t_α1/2时，接受H₀，PRM为F_(f+1)=0；否则拒绝H₀，PRM为众多备选模型其中一个，而要确定PRM，需重新进行假设检验。具体步骤为：

1) 将待选线性回归模型中参数最多的模型表示为式(4)形式，以参数最小二乘解为初值，根据式(6)和式(7)迭代求解无约束线性回归模型的系数矩阵误差E_A(m)、待估参数ξ_(m)和验后单位权中误差$ {\hat \sigma _0}$。

2) 将其他待选模型统一为式(8)形式，并在(ξ_(m)，e_As(m))处采用泰勒级数展开，得到式(9)形式。

3) 设原假设H₀和备选假设H_A。H₀：待选模型中参数个数最少者为PRM；H_A：其他模型为PRM。

4) 对于由参数约束最多的模型，根据式(13)和式(14)得到拉格朗日算子向量K及其对应的方差协方差阵D_KK。

5) 根据式(16)和式(21)构造假设检验统计量T并进行t检验，若T < t_α1/2，H₀成立，得到PRM，算法结束；否则执行步骤6)。

6) 从待选模型组合中删除原假设对应的线性回归模型，得到新的待选模型组合，并重复步骤3)~5)，直到H₀成立，算法结束。

2.2 最佳模型选择

要从众多待选模型中得到最佳模型，需要对假设检验得到的PRM与其他待选模型作进一步分析。为此，对比分析所有待选模型平差后的参数精度。当各模型参数对应的中误差接近时，取参数个数最少的模型为最佳模型；而当全部或部分模型参数对应的中误差相差较大时，取最小中误差对应的模型为最佳回归模型(ORM)。

当无约束线性回归模型取得最优解时，待估参数的方差协方差阵为：

$ {{\boldsymbol{D}}_{{{\hat \beta }_{{\rm{(}}m{\rm{)}}}}{{\hat \beta }_{{\rm{(}}m{\rm{)}}}}}} = \hat \sigma _0^2{{\boldsymbol{Q}}_{{{\hat \beta }_{{\rm{(}}m{\rm{)}}}}{{\hat \beta }_{{\rm{(}}m{\rm{)}}}}}} = \hat \sigma _0^2{\boldsymbol{N}}_{BB{\rm{(}}m{\rm{)}}}^{ - 1} $

(22)

式中，$ {\hat \sigma _0}$为采用无约束线性回归模型求得的验后单位权中误差。

当第j个附有参数约束的待选模型取得最优解时，待估参数的方差协方差阵为：

$ \begin{array}{l} {{\boldsymbol{D}}_{{{\hat \beta }_j}{{\hat \beta }_j}}} = \hat \sigma _j^2{{\boldsymbol{Q}}_{{{\hat \beta }_j}{{\hat \beta }_j}}} = \hat \sigma _j^2\left( {{\boldsymbol{N}}_{BB(m)}^{ - 1} - {\boldsymbol{N}}_{BB(m)}^{ - 1}{{{\boldsymbol{G'}}}_j}{{\left( {{{\boldsymbol{\xi }}_{(m)}}} \right)}^{\rm{T}}}} \right.\\ \left. {\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{N}}_{GG(m)}^{ - 1}{{{\boldsymbol{G'}}}_j}\left( {{{\boldsymbol{\xi }}_{(m)}}} \right){\boldsymbol{N}}_{BB(m)}^{ - 1}} \right) \end{array} $

(23)

式中，$ {\hat \sigma _j}$为采用第j个待选模型求得的验后单位权中误差，其求解公式为：

$ {\hat \sigma _j} = \sqrt {{\boldsymbol{v}}_{j{\rm{(}}m + 1{\rm{)}}}^{\rm{T}}{\boldsymbol{Q}}{{\boldsymbol{v}}_{j{\rm{(}}m + 1{\rm{)}}}}/\left( {n - t + {{c'}_j}} \right)} $

(24)

式中，v_j(m+1)为第j个待选模型取得最优解时求得的随机误差向量，Q为其对应的协因数阵，c′_j为第j个待选模型中参数约束条件的个数。

为叙述方便，称以上包含模型初选和最佳模型选择等环节的线性回归模型优选方法为引入拉格朗日算子的最佳线性回归模型选择算法(optimum linear regression model selection algorithm with Lagrange multipliers)，简称OLRS-LM法。OLRS-LM算法步骤为：首先将众多待选模型统一为附有参数约束的线性回归模型，然后进行最佳线性回归模型的筛选分析。由此可见，当不考虑回归模型中自变量的测量误差时，该方法同样可行，只是无需进行模型线性化及迭代求解。

3 算例分析

设计坐标转换模型优选实验对OLRS-LM算法的应用效果进行分析。选取某GPS网中19个控制点(点号为1~19)在工程独立坐标系下的坐标作为Ⅰ套坐标系下的模拟真值，并设计坐标转换参数：平移参数X₀=10 m、Y₀=20 m，尺度参数κ_x=1、κ_y=1，旋转角ω_x=5.00″、ω_y=5.00″，得到控制点在Ⅱ套坐标系下的坐标真值。控制点两套坐标真值见表 1。

在控制点的2套坐标中添加期望u=0、中误差σ=5.0 mm的正态分布随机误差，分别采用OLRS-LM算法和线性假设法进行坐标转换模型优选分析。由于经典的线性假设法只能用于传统线性回归模型，而且要求被检验的2个模型中至少有一个成立，因此，本文在应用时对经典的线性假设法进行适当改进：1)逐个选择各待选模型与参数最多的模型，进行单个备选假设检验；2)对检验通过的待选模型两两组合进行检验，由此得到最佳的平差模型。该过程可确保每次检验的2个模型中至少有一个符合观测数据实际，以保证线性假设法检验结果的可靠性。

在平面坐标转换模型中，尽管两参数和三参数模型实际应用较少，却是两种重要的模型，如在定量分析不同坐标系之间的系统误差时，应尽可能地考虑不同坐标系之间的差异可能性。因此，将待选转换模型均设为仿射变换模型、相似变换模型、三参数模型和两参数模型。其中，后3种转换模型均可由仿射变换模型附加相应参数约束得到，具体如下。

1) 仿射变换模型：

$ \left\{ \begin{array}{l} {x_{{\rm{I}}{{\rm{I}}_i}}} = {X_0} + {a_1}{x_{{{\rm{I}}_i}}} + {a_2}{y_{{{\rm{I}}_i}}}\\ {y_{{\rm{I}}{{\rm{I}}_i}}} = {Y_0} - {b_1}{x_{{{\rm{I}}_i}}} + {b_2}{y_{{{\rm{I}}_i}}} \end{array} \right. $

(25)

式中，(x_Ⅰi, y_Ⅰi)、(x_Ⅱi, y_Ⅱi)分别为第i个点在Ⅰ、Ⅱ两套坐标系中的二维坐标，a₁=κ_xcosω_x、a₂=κ_ysinω_y、b₁=κ_xsinω_x、b₂=κ_ycosω_y，(X₀，Y₀)、(κ_x，κ_y)、(ω_x，ω_y)分别对应两个坐标轴的平移、尺度和旋转参数。

坐标变换过程中图形保持正形特点需满足柯西-黎曼微分方程：

$ \left\{ \begin{array}{l} \frac{{\partial \Delta {x_i}}}{{\partial {x_{{{\rm{I}}_i}}}}} = \frac{{\partial \Delta {y_i}}}{{\partial {y_{{{\rm{I}}_i}}}}}\\ \frac{{\partial \Delta {x_i}}}{{\partial {y_{{{\rm{I}}_i}}}}} = - \frac{{\partial \Delta {y_i}}}{{\partial {x_{{{\rm{I}}_i}}}}} \end{array} \right. $

(26)

式中，Δx_i=x_Ⅱi－x_Ⅰi，Δy_i=y_Ⅱi－y_Ⅰi。

由式(26)可得满足正形变换条件的相似变换模型，附加参数约束为^[1]：

$ \left\{ \begin{array}{l} {a_1} - {b_2} = 0\\ {a_2} - {b_1} = 0 \end{array} \right. $

(27)

在相似变换模型基础上，若要得到尺度参数为1的三参数模型，附加参数约束为：

$ a_1^2 + a_2^2 = 1 $

(28)

若2套坐标系之间仅有沿坐标轴的平移，则为两参数模型，附加参数约束为：

$ {a_2} = 0 $

(29)

根据式(25)~(29)可得顾及自变量和因变量测量误差的仿射变换模型，将其他待选模型表示为附有参数约束的线性回归模型。OLRS-LM算法的模型参数优选分析结果见表 2和表 3，线性假设法结果见表 4。其中，2种方法的显著水平均分别取α₁=0.05和α₂=0.01。

从表 2和表 3可以看出，三参数模型为初选平差模型，且为最佳平差模型。具体来看，无论取显著水平α₁=0.05还是α₂=0.01，对于两参数模型，根据拉格朗日算子构造的假设检验统计量远大于阈值，表明两参数模型中附加的部分(或全部)参数约束与观测数据之间存在不兼容，同时也说明旋转参数在坐标转换模型中具有显著作用，必须引入；对于三参数模型，检验统计量小于阈值，说明该模型中附加的全部参数约束与观测数据之间兼容性较好，旋转角在X、Y两个方向的差异不显著。分别采用三参数模型、相似变换模型和仿射变换模型进行平差，两个平移参数的中误差存在明显差异，其中，三参数模型和相似变换模型求得的模型参数精度明显高于仿射变换模型，且三参数模型的解算精度高于相似变换模型。由此可以说明，三参数模型更符合观测数据的实际情况，为最佳平差模型，这与仿真数据的实际情况相符。

由表 4可知，两参数模型、三参数模型和相似变换模型分别与仿射变换模型组合进行差异性检验，无论取显著水平α₁=0.05或α₂=0.01，两参数模型对应的检验统计量均远超过阈值，相似变换模型对应的检验统计量均小于阈值；三参数模型对应的检验统计量大于阈值F_α1，小于阈值F_α2。由此说明，Ⅰ、Ⅱ两坐标系之间无显著的尺度差异，但可能存在旋转角。需要注意的是，此时并无法确定三参数模型即为参数解算精度最高的最佳模型，实际上其仅仅是完成最佳平差模型的初选。

分别采用4种待选坐标转换模型计算坐标转换参数，并用模型参数的均方误差M^[21]评定参数估计结果的准确度(参数均方误差越小，参数估值准确度越高)，以进一步验证以上分析结果的正确性，结果见表 5。

由表 5可知，两参数模型求得参数的均方误差较大，且解算得到的坐标转换参数估计值严重偏离设计真值；而另外3种模型的参数估计结果较好，参数估值与设计值均相差较小，所得参数的均方误差也都较小。三参数模型、相似变换模型和仿射变换模型求得的坐标转换参数虽然均与设计值相近，但三参数模型的均方误差最小，而且求得的平移、旋转和尺度参数均与设计值最接近。由此可见，三参数模型较相似变换模型和仿射变换模型准确度更高，结果更优。为说明忽略自变量误差建立回归模型的不足，进一步建立不考虑自变量误差的三参数坐标转换模型，并采用最小二乘法求解坐标转换参数，计算其模型参数均方误差M，结果见表 6。

结合表 5和表 6可知，忽略自变量误差后，所求参数的均方误差变大，解算的坐标转换参数除平移参数Y₀基本未变化外，其他转换参数均与设计值存在较大差异。

综上分析可知，只有采用合理的平差模型，才能准确表达观测量之间的物理或几何关系。在实际应用中，应首先确定最符合观测数据实际的平差模型，然后再进行求解，否则可能得到错误的结果。同时可以发现，采用OLRS-LM算法和改进的线性假设法均能对顾及自变量和因变量误差的线性回归模型选择作出客观、量化的评判，且OLRS-LM算法更为简便。此外，在相同显著水平下，采用拉格朗日算子构造假设检验统计量可得到正确的判断结果，而线性假设法会出现“弃真”错误，表明以拉格朗日算子构造假设检验统计量较线性假设法的检验功效更高，得到的结果更可靠。

4 结语

对线性回归分析、坐标转换和自回归分析问题进行建模时，应首先进行模型优选分析，以建立既符合观测数据实际、参数解算精度又高的最佳模型。本文基于附有参数约束的测量平差理论和含有多个备选假设的假设检验理论，以拉格朗日算子构造假设检验统计量，提出最佳线性回归模型选择算法，推导其求解公式，并设计具体求解算法，通过实验验证算法的有效性，得到以下结论：

1) 最佳模型的形式完全由观测数据的实际情况决定，将众多待选模型统一为附有参数约束的线性回归模型，该观点是线性回归分析、坐标转换和自回归分析等问题建模时获得最佳模型的依据，由此可将模型的优选问题转化为含有多个备选假设的假设检验问题。

2) 对于合理的参数约束，顾及参数约束的线性回归模型，其参数解算精度较无约束的线性回归模型可得到一定程度提高，因此在实际建模时，对模型进行优选分析具有重要意义。

3) OLRS-LM算法可准确找出既符合观测数据实际、参数解算精度又高的最佳平差模型，其以拉格朗日算子构造假设检验统计量，能够客观、量化地诊断参数约束与观测数据之间的兼容性，较以残差平方和构造假设检验统计量的线性假设法的检验功效更高，结果更可靠。