摘要 :針對(duì)直線間平行、垂直和斜交的3種位置關(guān)系,提出將直線間的位置關(guān)系轉(zhuǎn)化為等式約束條件,結(jié)合加權(quán)整體最小二乘法進(jìn)行直線擬合。首先,將直線間的3種位置關(guān)系表達(dá)為直線斜率的關(guān)系,對(duì)其線性化,得到3種約束模型;其次,根據(jù)加權(quán)整體最小二乘準(zhǔn)則確定附有等式約束的加權(quán)整體最小二乘法平差模型,推導(dǎo)其解算方程,并建立平差算法;最后,針對(duì)3種約束模型進(jìn)行擬合實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明:新算法具有可行性,擬合精度較高,且能確保直線間的平面位置關(guān)系正確。
關(guān)鍵詞 :平面位置關(guān)系;直線擬合;加權(quán)整體最小二乘;等式約束
中圖分類號(hào):P207"" 文獻(xiàn)標(biāo)志碼:A"" 文章編號(hào):1004-0366(2025)01-0001-08
直線擬合在工程實(shí)踐中有著廣泛的應(yīng)用[1-2]。傳統(tǒng)上,常用最小二乘法(LS,least squares)估計(jì)直線參數(shù),該方法易于理解與計(jì)算,但只考慮了因變量誤差,沒有考慮自變量誤差,導(dǎo)致估計(jì)結(jié)果只是某個(gè)坐標(biāo)軸上的最佳[3-4]。因此LS直線擬合結(jié)果會(huì)隨自變量選取的不同而不同,不是全體觀測(cè)誤差下的最優(yōu)估計(jì)結(jié)果。整體最小二乘法(TLS,total least squares)同時(shí)考慮自變量和因變量誤差,參數(shù)估計(jì)理論更嚴(yán)密,能夠得到全體觀測(cè)誤差平方和最小條件下的最優(yōu)估計(jì)結(jié)果。學(xué)者們針對(duì)TLS進(jìn)行了大量的研究,先后提出了奇異值分解算法[5-6](SVD,singular value decomposition)、TLS迭代算法和遞推的等價(jià)算法[7-11]。丁克良等[12-13]分別采用LS和TLS進(jìn)行了直線擬合,并對(duì)比分析了LS與TLS擬合精度,指出TLS能有效提高擬合精度,獲得了更好的擬合結(jié)果。針對(duì)非等精度的直線擬合問題,文獻(xiàn)[14-16]中提出了嚴(yán)密的加權(quán)整體最小二乘(WTLS,weighted total least square)直線擬合模型;王繼剛等[17]從測(cè)量平差角度入手,基于附有參數(shù)的條件平差理論,建立了一種簡(jiǎn)單的WTLS直線擬合模型。
當(dāng)先驗(yàn)約束條件已知時(shí),國內(nèi)外學(xué)者研究了附約束條件的WTLS平差問題[18-19],并導(dǎo)出了解算的通用公式[20]。先驗(yàn)約束條件有線性和非線性之分,王樂洋[21]對(duì)約束條件為線性的WTLS進(jìn)行了研究,胡川等[22]對(duì)非線性約束的WTLS進(jìn)行了研究。現(xiàn)有的研究多數(shù)是針對(duì)單一直線擬合問題進(jìn)行,然而在工程實(shí)踐過程中存在擬合直線與相鄰直線具有位置約束關(guān)系的情況,此時(shí),采用線性或非線性等式約束的加權(quán)整體最小二乘法(ECWTLS,equality constrained weight total least squares)來估計(jì)直線參數(shù)更為合理。本次研究以ECWTLS為基礎(chǔ),將平面直線的相對(duì)位置關(guān)系轉(zhuǎn)換成直線間斜率的約束關(guān)系,建立對(duì)應(yīng)的等式約束方程,構(gòu)建位置關(guān)系約束下的直線擬合模型。
1 位置關(guān)系約束下的直線參數(shù)整體最小擬合模型
1.1 顧及直線位置關(guān)系的直線擬合模型
常規(guī)WTLS的函數(shù)模型和隨機(jī)模型[6,11]為
y-ey=(A-EA)ξ, (1)
eyeA:=ey vec (EA)~00,σ20Qy00Q0Qx,(2)
其中:y和ey分別是n×1的觀測(cè)向量和誤差向量;A和EA分別是n×m的系數(shù)矩陣和誤差矩陣,且有 rank (A)=mlt;n;ξ為m×1的參數(shù)向量; vec 表示對(duì)矩陣的拉直運(yùn)算;
σ20表示未知的單位權(quán)方差;Qy和Qx都是n×n的協(xié)因數(shù)陣;Q0與A的結(jié)構(gòu)有關(guān);符號(hào)表示克羅內(nèi)克積。
對(duì)式(1)稍作變換,整理后[22]可得
y-ey=A(x) -(ξTIn)Nex, (3)
其中:A(x)表示A是關(guān)于x的矩陣函數(shù);In表示n階單位矩陣;ex表示n×1的誤差向量;N表示nm×n的轉(zhuǎn)換矩陣[23],根據(jù)向量求導(dǎo)法得到N的計(jì)算公式[24]為
N= d { vec (A(x))} d xT。 (4)
為將直線間的平面位置關(guān)系融入WTLS,首先將直線位置關(guān)系轉(zhuǎn)化為可以執(zhí)行的數(shù)學(xué)表達(dá)式。以兩條直線為例,根據(jù)兩直線間的位置關(guān)系,得到的數(shù)學(xué)表達(dá)式為
f(a1,b1,a2,b2)=0, (5)
其中:a1和a2分別表示兩條直線的斜率;b1和b2分別表示兩條直線的截距。
當(dāng)兩條直線為平行關(guān)系時(shí),兩直線的斜率相等,此時(shí)平行關(guān)系約束模型表示為
a1-a2=0, (6)
并令
a1=1+Δa1, (7)
a2=2+Δa2, (8)
其中:1和2分別表示兩條直線斜率的預(yù)測(cè)值;Δa1和Δa2分別表示兩條直線斜率的誤差預(yù)測(cè)值。
將式(7)和式(8)代入式(6)可得
Δa1-Δa2+1-2=0。 (9)
當(dāng)兩條直線相交但不垂直時(shí),即為斜交關(guān)系,兩條直線間存在一個(gè)夾角θ,此時(shí)位置關(guān)系約束模型表示為
a1-a21+a1a2= tan "θ, (10)
其中:以a1為斜率的直線傾斜角大于以a2為斜率的直線傾斜角。將式(7)和式(8)代入式(10)中并進(jìn)行線性化,斜交關(guān)系約束模型可進(jìn)一步表示為
Δa1-Δa2-2 tan "θΔa1-1 tan "θΔa2+1-2-(1+12) tan "θ=0。(11)
相交關(guān)系中存在一種特殊關(guān)系,即兩直線垂直,兩直線斜率之積等于-1,此時(shí)正交位置關(guān)系約束模型表示為
a1a2=-1。 (12)
將式(7)和式(8)代入式(12)中并進(jìn)行線性化,正交關(guān)系約束模型可進(jìn)一步表示為
2Δa1+1Δa2+12+1=0。 (13)
然后將式(9)、式(11)和式(13)表示為相同形式的約束模型:
Kδ=C, (14)
其中:K為等式約束矩陣;C為約束向量;δ為直線參數(shù)的誤差預(yù)測(cè)值向量,具體表示為(Δa1" Δb1" Δa2" Δb2)T。將式(14)作為約束方程加入 WTLS中組成ECWTLS 平差的函數(shù)模型。3種直線位置關(guān)系K和C的表達(dá)式如表1所列。
1.2 ECWTLS平差模型的解算
將TLS看成是一個(gè)非線性的參數(shù)估計(jì)問題[14],令
=A-A, (15)
=ξ(0)+δ。 (16)
將非線性方程(3)線性化為
Y-ey=δ-ZTx, (17)
其中:Y=y-Aξ(0),Z=NT(ξ(0)In)。
另外, ECWTLS 的極小值條件式為
eTyQ-1yey+eTxQ-1xex= min 。 (18)
采用拉格朗日乘數(shù)法求解 ECWTLS, 根據(jù)線性化的方程將目標(biāo)函數(shù)定義為
Φ(ey,ex,δ,λ,μ)=eTyQ-1yey+eTxQ-1xex+2λT(Y-ey-δ+ZTx)+2μT(Kδ-C),(19)
其中:λ(n×1)和μ(l×m)為拉格朗日乘積因子。對(duì)上式各待解參數(shù)求偏導(dǎo)并令它們等于0,可得
12Φey=Q-1yy-=0, (20)
12Φex=Q-1xx-Z=0, (21)
12Φ=-AT+KT=0, (22)
12Φλ=Y-y-δ+ZTx=0, (23)
12Φμ=Kδ-C=0。 (24)
根據(jù)式(20)和式(21)可得
y=Qy, (25)
x=-QxZ, (26)
將式(25)和式(26)帶入式(23)可得
=(Qy+ZTQxZ)-1(Y-Aδ)=R1(Y-Aδ),
(27)
其中:R1=(Qy+ZTQxZ)-1。再將式(27)帶入式(25)和式(26)可得
y=QyR1(Y-Aδ), (28)
x=-QxZR1(Y-Aδ)。 (29)
將式(27)代入式(22)可得
δ=(ATR1A)-1(ATR1Y-KT)。 (30)
令R2=(ATR1A)-1,再將式(30)代入式(24)可得
=(KR2KT)-1(KR2ATR1Y-C)。 (31)
對(duì)應(yīng)的單位權(quán)方差為
20=eTyQ-1yey+eTxQ-1xex(n-m+l)。 (32)
2 ECWTLS算法
根據(jù)以上推導(dǎo)的計(jì)算公式,顧及平面位置關(guān)系的直線參數(shù)ECWTLS估計(jì)算法設(shè)計(jì)如下:
輸入觀測(cè)向量 y1和y2、系數(shù)矩陣A1和A2及對(duì)應(yīng)的協(xié)因數(shù)陣Qx1、Qy1、Qx2和Qy2、轉(zhuǎn)換矩陣N。
(1) 計(jì)算矩陣:
y=y1y2,
A=A100A2,
Qx=Qx100Qx2,
Qy=Qy100Qy2。
(2) 計(jì)算ξ的初值:(0)=(AQ-1yA)-1AQ-1yy,并令E(0)A=0。
(3) 依照表1,根據(jù)直線間的位置關(guān)系計(jì)算相應(yīng)的約束矩陣K和約束向量C。
(4)根據(jù)線性化公式依次計(jì)算下列各值:
(i)=A-(i)A,
Y(i)=y-A(i-1),
Z(i)=NT((i-1)In)。
(5) 按下列公式計(jì)算參數(shù)值:
R(i)1=(Qy+Z(i)TQxZ(i))-1,
R(i)2=((i)TR(i)1(i))-1,
(i)=(KR(i)2KT)-1(KR(i)2(i)TR(i)1Y(i)-C),
δ(i)=R(i)2(TR(i)1Y-KT(i)),
(i)y=QyR(i)1(Y(i)-A(i)δ(i)),
(i)x=-QxZ(i)R(i)1(Y(i)-A(i)δ(i)),
(i)=(i-1)+δ(i)。
(6) 如果δ(i)gt;10-10,則E(i+1)A= vec -1(Ne(i)x),將值代回至第(3)步繼續(xù)進(jìn)行迭代計(jì)算;反之,終止迭代,并以式(32)進(jìn)行精度評(píng)定計(jì)算結(jié)束,以迭代終止前一步計(jì)算的值作為最終結(jié)果輸出,包括參數(shù)估計(jì)值和精度指標(biāo)σ20。
3 算例與分析
為了驗(yàn)證本文提出的直線間位置關(guān)系的約束模型,現(xiàn)針對(duì)直線間的平行、斜交和垂直3種位置關(guān)系分別進(jìn)行了3個(gè)實(shí)驗(yàn),其中實(shí)驗(yàn)1和實(shí)驗(yàn)2使用仿真數(shù)據(jù),實(shí)驗(yàn)3使用實(shí)測(cè)數(shù)據(jù)。
實(shí)驗(yàn)1 驗(yàn)證平行線約束。假設(shè)兩條直線 l1和l2為平行線,其斜率a1=a2=0.45,截距b1=1.6,b2=3.2。為l1給定7個(gè)沒有觀測(cè)誤差的x坐標(biāo)值(1.39,2.32,3.93,5.24,6.04,7.33,8.93),為l2給定8個(gè)沒有觀測(cè)誤差的x坐標(biāo)值(9.69,8.23,7.37,6.39,5.26,4.29,3.36,1.97),并根據(jù)斜率和截距計(jì)算對(duì)應(yīng)的y 坐標(biāo),單位均為m。然后為直線 l1和直線l2的x 坐標(biāo)附上期望為0、方差均為0.01 m2的隨機(jī)誤差,為直線 l1和直線l2的y 坐標(biāo)附上期望為0、方差均為0.008 1 m2的隨機(jī)誤差,得到了含有隨機(jī)誤差的坐標(biāo)觀測(cè)量,坐標(biāo)的權(quán)用單位矩陣乘方差的倒數(shù)表示。選擇表1中平行的約束條件模型進(jìn)行1 000次擬合,每次都采用LS、TLS、WTLS和ECWTLS 4種方法進(jìn)行解算。
4種方法進(jìn)行參數(shù)估計(jì)后,采用參數(shù)估計(jì)值減去給定的參數(shù)真值,得到估計(jì)值與真值的差值,即估計(jì)值誤差,其結(jié)果如圖1所示,再計(jì)算所有誤差絕對(duì)值的平均值以及總和,結(jié)果如表2和表3所列。另外,分別計(jì)算4種方法1 000次實(shí)驗(yàn)中參數(shù)誤差的標(biāo)準(zhǔn)差,用來表示參數(shù)誤差的離散程度,結(jié)果如表4所列。
由圖1可知,每次采用ECWTLS估計(jì)兩條直線斜率的估計(jì)值與真值之差都相等,而采用另外3種方法的差值不是每次都相等,由于此次實(shí)驗(yàn)設(shè)置的兩條直線斜率真值相等,可得采用ECWTLS估計(jì)兩條直線的斜率一定相等,LS、TLS、WTLS 3種方法估計(jì)兩條直線的斜率不一定相等。由表2和表3可知,采用ECWTLS所得參數(shù)誤差絕對(duì)值的均值與總和均相等,同樣說明在平行線約束下,ECWTLS估計(jì)的直線斜率相等;在4種方法中ECWTLS得到的誤差絕對(duì)值的平均值、誤差絕對(duì)值的總和均為最小,可以認(rèn)為4種方法中ECWTLS的估計(jì)值更接近真值。由表4可知,LS、TLS、WTLS 3種方法所得參數(shù)誤差的標(biāo)準(zhǔn)差比較接近,ECWTLS計(jì)算所得參數(shù)誤差的標(biāo)準(zhǔn)差最小,認(rèn)為采用ECWTLS估計(jì)平行線的參數(shù)時(shí)結(jié)果更加穩(wěn)定。
實(shí)驗(yàn)2 驗(yàn)證斜交約束。假設(shè)兩條直線 l3和l4,其中直線l3 的傾斜角為120°,斜率 a3為-3,截距b3為3.4;直線l4的傾斜角為15°,斜率a4為2-3,截距b4為2.8。為l3給定7個(gè)沒有觀測(cè)誤差的x 坐標(biāo)值(9.46,7.15,6.41,5.35,4.89,3.23,1.26),為 l4給定8個(gè)沒有觀測(cè)誤差的x 坐標(biāo)值(0.86,1.64,2.16,3.49,5.94,6.48,7.56,8.23),并根據(jù)斜率和截距計(jì)算對(duì)應(yīng)的 y 坐標(biāo),單位均為m。然后分別為直線 l3和直線l4的x坐標(biāo)附上期望為0、 方差均為0.01 m2的隨機(jī)誤差,分別為直線 l3和直線l4的y 坐標(biāo)附上期望為0、方差均為0.04 m2的隨機(jī)誤差,得到了含有隨機(jī)誤差的坐標(biāo)觀測(cè)量,坐標(biāo)的權(quán)用單位矩陣乘方差的倒數(shù)表示。選擇表1中斜交的約束條件模型,分別用LS、TLS、WTLS、ECWTLS 4種方法進(jìn)行直線參數(shù)估計(jì)。計(jì)算得到的參數(shù)估計(jì)值及單位權(quán)方差如表5所列,參數(shù)估計(jì)值與真值之差如表6所列。最后對(duì)擬合的兩條直線的夾角進(jìn)行計(jì)算,結(jié)果如表7所列。
由表5和表6可知,LS、TLS、WTLS 3種方法得到參數(shù)估計(jì)值與真值的差值差距不大,均大于ECWTLS,采用ECWTLS估計(jì)的參數(shù)最接近真值。由表7可知,采用ECWTLS擬合的兩條直線的夾角基本等于105°,其他3種方法擬合的兩條直線的夾角與105°有明顯差異,本文的方法能夠確保直線間的斜交約束關(guān)系。通過實(shí)驗(yàn)2的結(jié)果,說明ECWTLS通過對(duì)直線間的夾角進(jìn)行約束,能夠提高直線參數(shù)估計(jì)的準(zhǔn)確性。
實(shí)驗(yàn)3 驗(yàn)證垂直約束。選取兩條設(shè)計(jì)夾角為90°的道路,采用全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS,global navigation satellite system)測(cè)量道路中心線上的點(diǎn),認(rèn)為是等精度觀測(cè),實(shí)測(cè)坐標(biāo)如表8所列。仍然采用LS、TLS、WTLS、ECWTLS 4種方法進(jìn)行直線參數(shù)估計(jì),參數(shù)估計(jì)結(jié)果如表9所列,4種方法擬合的兩條直線的夾角如表10所列。
由表9可知,由于是等精度觀測(cè),TLS和WTLS兩種方法得到的參數(shù)估計(jì)值完全相同。另外,由于LS、TLS、WTLS 3種方法是分開估計(jì)兩條直線的參數(shù),所得參數(shù)的單位權(quán)方差不相同,而ECWTLS是同時(shí)估計(jì)兩條直線的參數(shù),所得直線參數(shù)的單位權(quán)方差相同。由表10可知,采用ECWTLS擬合得到兩條直線的夾角基本等于90°,其他3種方法擬合出的兩條直線的夾角均與90°有一定的差距。實(shí)驗(yàn)3的結(jié)果說明本文提出的直線間垂直約束模型能夠確保擬合直線的正交關(guān)系。
4 結(jié)語
以直線間的平面關(guān)系為出發(fā)點(diǎn),分別以平行、斜交和垂直3種位置關(guān)系推導(dǎo)了不同約束條件的約束模型,并與WTLS相結(jié)合組成ECWTLS,針對(duì)平行約束和斜交約束進(jìn)行模擬數(shù)據(jù)實(shí)驗(yàn),針對(duì)垂直約束進(jìn)行實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn),得出以下結(jié)論:
(1) 本文給出的方法和算法具有可行性。
(2) 在已知直線位置關(guān)系時(shí),與LS、TLS、WTLS相比,使用ECWTLS能確保擬合直線位置關(guān)系正確。
(3) 擬合平行直線時(shí),ECWTLS計(jì)算的參數(shù)標(biāo)準(zhǔn)差最小,參數(shù)估計(jì)值最接近真值;擬合斜交直線時(shí),ECWTLS計(jì)算的參數(shù)估計(jì)值最接近真值。
參考文獻(xiàn):
[1] KRYSTEK M,ANTON M.A least-squares algorithm for fitting data points with mutually correlated coordinates to a straight line[J].Measurement Science and Technology,2011,22(3):035101.
[2] PETROLINI A.Linear least squares fit when both variables are affected by equal uncorrelated errors[J].American Journal of Physics,2014,82(12):1178-1185.
[3] 李雄軍.對(duì)X和Y方向最小二乘線性回歸的討論[J].計(jì)量技術(shù),2005(1):50-51.
[4] 梁家惠.用最小二乘法進(jìn)行直線擬合的討論[J].工科物理,1995(3):11-15.
[5]HUFFEL S V,VANDEWALLE J.The total least squares problem:computational aspects and analysis[M].Philadelphia:Society for Industrial and Applied Mathematics,1992.
[6] GOLUBG H,VANLOAN C F.An analysis of the total least squares problem[J].Siam Journal on Numerical Analysis,1980,17(6):883-893.
[7] 魯鐵定,周世健.總體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2010,35(11):1351-1354.
[8] 邱衛(wèi)寧,齊公玉,田豐瑞.整體最小二乘求解線性模型的改進(jìn)算法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2010,35(6):708-710.
[9] 孔建,姚宜斌,吳寒.整體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2010,35(6):711-714.
[10] SCHAFFRIN B.A note on constrained total least-squares estimation[J].Linear Algebra and Its Applications.2006,417(1):245-258.
[11] SCHAFFRIN B,WISER A.On weighted total least-squares adjustment for linear regression[J].Journal of Geodesy,2008,82(7):415-421.
[12] 丁克良,沈云中,歐吉坤.整體最小二乘法直線擬合[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,29(1):44-47.
[13] 丁克良,歐吉坤,趙春梅.正交最小二乘曲線擬合法[J].測(cè)繪科學(xué),2007,147(3):18-19,192.
[14] SHEN Y,LI B,CHEN Y.An iterative solution of weighted total least-squares adjustment[J].Journal Geodesy,2011,85(4):229-238.
[15] NERI F,SAITTA G,CHIOFALO S.An accurate and straightforward approach to line regression analysis of error-affected experimental data[J].Journal Physics E:Scientific Instruments,1989,22(4):215-217.
[16] AMIRI-SIMKOOEI A,JAZAERI S.Weighted total least squares formulated by standard least squares theory[J].Journal of Geodetic Science,2012,2(2):113-124.
[17] 王繼剛,周立,蔣廷臣等.一種簡(jiǎn)單的加權(quán)整體最小二乘直線擬合方法[J].測(cè)繪通報(bào),2014(4):33-35.
[18] SCHAFFRIN B,F(xiàn)ELUS Y.An algorithmic approach to the total least-squares problem with linear and quadratic constraints[J].Studi Geophysica Camerti,2009,53(1):1-16.
[19] MAHBOUB V.On weighted total least-squares with linear and quadratic constraints[J].Journal of Geodesy,2013,87(3):279-286.
[20] FANG X.Weighted total least-squares with constraints:a universal formula for geodetic symmetrical transformations[J].Journal of Geodesy,2015(89):459-469.
[21] 王樂洋.附有等式約束的加權(quán)總體最小二乘平差方法[J].東華理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,36(2):245-248.
[22] 胡川,方興,趙立都.融合正交幾何信息的非線性等式約束整體最小二乘平差及迭代算法[J].測(cè)繪學(xué)報(bào),2020,49(7):816-823.
[23] XU P,LIU J,SHI C.Total least squares adjustment in partial errors-in-variables models:algorithm and statistical analysis[J].Journal of Geodesy,2012,86(8):661-675.
[24] HU C,CHEN Y,ZHU W D.Generalized total least squares solution based on pseudo-observation method[J].Survey Review,2016,48(348):157-167.
An algorithm for total least squares of line parameters
considering the position relationship
YAN Xiaoyu,LIU Yufan,LU Ziyang,ZHANG Chongyang
(School of Smart City,Chongqing Jiaotong University,Chongqing 400074,China)
Abstract
Position relationships between lines can be divided into parallel,perpendicular and intersecting,and these relationships can be transformed into constraints in line fitting.For taking advantage of this constraint,the weighted total least squares estimation method with equation constraint is introduced for line fitting.Position relationship between lines is expressed as the relationship between slopes of lines,and is linearized to obtain three constraint models,which together with the weighted total least squares criterion determine the equality constrained weighted total least squares model.Then the solution
formula of the model is derived and an iterative algorithm is proposed.Finally,three simulation experiments demonstrated the feasibility and effectiveness of the developed algorithm.The result shows that the new algorithm is feasible with a higher fitting accuracy,and can ensure the correct plane position relationship between straight lines.
Key words
Plane position relationship;Straight-line fiting;Weighted total least squares;Equality constraint
(本文責(zé)編:葛 文)