Algoritmus určování rovnice roviny pro laserové skenování Ing. Bronislav Koska, Ing. Martin Štroner, Ph.D., Doc. Ing. Jiří Pospíšil, CSc., ČVUT - Fakulta stavební, Praha.
1 Úvod V rámci řešení projektu GA ČR „Moderní optoelektronické metody topografie ploch“ je hledána vhodná metoda a postup bezkontaktního měření tvaru objektu s dostatečnou přesností a v krátkém čase. Navrhované optoelektronické metody jsou založeny na principu detekce optického signálu generovaného vhodným zdrojem optického záření a jeho následném počítačovém vyhodnocení. Na základě takto získaných informací bude možno vytvořit skelet měřeného objektu, který může být počítačově domodelován. V této etapě řešení je navrhován systém pro skenování prostorově členitého povrchu, který jako zdroje optického záření využívá laserového záření ve viditelné oblasti spektra, jako detekčního zařízení digitální barevné CCD kamery a dále zařízení pro otáčení zkoumaného předmětu. Pro určení tvaru povrchu bude využívána světelná rovina vytvořená laserem a optickým členem. Průnik roviny s předmětem vytvoří světelný profil. Tento profil bude snímán CCD kamerou a dalším vyhodnocením a zpracováním budou získány prostorové souřadnice bodů tohoto profilu. V tomto příspěvku je řešena problematika určení rovnice světelné roviny na základě 3D souřadnic bodů, které budou měřeny geodetickými metodami v nadbytečném počtu.
2 Základní používané pojmy Pro řešení byla použita metoda vyrovnání podmínkových měření s neznámými pro korelovaná měření viz. [1]. Vychází se z obecné rovnice roviny: A⋅ x + B ⋅ y + C ⋅ z + D = 0 ,
(1)
která se upraví na tvar s minimálním počtem jednoznačně definovaných neznámých: a ⋅ x + b ⋅ y + c ⋅ z +1 = 0 .
(2)
Nevýhodou tohoto tvaru je nemožnost vyjádření roviny procházející počátkem. Výchozím vztahem je výraz pro vzdálenost bodu od roviny (2), který je dle [2] definován takto: d=
a ⋅ x + b ⋅ y + c ⋅ z +1 a 2 + b2 + c2
a který je pro body ležící v rovině, tedy vyrovnané body, roven nule.
1
(3)
3n
Podmínka pro vyrovnání je Ω = ∑ vi2 = min , kde n je počet bodů. Formulaci můžeme i =1
rozšířit pro body zaměřené s různou přesností a zapsat maticově: Ω = v T ⋅ P ⋅ v = min .
(4)
3 Příprava a formát vstupních hodnot Vstupní hodnoty vyrovnání jsou souřadnice bodů. Ty se v našem případě neměří přímo, ale počítají se z měřených délek a směrů. K získání středních chyb hodnot vstupujících do vyrovnání je tedy nutno vycházet z těchto vztahů a aplikovat na ně zákon hromadění středních chyb. Dalším aspektem je podchycení matematické korelace mezi takto vypočtenými souřadnicemi. Vstupní veličinou vyrovnání tedy jsou souřadnice bodů a jejich kovarianční matice, resp. matice váhových koeficientů a střední chyba jednotková nebo matice vah a střední chyba jednotková. Z důvodů konvence použité v [1] použijeme pro naše měření (souřadnice) matici váhových koeficientů, protože se jedná o měření korelovaná. Výchozí vztahy jsou: X = X 0 + l ⋅ sin z ⋅ cos α Y = Y0 + l ⋅ sin z ⋅ sin α ,
(5)
Z = Z 0 + l ⋅ cos z kde
X0, Y0, Z0 .. jsou souřadnice počátku měření délek a úhlů, l .. šikmá délka, z .. zenitový úhel, α .. směrník.
Dle [1] je zákon hromadění středních chyb definován takto: Sh = H ⋅ M 2 ⋅ H T .
(6)
Podle (6) pro kovarianční matici souřadnic platí: 2 S XYZ = H ⋅ M mer ⋅ HT ,
kde
(7)
SXYZ .. je kovarianční matice souřadnic, H.. matice koeficientů lineární funkce přenosu matice skutečných chyb měření na matici skutečných chyb neznámých. Vzniká derivací vztahů (5) podle jednotlivých měření, M2mer .. matice variancí (na diagonále jsou kvadráty středních chyb měření).
V matici středních chyb měření vystupují střední chyby délek, směrníků a zenitových úhlů: ml, mα, mz. Ty jsou pro danou situaci uvažovány pro všechny body stejné. Pro další výpočty je nutno vyjádřit matici váhových koeficientů. K tomu je potřeba volit stření chybu jednotkovou, např. jako průměrnou hodnotu kvadrátů středních chyb všech souřadnic. Matice váhových koeficientů:
2
Q XYZ =
1 S XYZ , m02
(8)
kde mo2 volíme např. jako průměr hodnot diagonálních prvků matice SXYZ. Protože
∂X i = 0 pro i ≠ j , analogicky pro ostatní funkce a proměnné, mají funkce H a ∂l j
SXYZ tvar: ∂X 1 ∂l 1 ∂Y1 ∂l 1 ∂Z1 H = ∂l1 0 . 0
S XYZ
mX2 1 mY1 , X1 mZ , X = 1 1 0 . 0
∂X 1 ∂α1
∂X 1 ∂z1
0
0
.
0
∂Y1 ∂α1
∂Y1 ∂z1
0
0
.
0
∂Z1 ∂α1
∂Z1 ∂z1
0
0
.
0
0
0
∂X 2 ∂l2
∂X 2 ∂α 2
.
0
.
.
.
.
.
0
0
0
0
0 0 0 , 0 . ∂Z n ∂zn
. ∂Z n . ∂α n
mX1 ,Y1
mX1 , Z1
0
0
.
0
mY21
mY1 , Z1
0
0
.
0
0
0
.
0
.
0
mZ1 ,Y1
2 Z1
m
2 Y2
0
0
m
mY2 , X 2
. 0
. 0
. 0
. 0
. . . mZ n ,Yn
0 0 0 . 0 . mZ2n
(9)
(10)
4 Výpočet přibližných hodnot neznámých a,b,c Vypočet je možno provést řešením rovnice roviny (2) pro nutný počet bodů (tři body). Elegantnější řešení rovnice roviny v obecném tvaru (1) nabízí Bronštejn - Semenďajev v [3]. Řešíme rovnici s determinantem: x − x1 x2 − x1 x3 − x1
y − y1 y2 − y1 y3 − y1
z − z1 z2 − z1 = 0 , z3 − z1
(11)
a tedy A = ( y2 − y1 ) ⋅ ( z3 − z1 ) − ( y3 − y1 ) ⋅ ( z 2 − z1 ) ,
(12)
B = − [ ( x2 − x1 ) ⋅ ( z3 − z1 ) − ( x3 − x1 ) ⋅ ( z 2 − z1 )] ,
(13)
C = ( x2 − x1 ) ⋅ ( y3 − y1 ) − ( x3 − x1 ) ⋅ ( y2 − y1 ) ,
(14)
D = − x1 ⋅ A − y1 ⋅ B − z1 ⋅ C .
(15)
3
Ke stejným vztahům lze také dojít aplikací Cramerova pravidla. Při programovém řešení je vhodné počítat tyto přibližné koeficienty z několika různých kombinací bodů s vyloučením odlehlých měření. Je to vhodné z důvodu neřešitelnosti této rovnice pro body ležící na přímce (všechny koeficienty jsou nulové) a z důvodu možného zavádějícího řešení pro body přímce se blížící. Jiné jednodušší ověření vhodnosti volby bodů pro přibližné řešení je výpočet přímky ze dvou bodů a následný výpočet vzdálenosti třetího bodu od této přímky. Zde by mohlo dojít např. k porovnání této vzdálenosti bodu od přímky se vzdáleností dvou původních bodů. Toto pravidlo je ale oprávněné pouze při přibližně pravidelné síti měřených bodů. Matematická podmínka splněná pro vhodné body může být např.: d 3, p >
d1, 2 , kde 10 d3,p .. vzdálenost třetího bodu od přímky d1,2 .. vzdálenost prvního a druhého bodu
Dle [2] můžeme vzdálenost d3,p vyjádřit: r uuuur u × ( P1 P3 ) r d3, P = , kde u = P2 − P1 . r u
(16)
A d1,2 podle známého vztahu: d1,2 =
( x2 − x1 ) + ( y2 − y1 ) + ( z2 − z1 ) 2
2
2
.
(17)
= 0.
(18)
5 Linearizace a maticová formulace Výchozí podmínka pro vyrovnání je dána rovnicí: f ( x , y , z , a , b, c ) =
a ⋅ x + b ⋅ y + c ⋅ z +1 a 2 + b2 + c 2
V této rovnici představují, poněkud nezvykle, měření proměnné x,y,z a hledané neznámé koeficienty a,b,c. Počet podmínek odpovídá počtu zaměřených bodů. Všechny podmínky můžeme maticově zapsat takto: f(xT , hT ) = 0 ,
(19)
kde xT je matice všech souřadnic: x T = ( x1 , y1 , z1 , x2 ,...., z n ) a hT matice neznámých: hT = (a, b, c) . Podmínky můžeme linearizovat a vyjádřit maticově: AT ⋅ v + B ⋅ dh + u = 0 , kde:
4
(20)
∂a1 P Q R 0 0 . 0 0 ∂a 0 0 0 P Q . 0 0 2 T ,B = . A = . . . . . . . . . 0 0 0 0 0 . Q R ∂a n
(
v T = vx1
v y1
vz1
∂b1 ∂b2 . . ∂bn
∂c1 u01 u ∂c2 da 02 . , dh = db , u = . , (21) dc . . u ∂cn 0n
)
. vzn ,
vx2
kde: a0 b0 c0 ∂fi ∂f ∂f = = P, i = =Q, i = = R, ∂xi ∂yi ∂zi a02 + b02 + c02 a02 + b02 + c02 a02 + b02 + c02 ∂fi = ∂ai = ∂a ∂fi = ∂bi = ∂b ∂fi = ∂ci = ∂c
xi
(a
2 0
1 2 2 0
)
+b +c 2 0
yi
(a
2 0
1 2 2 0
+b +c 2 0
)
zi
(a
2 0
1 2 2 0
+b +c 2 0
u0 i =
)
−
a0 ⋅ xi + b0 ⋅ yi + c0 ⋅ zi + 1
(a
2 0
−
(a
3 2 2 0
+b +c 2 0
)
a0 ⋅ xi + b0 ⋅ yi + c0 ⋅ zi + 1
(a
2 0
3 2 2 0
+b +c 2 0
a0 ⋅ xi + b0 ⋅ yi + c0 ⋅ zi + 1
6 Hledání řešení podmínky
)
+b +c
a0 ⋅ xi + b0 ⋅ yi + c0 ⋅ zi + 1 2 0
−
3 2 2 0
2 0
a02 + b02 + c02
)
⋅ a0 ,
⋅ b0 ,
⋅ c0 ,
.
-1 Ω = v T ⋅ QXYZ ⋅ v = min
Použijeme tzv. „přímé řešení“ podle Lagrangeova postupu hledání minima: -1 Ω = v T ⋅ Q XYZ ⋅ v − 2 ⋅ k T ⋅ (AT ⋅ v + B ⋅ dh + u) = min ,
kde k je vektor zatím neurčených tzv. Lagrangeových koeficientů. Minimum této funkce získáme pokud položíme její derivaci podle v a dh rovnu nule. Maticově tyto rovnice můžeme zapsat takto: AT ⋅ Q XYZ ⋅ A B k u ⋅ + = 0 . T B 0 dh 0 Pokud řešení vyjádříme pomocí inverzní matice: −1
AT ⋅ Q XYZ ⋅ A B u k = − ⋅ . BT 0 0 dh Pro další formulace je vhodné vyjádřit tuto rovnici ve zkrácené formě: y = − N −1 ⋅ n a také provést další zjednodušení v zápise:
5
(22)
N 0 = AT ⋅ QXYZ ⋅ A . Podle pravidel pro inverzi blokové matice lze inverzi matice N-1 vyjádřit obecně takto: N
−1
N 0−1 − N 0−1 ⋅ B ⋅ (BT ⋅ N 0−1 ⋅ B)−1 ⋅ B T ⋅ N 0−1 = (B T ⋅ N 0−1 ⋅ B)−1 ⋅ B T ⋅ N 0−1
N 0−1 ⋅ B ⋅ (B T ⋅ N 0−1 ⋅ B)−1 −(BT ⋅ N 0−1 ⋅ B)−1
a opět zjednodušit značení: Q N −1 = kk Qhk
Qkh . Qhh
(23)
Pomocí těchto matic lze již přímo vyjádřit řešení neznámých: dh = −Qhk ⋅ u . Ale i hodnoty korelát a opravy měření: k = −Qkk ⋅ u , v = Q XYZ ⋅ A ⋅ k .
7 Rozbor přesnosti Pro vyrovnané neznámé a,b,c je nutno určit jejich střední chyby a také jejich vzájemné korelace, pro potřeby dalších výpočtů s nimi. Princip jejích určení je založen na zákonu hromadění vah na množině funkcí publikovaném např. v [1], kde je také odvozen výpočet matice váhových koeficientů pro neznámé u výše popsaného způsobu vyrovnání. Bez odvození uvádíme výsledný vztah: Qh = −Qhh . Výpočet Qhh je popsán v předcházející kapitole (vzorec (23)). Výpočet středních chyb jednotlivých vypočtených neznámých se tedy uskuteční podle vzorce: ma = m0 ⋅ qa , kde qa Qh = q b ,a q c ,a
qa ,b qb qc ,b
qa ,c q b ,c . qc
Výpočet aposteriorní střední chyby jednotkové m0 se provede podle vzorce: m0 =
−1 v T ⋅ Q XYZ ⋅v , r −k
(24)
kde r .. počet podmínek a k .. je počet neznámých. Výpočet kovarianční matice potřebné pro další výpočty s neznámými a,b,c se tedy provede podle vzorce: Sh = m0 ⋅ Qh .
6
(25)
8 Kontrola výpočtu Pro kontrolu linearizace a ukončení výpočtu bylo stanoveno pravidlo, aby 1/10 střední chyby byla větší než absolutní hodnota přírůstku příslušné neznámé. Tím by mělo být zabezpečeno, že chyba z linearizace, neboli z nepřesnosti přibližných hodnot neznámých příliš neovlivní výsledné neznámé. Symbolicky: ma m m > da ∧ b > db ∧ c > dc . 10 10 10
(26)
9 Algoritmus výpočtu 1) Výpočet souřadnic vstupujících do vyrovnání a jejich kovarianční matice podle kapitoly 3. 2) Výpočet přibližných hodnot neznámých způsobem uvedeným v kapitole 4. 3) Výpočet vyrovnaných neznámých podle vzorců v kapitole 6. 4) Výpočet středních chyb a kovarianční matice vyrovnaných neznámých podle kapitoly 7. 5) Kontrola výpočtu podle nerovnic v odstavci 8. Pokud jsou nerovnice splněny, ukončení výpočtu. Pokud ne, opakování výpočtů od bodu 2). Za přibližné hodnoty neznámých jsou přiřazeny hodnoty, které jsou výsledkem vyrovnání.
10 Příklad Pro použití uvedené metody v projektu bude vypracován program umožňující její rychlé a flexibilní použití.
10.1 Provedené měření Vstupní hodnoty do níže uvedeného výpočtu byly získány na základě měření. To bylo provedeno z několika důvodů. Zejména se jednalo o: 1. Odhalení technických realizované laserem.
problémů
v zamýšlené
metodě
zaměření
roviny
2. První přibližné zjištění skutečných středních chyb určovaných koeficientů, což je nezbytné pro reálný odhad přesnosti celého systému. 3. Ověření odhadů středních chyb vstupujících veličin na základě vyrovnání. Při měření byly simulovány podmínky, které jsou předpokládané pro vyvíjený systém. Jednalo se zejména o rozměr a prostorovou konfiguraci. Byly také použity přístroje, jejichž použití se v systému předpokládá. Byla zaměřena přibližně pravidelná matice bodů 5x6 (sloupce x řádky) o rozměru přibližně (1.1 x 0.6)m. Bližší informace o měření by měly být uveřejněny v některém z dalších článků.
10.2 Výpočet Nyní uvádíme příklad vypočtený podle výše uvedeného algoritmu v programu Mathcad 2001i Professional.
7
Bylo zaměřeno 30 bodů prostorovou polární metodou. Uvažovaná střední chyba měřeného směrníku je mα = 0.0030gon , zenitového úhlu mz = 0.0030gon a měřené délky je ml = 0.0006m . Po prvním kontrolním výpočtu byla zjištěna řádově větší vzdálenost některých bodů od vyrovnané roviny. Jednalo se o body číslo 5,9 a 10. U těchto bodů lze předpokládat hrubou chybu a proto byly z konečného výpočtu vyloučeny. Pro níže uvedené výpočty bylo tedy použito 27 bodů. Jelikož se počítá s maticemi velkých rozměrů (až 81x81), jsou v článku někdy uvedeny pouze zmenšené verze těchto matic, vždy tak, aby jasně vyjadřovaly tvar a trend matice. Naměřené veličiny po všech nutných redukcích jsou (podle sloupců: směrník [gon], zenitový úhel [gon], prostorová délka [m])(zobrazena je matice pro prvních 8 bodů): 80.6794 77.9014 74.2232 85.3904 92.2478 92.2342 86.0527 80.6832
103.2759 103.0358 102.7201 103.6396 104.1519 105.6659 105.0451 104.4805
2.2969 2.4670 2.7451 2.0638 . 1.8167 1.8177 2.0368 2.2989
Z nich byly vypočteny prostorové souřadnice těchto bodů podle (5) a podle zákona hromadění vah také jejich matice váhových koeficientů (zobrazeno prvních 8 bodů): 11.080 10.469 10.646 10.686 X = 10.838 10.220 10.441 10.685
12.521 12.006 12.152 12.189 12.317 11.797 11.982 12.188
9.883 9.882 9.685 9.882 . 9.882 9.838 9.839 9.838
Z důvodu velkého rozměru matice váhových koeficientů je uveden pouze její výřez pro šest souřadnic (první dva body). I z něj je zřejmá matematická korelace mezi souřadnicemi jednotlivých bodů:
Q XYZ
0.545999 0.969444 -0.045089 = 0.000000 0.000000 0.000000
0.969444 2.392490 -0.105205 0.000000 0.000000 0.000000
-0.045089 -0.105205 0.135645 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.214943 0.604793 -0.035543
0.000000 0.000000 0.000000 0.604793 2.662647 -0.152151
První výpočet: Přibližné hodnoty neznámých vypočtené z prvních tří bodů:
8
0.000000 0.000000 0.000000 . -0.035543 -0.152151 0.082846
a0 = 0.27765 b0 = -0.33015 . c0 = 0.00580
Matice AT (je zobrazena opět jenom část matice pro první tři body) : 0.643572 -0.765268 0.013439 A = 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 T
0.000000 0.643572 0.000000
0.000000 -0.765268 0.013439 0.000000 0.000000 0.000000 . 0.000000 0.000000 0.643572 -0.765268 0.013439
0.000000
0.000000 0.000000 0.000000
Matice B (prvních 8 bodů): 25.684 24.266 24.677 24.769 B = 25.123 23.691 24.202 24.769
29.023 27.830 28.168 28.253 28.551 27.345 27.774 28.251
22.908 22.906 22.450 22.906 . 22.907 22.805 22.806 22.805
Matice N (část matice): 0.673767 0.000000 0.000000 N = 0.000000 0.000000 0.000000 0.000000
0.000000 1.055167 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.895553 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.878873 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.784284 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 1.330378 0.000000
0.000000 0.000000 0.000000 0.000000 . 0.000000 0.000000 1.079124
Přírůstky a vyrovnané neznámé jsou: 2.277e - 3 dh = -2.922e - 3 a 1.155e - 3
0.279924 h = -0.333071 . 0.006953
Střední chyba jednotková aposteriorní a matice váhových koeficientů: m0 = 0.000337
a
2.254876 Qh = -2.396273 0.527269
-2.396273 2.584269 -0.608179
0.527269 -0.608179 . 0.184029
Střední chyby a jejich porovnání s přírůstky dh: ma = 0.000506 < da ⋅10 = 0.022777 , mb = 0.000541 < db ⋅10 = 0.029227 , mc = 0.000144 < dc ⋅ 10 = 0.011552 .
Protože požadovaná podmínka není splněna pro všechny neznámé je nutno výpočet opakovat.
9
Druhý výpočet: Jako vstupující přibližné hodnoty jsou převzaty výsledné hodnoty z předcházejícího kroku (prvního). Pro stručnost budou uvedeny pouze výsledky a porovnání s podmínkou: Přírůstky a vyrovnané neznámé jsou: 1.975e - 5 dh = -2.530e - 5 a 9.961e - 6
0.279944 h = -0.333096 . 0.006963
Střední chyba jednotková aposteriorní a matice váhových koeficientů: m0 = 0.000337
a
2.334777 Qh = -2.485306 0.551063
-2.485306 2.684080 -0.635460
0.551063 -0.635460 . 0.192094
Střední chyby a jejich porovnání s přírůstky dh: ma = 0.000514 > da ⋅ 10 = 0.000198 , mb = 0.000551 > db ⋅10 = 0.000253 , mc = 0.000148 > dc ⋅ 10 = 0.000100 .
Protože podmínka je nyní splněna pro všechny neznámé, jsou dosažené výsledky považovány za správné. Pro další výpočty budou tedy použity vyrovnané hodnoty a jejich kovarianční matice Sh: a = 0.27994 0.00017555 b = -0.33310 a Sh = -0.00019554 0.00005204 c = 0.00696
-0.00019554 0.00021944 -0.00005995
0.00005204 -0.00005995 . 0.00001787
10
Příspěvek vznikl v rámci řešení projektu GA ČR 103/02/0357. Literatura: [1] Böhm J., Radouch V. a Hampacher M.: Teorie chyb a vyrovnávací počet. Praha: GKP, 1990, 416 s. [2] Budinský B., Charvát J.: Matematika I. Skriptum. Praha: ČVUT, 1994. [3] Bronštejn L.N., Semenďajev K.A.: Príručka matematiky pre inžinierov a pre studujůcích. Bratislava: Slovenské vydavatelstvo technickej literatúry, 1964.
Anotace: V příspěvku je řešena problematika určení rovnice světelné roviny, vyvíjeného laserového skeneru, na základě 3D souřadnic bodů, které budou měřeny prostorovou polární metodou v nadbytečném počtu. Pro řešení je užita metoda vyrovnání podmínkových měření s neznámými pro korelovaná měření.
11
