王建民，倪福澤，趙建軍

（1. 太原理工大學礦業(yè)工程學院，山西太原030024；2. 成都理工大學地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室，四川成都610059；3. 中煤（西安）航測遙感研究院有限公司，陜西西安710100）

高斯?馬爾可夫模型在許多工程實踐中得到成功應用，通常認為模型中的設計矩陣是無誤差的。然而，在許多情景下，設計矩陣中的元素并不能保證是常數(shù)，其本身也是含有隨機誤差的觀測值，設計矩陣和觀測向量中含有隨機誤差的模型通常被定義為變量誤差模型（errors-in-variables，EIV）［1-2］。

最初，因為求解EIV模型參數(shù)的方法不成熟，通常忽略設計矩陣中的隨機誤差，而采用加權最小二乘法（weighted least-squares，WLS）近似估計EIV參數(shù)，這種方法是不嚴密的。1980 年，Golub 等［3］提出采用整體最小二乘法（total least-squares，TLS）估計EIV模型的未知參數(shù)。當設計矩陣和觀測向量均含有隨機誤差時，TLS 就是最小二乘（least-squares，LS）的擴展且更具一般性［4］。TLS 將設計矩陣中的所有元素按等權對待，但實際上設計矩陣中的元素并非都含有隨機誤差，且精度也不盡相同，因此TLS被進一步擴展為加權整體最小二乘法（weightedTLS，WTLS），WTLS 分別給觀測向量和設計矩陣附加了不同的權重，比TLS更具一般性。

一般來說，主要有兩類方法實現(xiàn)TLS和WTLS的求解，一類方法是基于奇異值分解（singular value decomposition，SVD）的 直 接 解 算 方 法［4–8］，Akyilmaz［6］在SVD分解中給設計矩陣和觀測向量賦予了相應的權矩陣，但引入的權陣并沒有真正起到調(diào)配誤差的作用，其本質(zhì)上是對每個參數(shù)和觀測值進行了縮放［9］。另一類常用的方法是基于拉格朗日輔助法構建極值函數(shù)分別對未知量進行求偏導數(shù)，許多文獻據(jù)此提出了WTLS的迭代算法［2，5，8，10-12］，建立在統(tǒng)計框架下的這些迭代算法估計EIV模型參數(shù)是嚴密的［13］。因為WTLS 同時顧及了設計矩陣和觀測向量的權重，在工程實踐中應用迭代法估計EIV模型參數(shù)更為普遍。

從現(xiàn)有的研究成果來看，雖然國內(nèi)外學者就EIV 模型提出了許多卓有成效的解算方法，但是與傳統(tǒng)的WLS 相比，算法比較復雜，不具有操作上的優(yōu)勢［14］，而且大多數(shù)迭代算法僅考慮了參數(shù)估計精度和迭代收斂速率，當面臨大數(shù)據(jù)集時計算效率有限。

為了提高WTLS 的效率，一些學者從不同側面改進了WTLS 算法。Shen 等［10］提出基于Newton-Gauss 的非線性WLS 平差方法，Jazaeri 等［15］導出了基于WLS準則下的迭代算法，這些方法都明顯提升了迭代收斂速率。但僅通過增強迭代收斂速率，計算效率提升有限。事實上，絕大多數(shù)EIV 模型的設計矩陣具有結構性特征，僅有部分元素是隨機的，其通常的解法是構造特定的權矩陣，使系數(shù)矩陣中固定元素位置的改正數(shù)為零，但這樣會導致部分矩陣維數(shù)增大，從而加大了矩陣的運算量，不利于計算效率的提高。為此，Xu 等［9］提出了partial EIV（PEIV）模型，并給出了相應的算法。PEIV模型一個顯著的優(yōu)點是未知數(shù)的總量明顯減小，然而多次迭代收斂拖慢了其計算效率［9，16-17］。Zhao［17］應用拉格朗日輔助法對參量求導增強了PEIV 算法的效率，王樂洋等［16］應用間接平差的原理進行改進，使算法更加簡潔、緊湊。這兩種方法改進的結果是迭代快速收斂，效率明顯提高，但不足之處是需要根據(jù)原設計矩陣的結構解析出一個固定矩陣和固定向量，而且在迭代過程仍需重構設計矩陣，不利于效率的提升。

WTLS在大地測量等眾多科學研究和工程領域中得到應用［18-19］，WTLS算法基本成熟，漸成理論體系。算法之間的差異主要體現(xiàn)在計算效率和算法原理的復雜程度上。在大數(shù)據(jù)的背景下，計算效率是決定算法能有效應用的一個重要因素。例如應用WTLS處理點云數(shù)據(jù)，數(shù)據(jù)量會顯著增加，算法性能就會顯得尤為重要。

本文針對EIV模型設計矩陣呈現(xiàn)出的結構性特征，將設計矩陣分成隨機和固定部分，僅給隨機部分賦予權重，并且不涉及拉格朗日輔助法對參量求偏導，在WLS 平差約束準則下，推證了一種估計EIV模型參數(shù)的迭代算法，該算法能夠顯著降低矩陣的運算量，極大提升計算效率。

1 基于WLS 準則的WTLS 高效算法原理

1.1 EIV模型及WTLS算法

EIV模型的函數(shù)模型和隨機模型通常表示為［8-9］

式（1）、（2）中：E為（n×t）階設計矩陣A的隨機誤差；e為n維觀測向量L的隨機誤差；（n×n）階矩陣Qe和（nt×nt）階矩陣QE分別為e 和vec（E）的協(xié)因數(shù)矩陣；vec（?）表示矩陣的拉直變換；θ為待估計的t個未知參數(shù)；σ2為單位權方差。

對于EIV 模型，可直接應用最小二乘準則eTQe-1e=min，其參數(shù)估值為［15］

僅對隨機誤差e 加以約束，難以估計E，一般的WTLS的平差準則為［10-11，15］

在式（4）的約束條件下，利用拉格朗日輔助法可以精確估計未知參數(shù)θ［8］，即

式（5）、（6）是WTLS 常用的迭代算式，一般的WTLS 方法在迭代過程中需要更新矩陣E?，不利于算法效率的提升。接下來，根據(jù)設計矩陣A 呈現(xiàn)出的結構性特征，推導一種基于WLS準則的EIV模型參數(shù)估計高效算法。

1.2 顧及設計矩陣結構的WTLS算法

理論上EIV模型中的設計矩陣的元素均可包含隨機誤差，但在實際應用中，絕大多數(shù)EIV模型的設計矩陣都會存在部分固定列并且不含有隨機誤差［9，17，20］。對于含有固定列的設計矩陣，通常是構建一個特定的協(xié)因數(shù)矩陣QE，在迭代的過程中使得固定元素的改正數(shù)為零，但這樣會導致QE的維數(shù)被固定列“擴大”，制約了計算效率。

設A中有t1列固定列和t2列隨機列，將A中的固定列和隨機列進行矩陣分塊，對應的參數(shù)θ 也分成兩部分，即

式（7）、（8）中：A1為由固定列組成的（n×t1）階矩陣；A2為由隨機列組成的（n×t2）階矩陣；與A1和A2對應的參數(shù)分別為θ1和θ2。

根據(jù)式（7）和式（8），EIV 模型（式（1）和式（2））重新改寫為

式（9）、（10）中：E2為A2的隨機誤差矩陣，其階數(shù)為（n×t2）；QE2為vec（E2）的協(xié)因數(shù)矩陣，其階數(shù)為（nt2×nt2）。

將式（9）改寫為以下誤差方程：

顧及式（11），vec(A2)表示為

將式（11）和式（12）組合形成如下的非線性方程，即

為進一步化簡式（14），需利用式（15）的矩陣求逆引理［21］，即

由式（16）可知，設計矩陣A2的隨機誤差的估值為

由式（7）可知，觀測向量L的殘差可表示為

由式（18）可知

將式（19）代入式（3），得

則

式中：λ=Q-21(L-Aθ?)，λ為n×1的向量。由式（21）可知

根據(jù)A的矩陣分塊（式（7）），式（22）可分解為兩部分，即

上述推導是基于WLS 最小二乘準則，并且式（32）在形式上與間接平差的法方程類似，但方程兩側都存在未知參數(shù)θ?，仍需迭代計算。為了加以區(qū)別和簡寫，上述算法稱為部分加權整體最小二乘法（partially weighted TLS，PWTLS），PWTLS 的算法流程步驟如下：

由上述算法可知，在迭代過程中，并不涉及對設計矩陣A拆分成A1和A2，僅需把隨機列A2對應的參數(shù)θ2分解出來參與迭代。只要EIV模型的設計矩陣存在固定列，根據(jù)列與參數(shù)的對應關系，很容易將A組合成固定和隨機部分。通常情況下，習慣上總是將EIV 模型的固定列排列在一起，而且絕大多數(shù)EIV 模型至少有一列是固定列，因而PWTLS 估計EIV參數(shù)具有普適性。另外，在實際應用中，設計矩陣的隨機列可能會含有少量固定元素，例如三維七參數(shù)坐標轉換模型中，有3列固定列組成A1，在隨機列中有一個固定零元素，對于隨機列中的個別固定元素賦予值為零的協(xié)因數(shù)后，仍可組成隨機列A2。

從算法流程可以看出，PWTLS 算法不同于一般的WTLS 算法，在每次迭代的過程中并不需要估計隨機誤差E；通過給設計矩陣的隨機列賦予權重，使得相關矩陣運算的時間復雜度由O（nt）降低到O（nt2），占用的內(nèi)存空間由O（nt2）減小到O（nt22），例如原參數(shù)矩陣X?的維數(shù)（nt × nt）降為X?2（nt2×nt2），基于這兩點，使PWTLS算法的效率顯著提高。

2 算例與分析

為了驗證PWTLS的性能，本文選擇了7種具有代表性的WTLS算法［10-11，15-17，20，22］與本文算法進行對比分析，這7 種算法從不同角度對WTLS 進行效率增強，具體原理見參考文獻，有的算法在相關文獻中也與其他算法進行過對比研究。

2.1 真實算例

直線擬合簡單而有效，經(jīng)常用于證明WTLS 算法的有效性。設有直線方程為y=ax+b，其中x 和y 是直線的坐標分量，a 和b 是直線的斜率和截距。實驗數(shù)據(jù)來自于文獻［23］，見表1，其中（xi，yi）和（wxi，wyi）表示坐標和相應的權值。

根據(jù)直線方程和坐標點數(shù)量，設計矩陣A，參數(shù)向量θ及觀測向量L具體形式為

表1 坐標觀測值和相應的權重Tab.1 Coordinate observations and weights

在式（35）中，設計矩陣A 可分成A1和A2兩部分，即

A2的協(xié)因數(shù)矩陣QE2的結構取自QE的一部分，即

設迭代收斂閾值ε=10-10，分別應用8種算法估計參數(shù)和單位權中誤差，實驗結果見表2。從表2中可明顯看出，8 種嚴密算法估計的參數(shù)相差在10-11量級且中誤差相同，表明PWTLS能夠獲得和另外7種嚴密算法相同的解。文獻［22］對PEIV 模型進行改進，其目的也是提升效率，從本例的結果來看，收斂速度高達454 次，其余算法的收斂速度都較快。文獻［16］同樣是對PEIV模型算法進行優(yōu)化，取得了理想的效果。比較而言，文獻［20］的算法收斂最快，收斂速度僅是影響計算效率一方面的因素。本例中直線擬合的設計矩陣結構簡單而且僅有10個點，主要用于PWTLS參數(shù)估計精度的驗證及收斂速度的初步分析，接下來通過大數(shù)據(jù)量應用仿真模擬實驗來對比分析各個算法的效率。

表2 8種不同算法的參數(shù)和中誤差的估計結果Tab.2 Estimated parameters and standard deviations obtained using eight algorithms

2.2 仿真算例

平面相似四參數(shù)轉換模型是常用的坐標基準轉換模型之一，通常由式（38）表示［24］

式中：ξ 和η 為平移參數(shù)；k為尺度參數(shù)；a為旋轉角；下標S和T分別表示原坐標系和目標坐標系。

假設有d 個重合點用于求解轉換參數(shù)，式（38）經(jīng)過變換可由EIV 模型描述，設計矩陣A，參數(shù)θ 及觀測向量L表示為

式中：u=k cos a；w=k sin a為輔助參數(shù)。

參數(shù)ξ 和η 對應的系數(shù)是1 和0 的常量，設計矩陣A 的第1 列和第2 列是固定列，組成A1；其他兩列組成A2形成隨機列，A1和A2表示為

對應于A1和A2，相應的參數(shù)可分成兩部分，即

原坐標系中的坐標由隨機函數(shù)產(chǎn)生d 個坐標點，目標坐標由給定的參數(shù)真值ξ=-27.366，η=-71.185，u=1.000 001 092 和w=6.400 15×10-7通過轉換模型反算得到。原坐標和目標坐標中的隨機誤差由隨機函數(shù)（式（42）、（43））隨機生成，并將隨機誤差附加到相應的坐標分量中。

式中：ΔS和ΔT分別表示原始坐標和目標坐標中的真誤差，mvnrnd（0，σ，2d）生成期望為0、標準差為σ=0.05 m、維數(shù)為2d×1的隨機向量。

由隨機誤差和標準差定義觀測量向量L的協(xié)因數(shù)陣QeT

根據(jù)設計矩陣A2的結構，vec（A2）的協(xié)因數(shù)陣QE2定義為

基于上述的模型和數(shù)據(jù)，隨機模擬1 000 次實驗，固定重合點數(shù)d=200，應用均方根誤差（RMSE）來衡量參數(shù)估計的精度，參數(shù)ξ、η、u 和w 的RMSE由式（46）、（47）估計

式（46）、（47）中：ξ?i、η?i、u?i和w?i為經(jīng)第i次模擬的對應參數(shù)的估計值。

在相同數(shù)據(jù)和計算平臺的實驗條件下（收斂閾值ε=10-10），分別應用7 種具有較高收斂速度的算法進行模擬，參數(shù)的RMSE的統(tǒng)計結果和1 000次模擬的累計耗時見表3，結果表明，嚴密算法中的RMSE差異為10×10-5mm，進一步證明了PWTLS與其他嚴密算法的參數(shù)估計精度相同。

從計算的總耗時可以看出，PWTLS 是所有嚴密算法中計算效率最高的，與其他算法計算效率相比，最少提高了161%，最多提高了近416%。其主要是因為PWTLS使用QE2代替了QE，維數(shù)由800×800 降到400×400；另外，在每次迭代過程中，PWTLS算法并不需要估計隨機誤差矩陣E。

為了突出PWTLS 計算效率的優(yōu)越性，從另一個角度執(zhí)行了10次模擬實驗，要求重合點數(shù)據(jù)以間隔為100個點，逐漸從d=100增加到1 000。計算耗時統(tǒng)計結果如圖1 所示，結果表明，PWTLS 的耗時曲線變化相對平緩，意味著隨著數(shù)據(jù)量的增大，PWTLS 的優(yōu)勢越來越明顯。

表3 1 000次模擬的RMSE和累計耗時Tab.3 Results of parameter estimation and details of the time consumed for different algorithms

圖1 7種算法的計算效率比較Fig.1 Computational efficiency of seven algorithms

3 結語

WTLS 已經(jīng)在眾多科學和工程技術中得到應用，當面臨處理大數(shù)據(jù)集時，其計算效率有限，降低算法的復雜度和提高計算效率對WTLS 算法將顯得越來越重要。

本文針對EIV模型的設計矩陣呈現(xiàn)出的結構性特征，提出了基于WLS 準則約束條件下的PWTLS算法，無需引入輔助量，迭代方程與經(jīng)典的間接平差法方程類似，便于理解，易于實現(xiàn)，估計的參數(shù)是嚴密的。PWTLS 通過給設計矩陣中的隨機列賦予權重，使得設計矩陣的協(xié)因數(shù)陣維數(shù)降低，從而在迭代過程中縮減了所有參與運算矩陣的維數(shù)，極大地減少了計算量，促進了效率的提升；另外，PWTLS 算法避免了重復估計設計矩陣隨機列的誤差，這兩方面的因素使矩陣的運算量降低，從而顯著提升了計算效率。最后，通過算例驗證了所提出的PWTLS算法與同類算法參數(shù)估計精度相當，但計算效率顯著提高。

作者貢獻說明：

王建民：提出了論文中PWTLS算法原理，完成了論文的整體寫作。

倪福澤：編寫了計算程序，完成了論文中的實驗計算和分析，繪制了圖表。

趙建軍：提供研究經(jīng)費的資助，修改了中英文摘要。