9. fejezet - A nem lineáris direkt feladat esete

Tartalom

9.1. Gravitációs mérés

Az előző fejezetben tárgyalt példákban (egyenes és másodfokú függvény illesztése) a direkt feladatok a paraméterek lineáris függvényei voltak. Ez számos esetben teljesül, például valamely mérési pontokra történő polinomfüggvények illesztésekor. Az általános geofizika területén számos olyan eset van, amikor egy függvényt sorfejtés alakjában írjuk fel, ekkor a függvény a sorfejtés együtthatóinak lineáris függvénye, ezért a sorfejtés együtthatóinak a mérési adatokból történő meghatározása az előző fejezetben tárgyalt módszer segítségével elvégezhető.

Az ilyen típúsú direkt feladat nagy előnyt jelentett az A alakmátrix számításához, mivel az alakmátrix a direkt feladatnak a paraméterek szerinti deriválásával jön létre, ezért ha a direkt feladat a paraméterek lineáris függvénye, akkor az alakmátrix nem tartalmazza a paramétert. Ez azt jelenti, hogy az alakmátrix független a paramétertől, így annak számításához nem kell ismernünk (vagy feltételeznünk) a paraméter (előzetes) értékét.

Ez a fenti feltétel már néhány viszonylag egyszerű esetben sem teljesül.

Tekintsük az alábbi példát! Illesszünk sinus függvényt, a 9.1. táblázatban található x-y pontpárokra!

x

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

y

0

2

3

6

9

11

12

11

9

6

2

0

2

5

8

10

11

10

9.1. táblázat. x-y pontpárok.

A direkt feladat az alábbi alakú lehet:

(9.1)

A paraméterek közül a0 egy konstans y irányú tolást jelent, a1 a sinusos jel amplitúdóját, a2 a körhullámszámot, a3 pedig a jel fázisát. Képezzük (9.1) egyenletnek a mérési helyeken vett deriválásával az alakmátrixot!

                                                                                                                               

(9.2)

Jól látszik, hogy az alakmátrixban szereplő mennyiségek a paraméterek függvényei.

A problémát úgy orvosolhatjuk, hogy a paramétereknek kezdőértéket (előzetes értéket vagy kiinduló értéket) adunk. Ezek a kezdőértékek lehetőség szerint a tényleges megoldáshoz közelítő értékek legyenek.

A sinus függvény illesztés esetén, ha a pontokra kézzel sinusgörbét rajzolnánk, meghatározhatnánk a paraméterek kezdőértékeit.  Leolvasható lenne az illesztett függvény átlagértéke, vagyis az a0 paraméter kezdőértéke. A rajzolt görbe amplitúdóját tekinthetjük az a1 paraméter kezdőértékének. 2π-t elosztva a leolvasott hullámhosszal kapjuk az a2 paraméter kezdőértékét. Megbecsüljük végül, hogy a függvényt vízszintes irányban a hullámhossz hányad részével kellene eltolni, hogy a sinushullám kezdőpontja (az inflexiós pontja) az x = 0 pont fölé essen, és ebből kiszámoljuk az a3 paraméter kezdőértékét. A továbbiakban a paraméterek kezdőértékeit a felső indexbe, zárójelbe írt 0-val jelöljük. Ezzel a jelöléssel a paramétervektor kezdőértéke: x(0).

Ezeknek a mennyiségeknek a segítségével már ki tudnánk számolni az A alakmátrixot.

Az előző fejezetben említettük, hogy az l tisztatag vektort úgy definiálhatjuk, hogy az a mért értékekből kivonjuk a kezdő paraméterekkel (előzetes értékekkel) számított direkt feladat megoldást. Ezt a definíciót most is alkalmazzuk.

(9.3)

Alkalmazzuk ezt a felírást a jelen feladatra! A direkt feladatot egy f(...) függvénnyel jelöljük. Keressük azt az x értéket (a paraméterek javítás vektorát) amelyre igaz az, hogy az értékeit hozzáadva a paraméterek kezdő értékeihez (x(0)+x), az f(x(0)+x) függvény értékei a méréseink kiegyenlített mennyiségeit szolgáltatják (U). Ekkor a direkt feladat megoldás a legjobban megközelíti a méréseinket (L). A mérésekhez hozzáadva a javítások vektorát (v) kapjuk a mérések kiegyenlített értékeit (U).

(9.4)

A függvény kiegyenlített értékét írjuk fel a függvény Taylor-sorfejtése segítségével:     

(9.5)

A függvény első deriváltja az A Jacobi mátrix. Ennek felhasználásával – kifejezve a javítások vektorát – kapjuk, hogy

(9.6)

ahol az l tisztatag vektor alakja:

(9.7)

Látható, hogy a kapott (9.6) egyenlet formailag megegyezik a lineáris esettel.

Visszatérve a sinusfüggvény illesztési feladathoz, a paraméterek kezdő értékeivel kiszámolt A alakmátrix és l tisztatag vektor segítségével a (8.29) képlettel kiszámoljuk az x vektor. Az x vektor itt nem a paraméterek konkrét értékét jelenti, hanem a paraméterek javítását. A paraméterek javított értéke az x(0) vektor és az x vektorok összege.

Amennyiben tehát a direkt feladat a paramétereknek nem a lineáris függvénye, úgy az Taylor-sor első két tagja nem adja meg pontosan a függvény értékét az (x(0)+x) helyen, ennek következtében a kapott x vektorral javított paraméterek még nem adják meg az eltérés négyzetek minimumát.

Ennek a problémának a megoldására iterációt végzünk. Ennek során az első iterációs lépésben a fenti műveletekkel felírjuk az A alakmátrixot és az l tisztatag vektort, amiből kiszámoljuk paraméterek x javításvektorát. A második iterációs lépésben a megjavított  (x(0)+x) paramétervektort tekintjük a paraméterek kezdőértékének, ezzel számoljuk ki az A alakmátrixot és az l tisztatag vektort, majd ezekből a paraméterek javítását leíró x vektort. Ezt a lépés többször megismételve általában azt látjuk, hogy a javítások x vektora egyre kisebb abszolút értékű számokat tartalmaz, az újabb lépésekkel a kiegyenlített paramétervektor elemei konkrét értékekhez konvergálnak. Az iterációt általában néhány lépés után leállítjuk. A leállítás feltétele általában az, hogy a kapott javításvektor normája az újabb iterációs lépés hatására már kevesebbet változik, mint egy előre megadott kicsiny küszöbérték.

Probléma lehet, hogy a normál egyenlet csak azt írja elő, hogy a függvény szélsőértékét keressük, ahol a függvény első deriváltja nulla. Az, hogy egy minimumhely, vagyis itt a reziduálok négyzetösszege minimális, a függvény második deriváltjának vizsgálatából derül csak ki.

A minimalizálandó függvényünk tehát a reziduálok négyzetösszege. Ez egy skalár mennyiség. Esetünkben ezt kritériumfüggvénynek tekinthetjük. (A kritériumfüggvény általánosabb definícióját az inverzió statisztikus elméleténél tárgyaltuk.) A kritériumfüggvény által visszaadott skalár mennyiség, valóban a paraméterek függvénye, és mi azt a paraméterkombinációt (paraméter vektor elemeinek azon értékeit) keressük, ahol ez a függvény minimális. Ezt a függvényt jelöljük g-vel.

(9.8)

Itt xi –k a paramétervektor elemeit – nem pedig a javításukat – jelenti. Legyen a g(x) függvény kétszer deriválható! Képezzük a g(x) függvénynek a paraméterek szerinti második deriváltjait tartalmazó mátrixot!

(9.9)

A H mátrix neve Hesse mátrix. A parciális deriválás sorrendjének felcserélhetősége miatt a Hesse mátrix mindig szimmetrikus. A Hesse mátrix vizsgálata mutatja meg, hogy az adott függvényérték minimum, maximum, vagy nyeregpont.

A minimum feltétele, hogy a Hesse mátrix pozitív definit legyen.

A pozitív definitség több megfogalmazása közül kiemelünk két definíciót. Az első szerint egy mátrix pozitív definit, ha valamennyi sajátértéke pozitív. Ezzel egyenértékű definíció, hogy valamennyi, a mátrix bal felső sarkából képzett négyzetes mátrix determinánsa pozitív legyen. (Konkrét számítás során bármelyik definíció felhasználható a pozitiv definitség vizsgálatára.)

A nemlineáris direktfeladat esetén láttuk, hogy a reziduálok négyzetösszegének minimumához tartozó paraméterértékek megtalálása csak több iterációs lépésben lehetséges. Ennek az az oka, hogy a direkt feladat függvényt a kiinduló paraméterek környezetében egy Taylor-sor első tagjáig bezárólag írtuk fel.

A kritériumfüggvényt felírhatjuk a kiinduló paraméterek környezetében a Taylor-sor második tagjáig bezárólag:

(9.10)

A fenti képletben a J Jacobi-mátrix a g skalár függvény gradiense (jelen esetben egy vektor), a H mátrix pedig az előbb megismert Hesse mátrix. Mivel ezt az egyenletet a kritériumfüggvényre írtuk fel, ennek a függvénynek a paraméterek (az x vektor elemei) szerinti minimalizálása adja a keresett paraméterek javításait. Az egyenletben szereplő J vektor és H mátrix a gyakorlati problémák esetén általában numerikus deriválással határozhatók meg a paraméterek kiinduló értékeinél (x(0)). Az egyenlet megoldása nagy mátrixok esetén általában a „Minimum kereső eljárások” c. fejezetben bemutatott iterációs eljárások valamelyikével történik.

9.1. Gravitációs mérés

Tekintsük az alábbi gravitációs mérési problémát, amelynek valamilyen változatát több inverziós könyv is felhasznál az inverzió alapelveinek bemutatására.

Egy megadott sík terület egy pontján egy nagytömegű gömb van elásva. Graviméterrel mérjük meg a nehézségi gyorsulást a területen egy (közel) egyenközű rácsban, és ebből határozzuk meg a gömb tömegét, felszín alatti mélységét és az elásás helyének síkkordinátáit!

Tételezzük fel, hogy egy jólképzett geofizikus elvégezte a mérést, tehát nekünk nem kell foglalkoznunk a mérés megtervezésével és végrehajtásával, rendelkezésünkre állnak a mérési helykoordináták és a graviméterrel mért nehézségi gyorsulás értékek.

Irányítsuk a koordinátarendszerünket úgy, hogy az X tengely észak felé, az Y tengely kelet felé, a Z tengely pedig lefelé irányuljon. A felszín a Z = 0 magasságban legyen. A felszín egy pontjának koordinátája így: (x, y, 0). A nehézségi gyorsulás háttér étéke – vagyis amit a gömbtől olyan távolságban mérnénk, ahol a gömb hatása már elhanyagolható – legyen g0 . A területen mért nehézségi gyorsulás értékekből vonjuk ki a nehézségi gyorsulás háttér értékét, az így kapott anomáliaértékeket tekinthetjük a gömb hatásának.  A mért adatainkat a (9.2.) táblázat foglalja össze.

Mérés sorszáma

X[m]

Y[m]

Z[m]

Δg[m/s2]

1

100

300

0

-1,68922e-014

2

100

310

0

2,08384e-013

3

100

320

0

9,11552e-014

4

100

330

0

1,3398e-013

5

100

340

0

1,35904e-013

6

100

350

0

7,48595e-014

7

110

300

0

1,07663e-014

8

110

310

0

1,71658e-013

9

110

320

0

3,90526e-013

10

110

330

0

9,64405e-013

11

110

340

0

6,26942e-013

12

110

350

0

3,54347e-013

13

120

300

0

-7,3977e-014

14

120

310

0

2,87387e-013

15

120

320

0

6,87852e-013

16

120

330

0

2,7552e-012

17

120

340

0

1,52322e-012

18

120

350

0

2,83205e-013

19

130

300

0

1,87324e-013

20

130

310

0

3,69917e-013

21

130

320

0

4,54991e-013

22

130

330

0

1,26672e-012

23

130

340

0

1,08389e-012

24

130

350

0

4,87435e-013

25

140

300

0

-1,34739e-013

26

140

310

0

-3,68395e-015

27

140

320

0

1,58061e-013

28

140

330

0

4,50285e-013

29

140

340

0

1,37541e-013

30

140

350

0

1,24567e-013

31

150

300

0

1,61628e-013

32

150

310

0

6,26579e-014

33

150

320

0

1,30561e-013

34

150

330

0

7,56786e-014

35

150

340

0

1,80268e-013

36

150

350

0

1,66113e-013

9.2. táblázat. Gravitációs anomáliaértékek (). Az anomáliaértékek szimulált adatok, azokat egy ismert helyű és tömegű ható terének számításával, majd a számított anomáliaértékekhez véletlen jellegű zaj hozzáadásával kaptuk.

A fenti táblázat adatait egy kétváltozós függvény értékeinek tekinthetők. Ennek vizualizációját – a mérési pontokra illeszkedő szintfelületet – az alábbi ábrán láthatjuk.

9.1 ábra. A „mért” gravitációs anomáliák () megjelenítése szintfelületként.

A 9.1 ábrán egy aszimmetrikus anomália csúcs látszik. Ebből következtethetünk, hogy a ható nem pontosan a rácspont alatt helyezkedik el.

Az inverzió végrehajtására, először írjuk fel a direkt feladatot! Egy tetszőleges felszíni, (x, y, 0) koordinátájú pontban az elásott gömb által okozott nehézségi gyorsulás (az ún. anomália), az alábbi képlettel számítható:

(9.11)

A képletben szereplő G, a gravitációs állandó, x és y a mérési pont koordinátái. A képletben szereplő m, x0, y0, z0 étékek azonban a feladat kiírásában szereplő meghatározandó ismeretlenek, vagyis a paraméterek.

A fenti képlet tehát azt fejezi ki, hogy ha ismernénk a paraméterek m,x0,y0,z0, értékeit, ki tudnánk számolni, hogy egy tetszőleges (x,y,0) koordinátájú helyen a gömb mekkora anomáliát okozna.

Az A alakmátrix meghatározásához tehát a (9.11) egyenletnek a mérési helyeken a paraméterek szerinti deriváltjait kell képezni. Esetünkben ezek a deriváltak analitikusan is számolhatók:

(9.12, 9.13, 9.14, 9.15)

Látszik, hogy ezekben a deriváltakban szerepelnek a paraméterek, tehát ahhoz hogy ezeket a deriváltakat számítani tudjuk, a paramétereknek valamilyen hihető kiinduló értékeket kell adnunk.

A paraméterek kiinduló értékeit tehát nekünk kell megadnunk. Válasszuk az anomália-térkép maximumhelyét a sikkordináták előzetes értékének: és  . A mélység előzetes értékét az anomáliakép félérték-szélességéből becsülhetjük: . A gömb tömegének kiinduló értékét válasszuk az 1 méter sugarú ólomgömb tömegének:  értéke: . (Az ólomgömb tömegének számításánál az ólom és a felszínközeli kőzetek sűrűségkülönbségét használtuk.)

A felvett előzetes értékek alapján már a mérési helyekre ki tudjuk számolni az A alakmátrixot. A mátrix 4 oszlopból (a paraméterek száma) és 36 sorból (a méréseink száma) áll. Az l tisztatag vektort a 36 pontban mért anomáliaértékek (9.2. táblázat) és felvett előzetes paraméterértékek alapján kiszámolt elméleti értékek (9.11 képlet) különbségéből képezzük.

A paraméterek javítását a

(9.16)

alapján számoljuk. A számolt javítások:

(9.17)

A javításokat hozzáadva az előzetes értékekhez megkapjuk a paraméterek kiegyenlített értékét. Újabb iterációs lépést végezve, a paramétereink tovább javulnak. Néhány iterációs lépés után a paramétereink már keveset változnak, ekkor leállítjuk az iterációt. A paraméterek becsült értékei:

(9.18)

A paraméterek kovariancia mátrixát a

(9.19)

képlet segítségével számolhatjuk, ahol N a mérések, P pedig a paraméterek száma. A számításhoz a PLL mátrixot egységmátrixnak tételeztük fel.

(9.20)

amiből a becsült paraméterek szórásai:

(9.21)

A gravitációs mérés feldolgozásához használt adatok nem valódi gravitációs adatok voltak, hanem szimulált adatok. Ezeket az adatok úgy álltak elő, hogy egy rögzített paramétervektor segítségével kiszámoltuk az anomáliaértékeket a feltételezett mérési pontokban, és az anomáliaértékeket 10-10 szórású, normális eloszlású, korrelálatlan zajjal szórtuk meg. A paraméterek rögzített értékei: x0 = 121, y0 = 133, z0 = 10 és m = 4733 kg.

Jelen esetben a paraméterek becsült értékei a „valódi” értékeket az egyszeres szórásnál jobban megközelítik. Amennyiben a szimulált mérésekhez nagyobb hibát adtunk volna, úgy a paraméterek szórásai is jelentősen megnövekedtek volna.

Ebben a feladatban az anomáliaértékeket – gravitációs gyakorlattól eltérően – SI mértékegységben adtuk meg. A geofizikában szokásos mértékegységek esetén a maximális hatás m/s2 nagyságrendbe esett, ami közel egyharmad mikrogal. Ez a legmodernebb – a szupravezetés elvén működő – graviméterekkel már jól mérhető pontosság.