• Wykład 4 Modele macierzy kowariancyjnej obserwacji w rachunku wyrównawczym

Test porównanie oceny „a priori” i „a posteriori

Ocena „a priori” tworzona jest na podstawie wiedzy przed wykonaniem obserwacji. Możemy tu wyróżnić dwa przypadki:

• Posiadamy tylko wiedzę teoretyczną pozwalającą na sformułowanie hipotezy badawczej określającej relacje parametrów, przy czym mogą wystąpić parametry o nieznanych wartościach. Eksperyment pomiarowy ma na celu zweryfikowanie hipotezy i dostarczenie oszacowań dla wszystkich parametrów.

• Posiadamy wiedzę teoretyczną i praktyczną pozwalającą na oszacowanie wartości parametrów obserwacji wykonywanych znanymi metodami w powtarzalnych warunkach. W szczególności inżynier nie tylko powinien potrafić określić dokładność pomiaru rutynowego ale dobrać metodę pomiaru do wymogów wynikających z przeznaczenia zleconych pomiarów.

Ocena „a posteriori” tworzona jest na podstawie wyników wykonanych obserwacji (próbki statystycznej) w postaci estymatora parametru.

Przykład

Mamy wykonać wielokrotny pomiar tej samej wielkości l₁, l₂, …, l_n. Dokładność każdej obserwacji charakteryzuje odchylenie standardowe o apriorycznej wartości σ = 0.01.

Wykonaliśmy 5 obserwacji 7.03, 7.02, 7.09, 7.04, 7.03

Estymator wartości oczekiwanej (wartość najprawdopodobniejsza) wynosi

$$\hat{m} = \overset{\overline{}}{l} = \frac{1}{n}\sum_{i = 1}^{n}l_{i} = 7.042$$

Estymator wariancji (wariancja z próby) wynosi

$${\hat{\mu}}_{2} = {\hat{\sigma}}^{2} = \frac{1}{n - 1}\sum_{i = 1}^{n}\left( l_{i} - \overset{\overline{}}{l} \right)^{2} = \frac{1}{n - 1}\sum_{i = 1}^{n}{\hat{v}}^{2} = 0.00077$$

Stąd ocena „a posteriori” odchylenia standardowego wynosi $\hat{\sigma} = 0.028$

Test statystyczny ma odpowiedzieć na pytanie czy z zadanym prawdopodobieństwem (np. 90%) taka rozbieżność w stosunku do apriorycznej wartości σ = 0.01 jest dopuszczalna.

Modele macierzy kowariancyjnej obserwacji w rachunku wyrównawczym

Modele charakterystyk dokładnościowych posiadają długą tradycję w geodezyjnym rachunku wyrównawczym. Związane jest to z wyrównaniem obserwacji niezależnych których dokładność charakteryzują wariancje tworzące diagonalną macierz kowariancyjną oraz uproszczeniami rachunkowymi wynikającymi z wprowadzenia pojęcia wag. Obecnie powszechne wykorzystanie komputerów w obliczeniach geodezyjnych pozwala na rezygnację z uproszczeń rachunkowych. Pozostaje jednak dalej rozróżnienie pomiarów badawczych i rutynowych. Pomiary badawcze, wykonywane często przy współpracy z naukowcami innych branż, powinny być opracowywane powszechnie uznanymi metodami statystycznymi – a te są bardzo zbliżone do klasycznego rachunku z wykorzystaniem wag obserwacji.

Model Gaussa-Markowa - współczynnik wariancji, macierz kofaktorów, macierz wag

Macierz kowariancyjna dla niezależnych obserwacji ma postać diagonalną

$$\mathbf{C} = diag\left( \sigma_{i}^{2} \right) = \begin{bmatrix} \begin{matrix} \sigma_{1}^{2} & 0 \\ 0 & \sigma_{2}^{2} \\ \end{matrix} & \begin{matrix} \ldots & 0 \\ \ldots & 0 \\ \end{matrix} \\ \begin{matrix} \vdots & \vdots \\ 0 & 0 \\ \end{matrix} & \begin{matrix} \ddots & \vdots \\ \ldots & \sigma_{n}^{2} \\ \end{matrix} \\ \end{bmatrix}$$

$$\text{gdzie\ \ \ \ \ \ }\mathbf{Q} = diag\left( q_{i} \right) = \begin{bmatrix} \begin{matrix} q_{1} & 0 \\ 0 & q_{2} \\ \end{matrix} & \begin{matrix} \ldots & 0 \\ \ldots & 0 \\ \end{matrix} \\ \begin{matrix} \vdots & \vdots \\ 0 & 0 \\ \end{matrix} & \begin{matrix} \ddots & \vdots \\ \ldots & q_{n} \\ \end{matrix} \\ \end{bmatrix}$$

Realizacja obliczeń wyrównawczych wymaga znajomości odwrotności macierzy kowariancyjnej obserwacji ale wystarczy macierz wag P proporcjonalna do odwrotności macierzy kofaktorów

$$\mathbf{C}^{\mathbf{- 1}} \approx \mathbf{P} = c\mathbf{Q}^{\mathbf{- 1}} = diag\left( \frac{c}{q_{i}} \right)$$

$$\mathbf{P} = diag\left( p_{i} \right) = \begin{bmatrix} \begin{matrix} p_{1} & 0 \\ 0 & p_{2} \\ \end{matrix} & \begin{matrix} \ldots & 0 \\ \ldots & 0 \\ \end{matrix} \\ \begin{matrix} \vdots & \vdots \\ 0 & 0 \\ \end{matrix} & \begin{matrix} \ddots & \vdots \\ \ldots & p_{n} \\ \end{matrix} \\ \end{bmatrix}$$

Wagi p_i są więc proporcjonalne do odwrotności kwadratów błędów średnich i najłatwiej je oszacować w przypadku gdy wszystkie obserwacje odnoszą się do tej samej wielkości geometrycznej np. przewyższeń w niwelacji geometrycznej wówczas wagi są odwrotnie proporcjonalne do długości ciągów niwelacyjnych a współczynnik wariancji jest kwadratem błędu średniego obserwacji jednostkowej (1 km niwelacji) – jego estymator (oszacowanie „a posteriori”) wynosi

${\hat{\sigma}}_{0}^{2} = \frac{1}{n - u}\ {\hat{v}}^{T}P\hat{v}$ gdzie u jest liczbą niewiadomych parametrów

Estymatory wariancji obserwacji wynoszą ${\hat{\sigma}}_{i}^{2} = {\hat{\sigma}}_{0}^{2}{*q}_{i}$

Analogicznie mnożymy przez estymator współczynnika wariancji ${\hat{\sigma}}_{0}^{2}$ wszystkie estymatory wariancji dla funkcji obserwacji np. parametrów i obserwacji wyrównanych – jest to charakterystyczna cecha modelu Gaussa-Markowa wynikająca modelowania dokładności obserwacji nieznanym współczynnikiem wariancji.

Refleksja dotycząca empirycznego określania dokładności obserwacji

Wariancja resztkowa, jako suma kwadratów k = n − u niezależnych składników o rozkładach normalnych N(0,σ₀), jest zmienna losową o rozkładzie χ² posiadającym k stopni swobody. Parametrami tego rozkładu są

m(χ²) = k σ₀² ; μ₂(χ²)=2k σ₀⁴.

Fisher (1925) wykazał, że gdy liczba stopni swobody k wzrasta do nieskończoności, to zmienna losowa $\sqrt{2\chi^{2}}$ ma rozkład asymptotyczny normalny $N\left( \sqrt{2k - 1},\ 1 \right)$; a dla k>30 można już korzystać z tego rozkładu granicznego. Wynika stąd, że estymator odchylenia standardowego

$${\hat{\sigma}}_{o} = \sqrt{\frac{\chi^{2}}{k}} = \sqrt{\frac{{\hat{v}}^{T}P\hat{v}}{n - u}}$$

ma asymptotyczny rozkład normalny $N\left( \sigma_{o}\sqrt{1 - \frac{1}{2k}},\frac{\sigma_{o}}{\sqrt{2k}} \right)$.

Pozwala to na stwierdzenie, że:

• dla małej liczby obserwacji estymator niedoszacowuje wartość odchylenia standardowego;

• uzyskanie „trzysigmowej” dokładności względnej określonej wskaźnikiem $q = \frac{100\%}{p\%}$ czyli

$$\frac{3\sigma}{m} = \frac{\frac{3\ \sigma_{o}}{\sqrt{2k}}}{\sigma_{o}\sqrt{1 - \frac{1}{2k}}} = \frac{3}{\sqrt{2k - 1}} < \frac{1}{q}$$

wymaga posiadania $k > \ \frac{9}{2}q^{2}$ obserwacji, a w szczególności uzyskanie

• p = 10%;  q = 10;    wymaga   450  obserwacji

• p = 20%;   q = 5;      wymaga   113  obserwacji

• p = 25%;   q = 4;      wymaga    72  obserwacji

Model Helmerta – standaryzacja poprawek (niemianowane zmienne losowe)

Model Helmerta ma zastosowanie w pomiarach rutynowych gdzie znana jest aprioryczna macierz kowariancyjna obserwacji C. Wprowadzamy standaryzację zmiennych losowych – obserwacji i ich poprawek w postaci $V_{i} = \frac{v_{i}}{\sigma_{i}}$

W ten sposób wszystkie poprawki nie mają jednostek miar obserwacji - są niemianowane i porównywalne między sobą. Wyrównanie obserwacji niezależnych sprowadza się do najprostszego wariantu MNK tak jak w przypadku wyrównania obserwacji jednakowo dokładnych. Charakterystyczną cechą modelu Helmerta jest obliczanie wariancji funkcji obserwacji bez współczynnika wariancji ${\hat{\sigma}}_{0}^{2}$.

Wykrywanie błędów grubych (obserwacji odstających)

Obserwacja w toku jest charakteryzowana poprawką v_i (błędem pozornym) stanowiącym estymator błędu prawdziwego ε_i. Jego wartością oczekiwaną, powinno być więc zero E(v_i) = 0. Testowanie wartości uzyskanych poprawek wyrównawczych ${\hat{v}}_{i}$ pozwala na wykrycie błędów grubych poszczególnych obserwacji a nawet wad parametrycznego modelu wyrównawczego.

Test globalny Helmerta (1876)

Pozwala na ocenę całego materiału obserwacyjnego i modelu wyrównawczego.

Taka szansa istnieje tylko w modelu Helmerta i dotyczy tzw. wariancji resztkowej czyli minimalnej wartości funkcji celu ${\hat{\mathbf{v}}}^{\mathbf{T}}\mathbf{C}_{\mathbf{\text{ll}}}^{\mathbf{- 1}}\hat{\mathbf{v}}$ która w założeniu rozkładów normalnych obserwacji ma rozkład χ²_n-u (chi-kwadrat o n-u stopniach swobody) – dla prawdopodobieństwa 95% w Excelu test sprowadza się do nierówności ${\hat{\mathbf{v}}}^{\mathbf{T}}\mathbf{C}_{\mathbf{\text{ll}}}^{\mathbf{- 1}}\hat{\mathbf{v}}$ <Rozkł.Chi.Odwr(0.95,n-u)

W przybliżeniu dla pierwiastka ze współczynnika wariancji (błąd spostrzeżenia typowego o wartości oczekiwanej równej jedności) do nierówności (Nowak 2000)

$\delta_{0} = \sqrt{\frac{{\hat{\mathbf{v}}}^{\mathbf{T}}\mathbf{C}_{\mathbf{\text{ll}}}^{\mathbf{- 1}}\hat{\mathbf{v}}}{n - u}} < 1 + \frac{1.145}{\sqrt{n - u + 0.435}}$

Należy nadmienić, że pierwszy wzór $m_{m} = \frac{m}{\sqrt{2\left( n - u \right)}}$na błąd średni błąd średniego błędu podał Gauss.

Test szczegółowy Baarda’y (1968)

Pozwala na wykrycie (i zlokalizowanie jeśli są spełnione wymogi niezawodności wewnętrznej) błędu grubego obserwacji. Taka szansa istnieje zarówno w modelu Helmerta jak i Gaussa-Markowa. Test dotyczy bezwzględnych wartości poprawek które w założeniu rozkładów normalnych obserwacji mają rozkłady normalne n(0,1) – dla prawdopodobieństwa 95% test sprowadza się do nierówności $\left| {\hat{v}}_{i} \right| < 1.96*\sigma_{v_{i}} = 1.96\sqrt{\left( C_{\hat{v}\hat{v}} \right)_{\text{ii}}}$

W modelu Gaussa-Markowa do obliczenia odchyleń standardowych poprawek wyrównanych σ_vi wykorzystywana jest wartość estymatora współczynnika wariancji ${\hat{\sigma}}_{0}^{2}$

Estymacja przedziałowa

Wyznaczamy przedział który z zadanym prawdopodobieństwem p pokrywa nieznaną (prawdziwą) wartość parametru. Do jego określenia wykorzystamy model statystyczny i wartości obserwacji

• Estymacja wartości oczekiwanej na podstawie niezależnych obserwacji zmiennej losowej X o rozkładzie n(m,σ). Estymatorem wartości oczekiwanej jest średnia $\overset{\overline{}}{x} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}$

  • Znane odchylenie standardowe Estymator określony średnią jest liniową funkcją zmiennych o znanym rozkładzie normalnym – posiada więc również rozkład normalny o tej samej wartości oczekiwanej i odchyleniu standardowym σ/n a po standaryzacji zmienna $u = \frac{\overset{\overline{}}{x} - m}{\sigma/n\ }$ ma rozkład normalny n(m,σ). Dla poziomu ufności $\alpha = \frac{1 - p}{2}$ określimy wartość krytyczną u_α dla której P(u < u_α) = α a ponieważ rozkład normalny jest symetryczny to przedział ufności ma postać $\overset{\overline{}}{x} - \frac{u_{\alpha}\sigma}{n} < m < \overset{\overline{}}{x} + u_{\alpha}\sigma/n$

  • Nieznane odchylenie standardowe wyznaczymy przy pomocy estymatora ${\hat{\sigma}}^{2} = \frac{1}{n - 1}\sum_{i = 1}^{n}\left( x_{i} - \overset{\overline{}}{x} \right)^{2}$ . Zmienna losowa $t = \frac{\overset{\overline{}}{x} - m}{\hat{\sigma}\ }$ ma rozkład t-Studenta o n-1 stopniach swobody – na jego podstawie wyznaczymy wartość krytyczną t_α z zależności P(t<t_α) = α i określimy przedział ufności $\overset{\overline{}}{x} - t_{\alpha}\hat{\sigma} < m < \overset{\overline{}}{x} + t_{\alpha}\hat{\sigma}$

• Estymacja wariancji. Estymatorem wariancji jest ${\hat{\sigma}}^{2} = \frac{1}{n - u}\sum_{i = 1}^{n}\left( x_{i} - \overset{\overline{}}{x} \right)^{2}$, występująca w nim suma kwadratów niezależnych zmiennych losowych o rozkładzie normalnym posiada rozkład χ² o n-u stopniach swobody (podany przez Helmerta w 1876 roku). Rozkład jest niesymetryczny dlatego określimy dwie wartości krytyczne P(χ²<χ₁²)=α i P(χ²> χ₂²)=α określają przedział stąd ostatecznie poszukiwany przedział

• Obszar wahań prostej regresji y|x :   y = ax + b Na podstawie szeregu obserwacji mamy wyznaczone MNK wartości estymatorów współczynników prostej regresji $\hat{a},\hat{b}$ oraz ich macierz kowariancyjną $\mathbf{C} = \begin{bmatrix} C_{\text{aa}} & C_{\text{ab}} \\ C_{\text{ab}} & C_{\text{bb}} \\ \end{bmatrix}$ dla dowolnej wartości x możemy wyznaczyć wartość oczekiwaną funkcji regresji $\hat{y} = \hat{a}x + \hat{b}$ (wartości te tworzą obraz prostej regresji) oraz przedział ufności dla prawdopodobieństwa 90% z rozkładu normalnego $< \hat{y} - 1.96*\sigma_{\hat{y}},\ \hat{y} + 1.96*\sigma_{\hat{y}} >$ gdzie $\sigma_{\hat{y}}^{2} = C_{\text{aa}}x^{2} + 2C_{\text{ab}}x + C_{\text{bb}}$. Granice przedziału ufności tworzą hiperbole wyznaczające obszar wahań prostej regresji. Prosta regresji przechodzi przez środek ciężkości zbioru punktów obserwowanych i w tym punkcie obszar jest najwęższy.

Przykład – wykorzystamy elementy przykładu wyrównania sieci niwelacyjnej z poprzedniego wykładu

Wyrównanie przebiegało według modelu Helmerta – poszczególne obserwacje były traktowane jako niezależne o znanych wartościach odchyleń standardowych σ_i

Czynność 11. Obliczenie ${\hat{\sigma}}_{0}^{2} = \frac{1}{n - u}{\hat{\mathbf{v}}}^{\mathbf{T}}\hat{\mathbf{v}}$=0,96

dostarcza danych do testu globalnego Helmerta ${\hat{\mathbf{v}}}^{\mathbf{T}}\mathbf{C}_{\mathbf{\text{ll}}}^{\mathbf{- 1}}\hat{\mathbf{v}}$ <Rozkł.Chi.Odwr(0.95,n-u)

w naszym wypadku n-u=1 więc test 0,96<3,841 jest spełniony

Czynność 7. Obliczenie macierzy odwrotnej i rozwiązanie układu równań normalnych $\begin{bmatrix} \left( \mathbf{A}^{\mathbf{T}}\mathbf{A} \right)^{\mathbf{- 1}} & \mathbf{x} \\ \end{bmatrix}$

0,008333	0,001667	-0,040
0,001667	0,008333	-0,200

dostarczyła danych do określenia niedokładności wyrównanych współrzędnych

$$\mathbf{\sigma}^{\mathbf{2}}\left( {\hat{\mathbf{x}}}_{\mathbf{i}} \right)\mathbf{=}{\left( \mathbf{A}^{\mathbf{T}}\mathbf{A} \right)^{\mathbf{- 1}}}_{\mathbf{\text{ii}}}$$

stąd $\mathbf{\sigma}^{\mathbf{2}}\left( {\hat{\mathbf{H}}}_{\mathbf{11}} \right)\mathbf{= 0,008333}$ $\mathbf{\sigma}\left( {\hat{\mathbf{H}}}_{\mathbf{11}} \right)\mathbf{= 0,091}$

$\mathbf{\sigma}^{\mathbf{2}}\left( {\hat{\mathbf{H}}}_{\mathbf{12}} \right)\mathbf{= 0,008333}$ $\mathbf{\sigma}\left( {\hat{\mathbf{H}}}_{\mathbf{12}} \right)\mathbf{= 0,091}$

Czynność 9. Obliczenie $\begin{bmatrix} {\mathbf{A}\left( \mathbf{A}^{\mathbf{T}}\mathbf{A} \right)}^{\mathbf{- 1}} & \mathbf{\text{Ax}} \\ \end{bmatrix}\mathbf{= A}\begin{bmatrix} \left( \mathbf{A}^{\mathbf{T}}\mathbf{A} \right)^{\mathbf{- 1}} & \mathbf{x} \\ \end{bmatrix}$ oraz A_i(A^TA)⁻¹A_i^T

0,083333	0,016667	-0,40	0,833333
-0,03333	0,033333	-0,80	0,333333
-0,01667	-0,08333	2,00	0,833333

dostarcza danych do obliczenia wariancji poprawek

$\mathbf{\sigma}^{\mathbf{2}}\left( {\hat{\mathbf{v}}}_{\mathbf{i}} \right)\mathbf{=}{\mathbf{1 -}\mathbf{A}_{\mathbf{i}}\left( \mathbf{A}^{\mathbf{T}}\mathbf{A} \right)}^{\mathbf{- 1}}\mathbf{A}_{\mathbf{i}}^{\mathbf{T}}$

W naszym wypadku dla wszystkich obserwacji

$\mathbf{\sigma}^{\mathbf{2}}\left( {\hat{\mathbf{v}}}_{\mathbf{i}} \right)\mathbf{= 0,166667}$ oraz $\mathbf{\sigma}\left( {\hat{\mathbf{v}}}_{\mathbf{i}} \right)\mathbf{= 0,408}$

stąd testy szczegółowe Baarda’y dla wszystkich obserwacji

Czynność 10. Obliczenie obserwacji wyrównanych, kontrola ostateczna,

obliczenie odchyleń standardowych poprawek, testy szczegółowe

Lp	P	K	ΔΗ	σ	Hp-Hk	V	v	obs+v	Hk-Hp	Kontrola	σV	|V|/σV	<1,96
1	101	11	-0,32	0,1	-0,32	0,40	0,04	-0,28	-0,28	0,00	0,408	0,98	tak
2	11	12	0,04	0,2	0,04	0,80	0,16	0,20	0,20	0,00	0,816	1,96	NIE
3	12	102	0,15	0,1	0,39	0,40	0,04	0,19	0,19	0,00	0,408	0,98	tak

Wyrównanie sieci niwelacyjnej algorytmem Banachiewicza z kontrolami i testami

Sieć niwelacyjną tworzą punkty charakteryzowane wysokościami H_i oraz pomierzone przewyższenia H_jk = H_k − H_j. Parametrami są wysokości punktów nowo wyznaczanych. Równania obserwacyjne sieci niwelacyjnej są liniowe. Punkty nawiązania (repery) nie podlegają wyrównaniu, są stałymi a nie wyznaczanymi niewiadomymi. Przykładowo - ciąg niwelacyjny o punktach wyznaczanych 11 i 12, nawiązany do reperów H₁₀₁ = 10.111 i H₁₀₂ = 10.222 składa się z przewyższeń H_101, 11 = −0.10;  H_11, 12 = −0, 01;   H_12, 102 = 0.25 zaobserwowanych z jednakową dokładnością σ_H = 0.01. Przyjmując wysokości przybliżone H₁₁^przybl = H₁₂^przybl = 10.000 otrzymamy układ równań poprawek

Oznaczenie przewyższenia	Równanie poprawki	Standaryzowane równanie poprawki
	dh 11	dh 12
101-11	1	0
11-12	-1	1
12-102	0	-1

Równania normalne z kontrolą sumową

dh11	dh12	L	S	Kontrola = suma-S
20000	-10000	-210	9790	0,00
-10000	20000	380	10380	0,00
-210	380	10,05	180,05	0,00

Rozszerzona tabela wtórna Choleskiego-Banachiewicza

dh11	dh12	L	S	kontr	dh11	dh12	101-11	11-12	12-102	s	kontr
141,42	-70,711	-1,485	69,226	0,000	0,007	0,000	0,707	-0,707	0,000	138,459	0,000
	122,474	2,245	124,720	0,000	0,004	0,008	0,408	0,408	-0,816	249,452	0,000
		2,803	2,803	0,000	0,001	-0,018	-0,967	-0,967	-0,967	2,690	0,000
						v	-0,010	-0,010	-0,010	Otrębski
				kwadr	6,67E5	6,667E5	0,66667	0,66667	0,666667	2
				błędy	0,0082	0,0082	0,8165	0,8165	0,8165

Wysokości wyrównane

Nazwa punktu	H^przybl	dH	H^wyr
11	10,000	0,001	10,001
12	10,000	-0,018	9,982

Test globalny SPEŁNIONY ponieważ

$$\delta_{0} = \sqrt{\frac{{\hat{\mathbf{v}}}^{\mathbf{T}}\mathbf{C}_{\mathbf{\text{ll}}}^{\mathbf{- 1}}\hat{\mathbf{v}}}{n - u}} = \sqrt{\frac{2,803}{3 - 2}} = 1,67 < 1 + \frac{1.145}{\sqrt{n - u + 0.435}} = 1,96$$

Kontrola ostateczna	Testy szczegółowe
Oznaczenie	L^wyr=Hk-Hp
101-11	-0,110
11-12	-0,019
12-102	0,240