王建民，倪福澤，趙建軍

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

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

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

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

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

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

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

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

1 基于WLS 準(zhǔn)則的WTLS 高效算法原理

1.1 EIV模型及WTLS算法

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

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

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

僅對隨機(jī)誤差e 加以約束，難以估計(jì)E，一般的WTLS的平差準(zhǔn)則為［10-11，15］

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

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

1.2 顧及設(shè)計(jì)矩陣結(jié)構(gòu)的WTLS算法

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

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

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

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

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

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

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

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

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

由式（16）可知，設(shè)計(jì)矩陣A2的隨機(jī)誤差的估值為

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

由式（18）可知

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

則

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

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

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

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

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

2 算例與分析

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

2.1 真實(shí)算例

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

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

表1 坐標(biāo)觀測值和相應(yīng)的權(quán)重Tab.1 Coordinate observations and weights

在式（35）中，設(shè)計(jì)矩陣A 可分成A1和A2兩部分，即

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

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

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

2.2 仿真算例

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

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

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

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

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

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

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

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

由隨機(jī)誤差和標(biāo)準(zhǔn)差定義觀測量向量L的協(xié)因數(shù)陣QeT

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

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

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

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

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

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

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

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

3 結(jié)語

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

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

作者貢獻(xiàn)說明：

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

倪福澤：編寫了計(jì)算程序，完成了論文中的實(shí)驗(yàn)計(jì)算和分析，繪制了圖表。

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