-
Verfahren zum Bestimmen der Intervallgeschwindigkeit von reflektierenden
Zonen in unterirdischen Formationen.
-
Die Erfindung bezieht sich auf geophysikalische «rforschung bzw.
Untersuchungen und insbesondere auf ein Verfahren zum Umwandeln von Feldseismogrammen
zu Intervallgeschwindigkeit und Neigung (dig) jeder unterirdischen Schicht.
-
bei seismischen Untersuchungen ist die Messung der treschwindigkeit
allgemein als Hauptparameter beim Verarbeiten und Interpretieren von Seismogrammen
anerkannt. Die bestimmung der Charakteristiken der Schallgeschwindigkeit von Seismogrammen
ist beschrieben in "Seismic Velocities From Subsurface Measurements", O.il. Dix,
Geophysics, Vol. 20, Seiten 68 - 86, 1965.
-
Kurz gesagt, wird die Schallgeschwindigkeit von Seismogrammen üblicherweise
bestimmt durch Verwendung des Verhältnisses
In dieser Gleichung bedeutet TR die Laufzeit einer seismischen Welle, die von einer
Quelle zu einer Zwischenfläche und zurück zu einem Oberflächenempfänger wandert.
Dies ist die Zeit der Reflektion auf dem Seismogram. 20 ist die ideale senkrechte
Wanderzeit von der Erdoberfläche zur Reflektionsstelle. (Diese Zeit wird auch als
Reflektionszeit bei einer
waagerechten Versetzung von Null bezeichnet).
x ist die waagerechte Versetzungsstrecke von der Quelle zu dem mpfänger, und Va
ist die Schallgeschwindigkeit. Während üblicherweise die vezeichnung "mittlere Geschwindigkeit"
oder Durchschnittsgeschwindigkeit2 verwendet wird, ist es genauer, festzustellen,
was die RMS-Schallgeschwindigkeit ist.
-
Hierzu wird auf Dynamics Correlation Analysis" verwiesen, eine Veröffentlichung
von W.A. Schneider und £ilo backus, dargeboten auf der 36-ten annual SSS Convention,
Houston, Texas, 1966.
-
zinke Technik, die neuerdings verwendet worden ist, um mittlere Geschwindigkeit
von einer folge von Seismogrammen zu erhalten, ist in der USA-Patentschrift 3 417
370 beschrieben. bei einer Technik wie dieser umfassen die Felddaten eine folge
von Seismogrammen, die entlang einer Untersuchungslinie einen waagerechten Abstand
bzw. eine waagerechte Versetzung mit bezug aufeinander haben. Diese Spuren sind
mit bezug aufeinander zeitverschoben, basierend auf verseniedenen angenommenen Werten
der mittleren Geschwindigkeit oder Durchschnittsgeschwindigkeit. i'ür jede reflektion
zeigt eine Korrelationsoperation die Signalenergie zwischen den Spuren für die verschiedenen
Werte angenommener Durchschnittsgeschwindigkeit an.
-
bin Maximum in der Signalenergie zeigt richtige burchscnnittsgeschwindigkeit
oder mittlere Geschwindigkeit für diese fteflektion an. Dies kann für jede reflektion
auf der folge von Seismogrammen ausgeführt werden, so dass die mittlere Geschwindigkeit
zu jedem dieser Reflektoren bestimmt wird.
-
In der schwebenden Patentanmeldung US Serial Nr. 769 590 vom 22.
Oktober 1968 wird eine Verbesserung dieser Art von Arbeitsweise beschrieben, wobei
eine Mehrzahl von Neigungswerten (dip) für jeden Reflektor ausgewertet wird. Das
heisst, für jeden Reflektor erfolgt ein Durchprüfen aller möglichen Neigungswerte
und Geschwindigkeitswerte, um die maximale Signalenergie zwischen den Spuren zu
bestimmen. Dies schafft Ausgänge, welche die Neigung jedes Reflektors und die mittlere
Geschwindigkeit
zu jedem Reflektor darstellen. Arbeitsweisen wie diese werden nachstehend als VIP-Verfahren
bezeichnet, ein Akronym für Geschwindigkeit durch integrierte energie binde andere
Verarbeitungstechnik für Seismogramme umfasst sianderung oder Verlagerung jedes
Reflektionspunktes zu seiner korrekten Tiefe und waagerechten Stellung. Gegenwärtig
verwendete Verlagerungstechniken werden benutzt, um Refraktion oder brechung des
Strahlenpfades in geschichteter Erde zu berücksichtigen, beispiele von aJanderungstechniken
oder Verlagerungstechniken (migration) sind in "Exploration Geophysics, beschrieben,
veröffentlicht von Trija Publishing Company, zweite Ausgabe, 1950, Los Angeles,
Kalifornien, J.J. Jakosky, Seiten 670 und 688. begenwärtig angewendete Verlagerungsverfahren
verwenden eine Geschwindigkeit, die sich nur mit der Tiefe ändert. Solche Verfahren
behandeln Geschwindigkeitsgraaienten, wie sie bei divergierenden Schichten auftreten,
nicht richtig.
-
bei dem Verisj-ren gemäss der -rfindung werden Intervallgeschwindigkeiten
durch ein Iterationsverfahren bestimmt, durch welches ein willkürliche Neigungen
zeigendes Schichten modell nach dem Snell'schen Gesetz an die telddaten angepasst
wird.
-
Reflektionen auf den beismogrammen sind durch drei Parameter gekennzeictmet,
nämlich durch ihre Ankunftszeit T0 bei waagerechter Versetzung Null, ihre scheinbare
mittlere Geschwindigkeit Va und durch ihre Zeitneigung (time dip) #t.
-
Das heisst, ein Satz von Reflektionen 1, 2, 3 ... i kann oharakterisiert
werden durch den batz von Ankunftszeiten Ti(X) und Neigungen #Ti. Diese Parameter
können direkt von den xeldseismogrammen bestimmt werden, und zwar durch mehrere
Verfahren, jedoch sind die vorgenannten VIP-Verfahren für diesen Zweck besonders
geeignet. bei der Erfindung wird die Intervallgeschwindigkeit erzeugt durch Anpassung
eines Schrägschichtenmodells an diese Daten.
-
Zusätzlich wird das verlagerte unterirdische Verhalten jedes Reflektors
erhalten0 Die Verwendung des Weschwindigkeitsmodells ermöglicht das Verlagern der
ursprünglichen Daten in ihre richtige Kaumstellung.
-
Die Erfindung wird nachstehend an Hand der Zeichnung beispielsweise
erläutert.
-
zeigt 1 ist ein Fliessdiagramm des Verfahrens gemäss der Erfindung.
-
Fig. 2 ist eine Feldaufstellung, die zum Erhalten der Felddaten geeignet
ist.
-
Fig. 3a zeigt StrahlendiagrammeO Fig. 3b zeigt Strahlendiagramme.
-
Fig 4 zeigt den Vergleich beobachteter Laufzeiten zu berechneten
Laufzeiten für einen Geschwindigkeitswert.
-
Fig. 5 zeigt die Fehlerfunktion als Funktion der Geschwindigkeit
Fig. 6 zeigt eine Folge von Felddiagrammen.
-
Fig. 7 zeigt Daten, die mittels eines VIP-Verfahrens erhalten worden
sind.
-
Fig. 8 zeigt Daten, die von oder mittels eines VIP-Verfahrens erhalten
worden sind und die eine Neigungsuntersuchung enthalten.
-
Fig. 9 zeigt den Ausgang des Verfahrens gemäss der ttfindung.
-
Bevor mit einer Beschreibung des Verfahrens gemäss der Erfindung
in Verbindung mit dem Fliessdiagramm gemäss Fig. 1 fortgefahren wird, wird zunächst
beschrieben, was bei dem Verfahren ausgeführt wird Fig. 2 zeigt eine typische ieldaufstellung
zum erzeugen einer Folge von beismogrammen. Gemäss iig. 2 wird seismische Energie
an den Schusspunkten S1, S2 und 53 erzeugt Die seismischen Spuren werden an den
Empfängern X1, k2 und R3 aufgezeichnet. Gewöhnlich sind selbstver-' ständlich viele
Empfänger vorhanden, um eine Folge von Seismogrammen zu erzeugen, wobei die empfänger
mit bezug aufeinander
entlang einer Untersuchungslinie eine waagerechte
Versetzung bzw. einen waagerechten Abstand haben. Die besondere Art von Feldaufstellung
gemäss Fig. 2 ist ein Aufzeichnungssystem mit gemeinsamem Oberflächenpunkt. Die
seismischen Reflektionen, die auf der Folge von Seismogrammen erscheinen, können
durch hyperbolische Kurven quer über die Solge angenähert werden, definiert durch
X2 go2 - R2 2 0 - R V a In der vorstehenden Formel ist TR die Ankunftszeit der reflektion
auf einer besonderen Spur in der Seismogrammfolge, X die waagerechte Versetzung
der besonderen Spur, 20 die Reflektionszeit bei waagerecht er Versetzung Null, und
Va die mittlere oder anscheinende Geschwindigkeit.
-
Es gibt zahlreiche Techniken, um anscheinende mittlere Geschwindigkeit
und die Reflektionszeit bei waagerechter Versetzung Null zu erhalten. Diese beiden
Parameter werden benötigt, um an der Folge von Seismogrammen eine Korrektur hinsichtlich
des Unterschiedes der normalen Ankunftszeit (normal moveout correction) vorzunehmen.
Das VIP-Verfahren zum Erhalten der anscheinenden Geschwindigkeit und der Reflektionezeit
bei waagerechter Versetzung Null wird nachstehend an Hand der Fig. 7 und 8 kurz
beschrieben.
-
bine Folge von Spuren, die in ihrem Aussehen ähnlich der Spur gemäss
Fig. 6 sind, ist von einer Feldaufstellung erhalten. (Fig. 6 enthält eine Anzahl
von Spuren die jeweils waagerecht voneinander in einem Abstand liegen. Die Ordinate
gemäss Fig. 6 ist die Aufzeichnungszeit von 0 bis etwa 3,0 sekunden; Die Abszisse
stellt die waagerechte Versetzung entlang der Untersuchungslinie dar.) bei jedem
Inkrement der Aufzeichnungszeit werden die @ Spuren mit Bezug aufeinander zeitverschoben.
Diese Zeitverschiebung (bzw. Korrektur hinsichtlich des normalen Unterschieds
in
der Ankunftszeit) basiert auf einem angenommenen Geschindigkeitswert. Dann wird
die Signalenergis aller zeitverschobener Spuren für diesen angenommenen Geschwindigkeitswert
erhalten. Eine tíecnnik, um dieses auszuführen, besteht darin, die Quer-Korrelationsfunktion
aller Spuren für eine Verzögerung oder Nacheilung von Null zu erhalten. Der Wert
dieser Quer-Korrelationsfunktion für Verzögerung oder Nacheilung Null wird für jede
angenommene Geschwindigkeit erhalten. sAie in lig. 7 dargestellt, ist der Wert der
mittleren Geschwindigkeit über eine Mehrzahl von Werten zwischen 5.000 und 13.000
Fuss je Sekunde iteriert worden. itig. 7 zeigt eine Mellrzahl von Signalenergiekurven,
und zwar eine Signalenergiekurve für jede einer Mehrzahl von Aufzeichnungszeiten.
bpitzen in der Signalenergiekurve zeigen den richtigen wert der angenommenen scheinbaren
Geschwindigkeit für einen besonderen Reflektor an. Anhand Fig. 7 ist zu bemerken,
dass verschiedene Reflektoren vorhanden sind, die eine scheinbare Geschwindigkeit
haben, die gemäss der Darstellung zwischen 5.000 und 7.000 Fuss je Sekunde liegt.
-
Die einzelne Kurve gemäss Fig. 7 ist eine sarstellung des Maximuns
in den Signalenergiekurven. Dies ist eine gute Anzeige der Reflektionszeit bei waagerechter
Versetzung Null einer Mehrzahl von reflektierenden Schichten in den unterirdischen
Formationen.
-
Zusätzlich zum Durchprüfen aller Geschwindigkeitswerte und senkrechten
Wanderungszeiten, wie es durch die Ausgangsdaten gemäss iig. 7 angezeigt ist, können
mit einem Verfahren dieser art auch eine heihe angenommener Neigungswerte durchgeprüft
werden, und zwar sowohl positive als auch negative. Diese Neigungskomponenten können
entweder in Form des angenommenen Neigungswinkels., d.h. verschiedener Ainkelwerte,
oder als Zeitverschiebung zwischen der ersten und der letzten Spur in der Folge
ausgedruckt werden Fig. 8 ist eine Darstellung der Gesamtenergiekurven, deren jede
der einzelnen Kurve gemäss Pig. 7 ähnlich ist, und
zwar für jeden
einer Mehrzahl angenommener Neigungswerte. In diesem Fall ist die Neigung ausgedrückt
als Gesamtzeitverschiebung zwischen benachbarten Punkten gemeinsamer Tiefe.
-
Aus dieser Information kann die annähernde Neigung jeder seismischen
neflektion bestimmt werden.
-
Während der Ausgang eines VIP-Verfahrens in den Fig. 7 und 8 zeichnerisch
dargestellt ist, ist zu verstehen, dass die Ausgänge gewöhnlich in digitale i'orm
gebracht werden und verfügbar sind, urn als eingang für ein Verfahren wie das Verfahren
gemäss der Erfindung verwendet zu werden. Von einem Verfahren der oben beschriebenen
Art sind seismische Daten verfügbar, welche die Ankunftszeiten T0 bei waagerechter
ersetzung Null von Reflektionen, die scheinbare Geschwindigzeit Va für jeQe -eflektion
und die Zeitneigung #tn jeder seismischen Reflektion umfassen. Aus diesen Daten
erzeugt das Verfahren gemäss der erfindung die Intervallgeschwindig-Keit für jede
Schicht und eine genauere bestimmung der Neigung. Zusätzlich wird die verlagerte
waagerechte und senkrechte Versc@iebung jedes heflektors erhalten. Dies wird erzielt
durch Erzeugung eines Modells des Strahlenpfades, genommen durch: die aufeinanderfolgenden
Schichten, und aurch Vergleich von ihm mit den beobachteten Felddaten, um das beste
Modell zu erhalten. Dies wird Schicht für Schicht ausgeführt, bs wird zunächst die
erste Schicht in Verbindung mit den Fig. 3a und 3b betrachtet. Die IntervalIgeschwindigkeit
V1 der ersten nicht wird genau bestimmt aus der Gleichung V1 =Va cos Die mittlere
Geschwindigkeit Va zu dem ersten Reflektor ist aus dem Ausgang des vorhergehenden
VIP-Verfahrens bekannt0 Dieses Verfahren spezifiziert gleichfalls die Neigung #l
(oder #t) für die erste Schicht. Aus der Neigung und der ersten Schicht wird die
Richtung des Strahlenpfades in der ersten Schicht bestimmt. Die Zeit T0 für waagerechte
Versetzung
Null aus dem VIP-Verfahren bestimmt die Dicke der Schicht.
-
Wenn diese Informationen für die erste Schicht direkt erhalten sind,
nimmt das Verfahren nunmehr seinen Fortgang, um die Informationen für die nachfolgenden
Schichten durch ein Modellverfahren zu erhalten. Zuerst werden verschiedene Geschwindigkeitswerte
für die zweite Schicht angenommen, bas heisst, die Geschwindigkeit wird aufeinanderfolgend
iteriert zu V21, V22, V23 usw. Für jede Geschwindigkeit bestimmt das Snell'sche
gesetz die Strahlenrichtung in der zweiten Schicht aus dem Ausdruck sinα1
= sinα2 V1 V2 -us der Strahlenrichtung in der zweiten Schicht kann die Neigung
der zweiten Schicht berechnet werden, weil der Strahl auf die zweite Schicht normal
auftrifft. Die Zeit 'i für waagerechte Versetzung Null aus den @elddaten bestimmt
zusammen mit der iterierten Intervallgeschwindigkeit, wie dick die zweite Schicht
ist. Nunmenr können die Laufzeiten für jeden wert von X für jeden iterierten Geschwindigkeitswert
bereciiiiet werden.
-
Zum Beispiel sind für die Geschwindigkeit von V21 die berechneten
Reflektionszeiten für jeden Wert waagerechter Versetzung durch die ausgezogene Linie
in Fig. 4 dargestellt0 Diese sind aus duu bekannten Migrationsgleichungen berechnet;
d.h. für jede Schicht kann die waagerechte Verschiebung des Strahles zwischen dem
oberen xtlue des Intervalls und dem unteren @nde des Intervalls bestimmt werden.
Die Länge des Strahlenpfades iii jeder Schicht wird berechnet und die Laufzeit über
die Schicht wird bestimmt unter Verwendung des angenommenen wertes deer Intervallgeschwindigkeit.
-
Die berechneten Wert von laufzeiten TNi für verschiedene werte von
A sinu in Fig. 4 durch die ausgezogene Linie dargestellt. Die beobachteten Werte
von Reflektionszeit an jeder der Spuren sind durch die Kreise in Fig. 4 angedeutet.
-
Die bumme der Unterschiede zwischen den berechneten und den beobachteten
Werten von Reflektionszeit für jeden ert von X werden summiert, um eine Fenlerfunktion
zu erzeugen. Das heisst, der fehler ist
worin i der Spurenindex, m die Anzahl von Spuren in der Folge, TNi die beobachtete
Zeit der Reflektion N auf der Spur i und TNi die berechnete Zeit der Reflektion
N auf der Spur i ist. Dieser Fehler wird für jede der iterierten Geschwindigkeitswerte
bestimmt. ns wird ein-e Reihe von Werten der Fehlerfunktion für jeden der iterierten
Geschwindigkeitswerte erzeugt. Die Fehlerfunktion als eine Funktion der Geschwindigkeit
ist in Fig. 5 dargestellt. Die den minimalen Fehler erzeugende Geschwindigkeit ist
die richtige Intervallgeschwindigkeit. Wie nachfolgend erläutert, wird das Minimum
durch Newton'sche Iteration ausgewählt. Das Minimum bestimmt die Intervallgeschwindigkeit
für die Schicht.
-
Nachdem die Intervallgeschwindigkeit V2 für die zweite Schicht bestimmt
ist, wird das Verfahren für die dritte Schicht wiederholt. Zusammenfassend sind
nunmehr unter Bezugnahme auf Fig. 3b die ueschwindigkeiten V1 und V2 genau bekannte
Die Neigung der ersten und der zweiten Schicht ist nunmehr genau bekannt. Die Zeit
T0 bei waagerechter Versetzung Null für die dritte schicht zeigt die Dicke der dritten
Schicht an. Für einen angenommenen Geschwindigkeitswert ist die Strahlenrichtung
in der dritten Schicht gegeben durch das Snell'sche Gesetz durch sinα4 = sinα3
V3/V2 Nunmehr können die iciufzeiten TNi durch die dritte Schicht für jeden Wert
von X berechnet werden. Die Berechnung
der Laufzeit durch die dritte
nicht kann in Übereinstimmung mit der nachstehenden Gleichung erfolgen
worin D. die Länge des Strahles in der i-ten Schicht, und V. die Intervallgeschwindigkeit
in der i-ten Schicht isto bei der errechnung von TNi werden die werte der Intervallgeschwindigkeit
und der Neigung e verwenuet. Die Leufzeiten für verschiedene Werte von X werden
fLir jeden einer Mehrzahl iterierter Geschwindigkeitswerte V31, V32, V33 usw.
-
berechnet.
-
Die berechneten Werte von TNi werden mit beobachteten inerten von
TNi verglichen, um eine Fehlerfunktion zu erzeugen. Das Minimum in dieser Fehlerfunktion
wählt oder bestimmt den neuen Wert der Intervallgeschwindigkeito An dieser Stelle
ist zu bemerken, dass das tatsächliche nachstehend zu beschreibende Verfahren eine
Variation des vorstehenden Verfahrens ist. Insbesondere werden die beobachteten
Reflektionszeiten TNi aus dem Seismogramm nicht verwendet, obwohl dies die ideale
Technik sein würde. Stattdessen wird ihr hquivalent
verwendet. Dieser Ausdruck ist für jeden Wert von X als Ausgang des vorgenannten
VIP-Verfahrens verfügbar. bs ist ersichtlicn, dass der Ausdruck zur Reflektionszeit
durch folgende Gleichung in beziehung gesetzt werden kann
Jedoch ist es, anstatt alle Reflektionen von der iolge von Seismogrammen aufzunehmen,
leichter, die mittlere Geschwindigkeit
zu verwenden, die aus dem
VIP-Verfahren verfügbar ist.
-
Die Erfindung kann anhand des Blockdiagramms gemäss Fig. 1 besser
verstanden werden. Der Eingang des Verfahrens, bei 1 angedeutet, umfasst Reflektionszeiten
TNi, worin N die Nummer der Schicht und i der Superindex sind. Wie oben erwähnt,
kann an Stelle der Reflektionszeiten der Eingang des Verfahrens die mittlere Geschwindigkeit
aus dem VIP-Verfahren sein. Andere Eingänge zu Fig. 1 umfassen #tN, die Neigung
aus @em VIP-Verfahren, THRES, der Schwellenwert für die @ehlerfunktion und VMIN,
den minimalen Wert der Intervallgeschwindigkeit, der als gültig angenommen wird0
@ie bei 7 angedeutet, wird ein Zähler schrittweise fortgeschaltet, so dass die Intervallgeschwindigkeiten
und Neigungen aufeinanderfolgend für jede Schicht bestimmt werden.
-
@ie oben erwähnt, ist es sehr wichtig, dass eine genaue Bestimmun,
der Intervallgeschwindigkeit für die erste Schicht erhalten vird. ie bei 3 angedeutet,
wird die Intervallgeschwindigkeit für die erste Schicht aus der bekannten mittleren
@eschwindigkeit und der Neigung der ersten Schicht bestimmt.
-
ie bei 4 angedeutet, wird ein Modell erzeugt, wobei eine nachfolgende
Schicht dem Modell während jedes Ansführens einer Schleife ninzugefügt wird. Für
jeden angenommenen wert der Intervallgeschwindigkeit V bestimmt das Snell'sche esetz
die Strahlenrichtung in der Daten Schicht durch den Ausdruck V2@ sinα2N =
sinα (2N-1).
-
V(2N-1) hierdurch wird die Neigung der Daten Schicht bestimmt, weil
der Strahlenpfad normal bzw. senkrecht auf die ll-te Schicht auftrifft Wenn das
Modell abwärts bis zur IJ-ten Schicht vervollständigt ist, ist das folgende bekannt:
Die Zeit T0 bei waagerechter
Versetzung Null zu der N-ten Schicht,
welche die Dicke der Schicht bestimmt. Weiterhin sind die Neigungen der Schicht
N-1 und aller vorhergehender Schichten bekannt.
-
Ausserdem sind die Intervallgeschwindigkeit VN der vorhergehenden
Schicht und aller weiteren vorhergehenden Schicht ten bekannt. Auch ist die Neigung
der l4-ten Schicht aus der berechnung nach dem Snell'schen Gesetz bekannt. Nunmehr
können die Laufzeiten TNi berechnet werden, wie es bei 5 angedeutet ist.
-
A Die berechneten Werte TNi werden von den entsprechenden Werten
TNi aus den Felddaten abgezogen, um eine Fehlerfunktion zu erzeugen. VN wird über
aufeinanderfolgende Werte iteriert, um eine Fehlerkurve zu erzeugen, und das Minimum
in dieser Kurve wird durch Newton'sche Iterationstechnik festgestellt, wie es bei
6 angedeutet ist. Wenn der minimale fehler festgestellt ist, wird er mit dem Schwellenwert
THRES verglichen, um zu bestimmen, ob er niedriger als der ochwellenwert ist. uies
ist bei 7 angedeutet. Wenn der fehler kleiner als der Schwellenwert ist, wird die
dem minimalen Fehler zu@eordnete Geschwindigkeit #N mit der minimalen annehmbaren
Geschwindigkeit verglichen, wie es bei 8 angedeutet ist. enn #N grösser als VMIN
ist, dann wird #N zusammen mit der Neigung der N-ten Schicht #@, wie sie aus den
Berechnungen nach dem Snell'schen gesetz bestimmt ist, gespeichert. Die Speicherung
ist bei 9 angedeutet0 In diesem hall wird der Zänler 2 einen Schritt fortgeschaltet
una die gleichen Funktionen werden für eine nachfolende Schicht ausgeführt0 Wenn
der Fehlerschwellenvergleich gemäss Schritt 7 oder der Minimalgeschwindigkeitsvergleich
gemäss Schritt 8 anzeigt, dass die bestimmung der Geschwindigkeit für diese Schicht
nicht zufriedenstellend gemacht ist (d.h., dass eine falsche Anzeige aus den Schritten
8 oder erhalten ist), wird diese Schicht fallengelassen. Dies ist bei 10 angedeutet.
Wenn beispielsweise die bingangsdaten aus dem VIP-Verfahren angezeigt haben, dass
vier Schichten 1, 2,
3, 4 vorhanden sind und wenn eine falsche
Anzeige in der Berechnung für die dritte Schicht vorhanden ist, wird diese Schicht
fallengelassen und die Schichten werden so numeriert, dass der letzte oder endgültige
Ausgang Geschwindigkeiten für die behichten 1, 2, 3 umfasst.
-
Aus Vorstehendem ist ersichtlich, dass das Verfahren gemäss der brfindung
unter Verwendung mehrerer bekannter Arten von Berechnungsvorrichtungen verwirklicht
werden kann.
-
Das Verfahren ist besonders geeignet zur Verwendung mit einem allgemeinen
ZwecKen dienenden Digitalrechner. Während die erfindung mittels verschiedener Programme
ausgeführt werden kann, wird ein geeignetes Programm, dargeboten in FORTRAN-Ausdrücken,
nachstehend angegeben.
-
PROGRAM STRATV (IM@UT,OU@@UT,PLOT,TA@@28=PLOT) DISEMSION TT (60,60),AT
(60),T(60),TIM(60),V(60),YAXIS(200),SCRAT@ 1000),D(60),VRAR(60),K@(60),XDIST(60),ZDIST(60),ALPHA(60),@D@L
(6@@ DIMEMSION VV(60).RAA@G(50),VVV(60).TT(60) DIMENSION XRAY (60).YRAY (60) ,DIPRAY
(60),DEPTH(60) DIMEMSION EE (60) .FV(60) DIMENSION IAREA (@),ISHOTP(9),TLINE (9)
EQUIVALENCE (V (1),VV (2)) EQUIVALENCE (SCRAT (1999),TIXL) SCRAT (2000) = 1.0 5
READ 10,NTR,LAYERS.VZFRO.DELR,@ZERO,EOL,ERROR,NPW,SPRD 10 FORMAT (2@10'5@@0.5,I5,F5.0)
PI=3.14159265 @ PII=PI/2.
-
15 TOL = (EOL**MP@) *NTR 16 READ 29,(TI(I),I=I,LAYERS) 17 READ 25.(VVV
(I)'I=1,LAYERS) 18 READ 29,(TDEL (I),I=I,LAYERS) 25 FORMAT (8F10.0) 26 NLAYERS=LAYERS
VVZERO=VZERO LAYEPS=NLAYERS VZERO=VVZERO DO 27 I=1,NTR 27 XDIST (I) = ((I - 1)*DELH)
/2.0 290 DO 24 II = I,LAYERS DO 24 I = 1,NTR TT (II,I) = TI (TI) * TI (II) *4.0*XDIST
(I)*XDIST (I)/(VVV(II)*VVV (II)) TT (TT.I) = SORTE (TT(II.I)) 24 CONTINUE 28 DO
281 I=1.LAYERS 281 AT (I)=TT (I,1) D(1)=SPRD RAANG (I) = TDEL (1) ALPHA (I) = PII
DO 333 I = 2,LAYERS RAANG (I) = TDEL (I) D(I) = 40000.
-
333 ALPHA (I) = PII V(LAYERS+1) = 15000.0 @ D (LAYERS+1)=40000.0 @
ALPHA(LAYERS+1) = PII 2290 NLAY=LAYERS LDRP=0 II = 1 T (1) = TT (1,1) /2.0 DO 30
I = 2.LAYERS 30 T (I) = (TT(I,I) - TT (I-1,1)) /2.0 31 IF (II.GT.LAYERS) 100,35
35 CONTINUE.
-
34 IF (FLAG.FQ.1.0) 36,39 36 SCRAT (2000)=1. 0 D (I)=SPRD.
-
RAANG (I)=TDEI (I) 39 CONTINUE DO 40 I=1.MTR 40 AT (I)=TT (TT . T)
HZFRO=SP@D 45 CALL VELFST (AT.T,TIM,V,VZERO, TOL,HZERO,DELH,NTR.ERROR,II,NPM.FLA@
1,XDIST,NCASE,ALPHA.D,SCRAT,VV,RAANG,SE,SV)
53 COMT@@@E IF (@L
AG.FQ@ I) 55,59 55 KK (LD@@+1)=I@+LDPR EE (LD@@+1)=5P EV (LDRP+1)=5V LAYERS=NLAY-(LDRP+1)
LD@P=LDRP@@ IF (KK(LD@@).@@.NLAY)100,54 54 DO 57 J=TI.L@YERS RAAN@(J)=@AA@G(J+I)
DO 57 JJ=1.MTR 57 TT (J.JJ)=TT (J+I,JJ) DO 5@ J=II.LAYERS 58 T (J) = (TT (J.1)-TT
(J-1,1)) /2.0 GO 10 35 59 COMTIMUE 659 XRAY (II)=SCRAT(1799+II) 60 YRAY (II)=SCRAT(1859+II)
65 DIPRAY(II)=(SCRAT(1919+II)/PI)*180.0 70 VZERO=V(II)-V(II)*.13 75 II=II=1 80 GO
TO 31 100 CONTINUE 110 PRINT 700, (XRAY(I),YRAY(I),DIPRAY(I),V(I),I=1,NLAYER) 700
FORMAT(4E16-8) END SUBROUTINF VELEST (TT,T,TIM ,V,VZERO,TOI ,HZERO,DELH,NTR.ERROR.LAYE
IS.@PW,FLAG.XDIST,NCASF,ALPHA.D,SCRAT,VD.RAANG,EE,VV) DIMENSION A@@HA (1).D(1) DIMENSION
TT (1).T(1),TIM (1),V (1).SF (5).SV (5),SCRAT (2000),XDIS@ (1) DIMENSION VD (1),PAANE
(1) 5 VMIN=4500.0 FLAG=0 JIM (1) =IT(1) ITERS=10 V (LAYERS)=VZFRO 10DO 30 II=1.2
CALL EPRVEL (J,V.II@@CELH.NTR.ITERS.ERROR,LAYERS,NL,HZERO,TT,SE(II) 1,XDIST,NCASE.ALPHA,D,SCRAT.VD.RAANG)
21
CONTINUE SV (II) = V (LAYERS) 20 V(LAYERS) =VZERO+0.001@VZERO 30 CONTINUE 40 SLOPE
= (SE (1) - SE (2))/(SV (1) - SV (2)) COR=SE (2) /SLOPE V (LAYERS) =V(LAYERS)-COR
IF (V(LAYERS).LT.0.0) 55,50 50 CALL ERRVE@ (T,V,TIM,CELH,NTR,ITERS,ERROR,LAYERS,NL,HZERO,TT,EE
1,XDIST,NCASE,ALPHA,D,SCRAT,VD,RAANG) 61 CONTINUE VV=V (LAYERS) 60 IF (EE.LT.SE
(2)) 66,55 66 CONTINUE SE (1) = SE (2) $ SV (1) = SV (2) SE (2) = EE $ SV (2) =
V (LAYERS) 67 GO TO 40 55 DUME=MAX@F (SE (1).SE (2)) IF (SF (1).FQ.DUME) 80,75 75
DUMV=SV (2) SE (2) = SE (1) $ SV (2) =SV(1) SE (1) =DUM@ $ SV (1) = DLMV 80 V (LAYERS)
= SV (2) @SIGNF ((0.0001 *SV(2)),SLOPE) CALL ERRVEL (T,V,TIM,@ELH,NTR,ITERS,ERROR,LAYERS,NL,HZERO,TT,EE
1,XDIST,NCASE,ALPHA,D,SCRAT,VD,RAANG) 81 COMTINUE VV = V (LAYERS) 85 IF (EE,LT.SE
(2)) 95,90 90 SLOPE = -SLOPE 95 IF (SV (2).LT.VMIN.AND,SLOPE.GT.0.0) 96,97 96 FLAG
= 1.0 GO TO 640 97 DELV=500,0 100 I@C=1 V (LAYEPS) = SV (2) 110 V(LAYEPS) =V (LAYERS)
-SIGNF (DELV,SLOPE) CALL ERRVEL (T.V,TIM,DELH.NIR,ITERS,ERROR,LAYERS,NL,HZERO,TT,SE(IN@
1+2),XDIST,NCASE,ALPHA.D,SCRAT,VK,RAANG) 111 CONTIN@E VV=V (LAYERS) 120 IF (SE(INC
+ 2).GT.SE (2).AND.INC.LT.2) 125.130 125 DELV=DELV/2.0 IF (DELV.LE.1.0) 126.100
126 EE=SE (3) SV (3) =VV GO TO 220 130 SV (INC+2) =V (LAYERS) INC=INC+1 140 IF (INC.GT.2)
145.110 145 IF (SE (3).LT.SE (2).AND.SE(4).LT.SE (3)) 150,160 150 IF (SV (2).LT.VMIN)
151,155 151 FLAG=1.0 GO @0 640
155 SE (2)=SE(3) $ SV(2)=SV (3)
SE (3)=SE (4) $ SV(3)=SV (4) INC=2 GO TO 110 160 IF (ARSF (SV (3)-SV (4)).LE.5.0)
161,165 161 IF (SE (3).LF.TOL) 166.163 166 FE=SF (3) $ V(LAYERS)=SV (3) GO TO 750
163 FLAG=1.0 GO TO 640 165 CALL FITPARB (SV (2).SV (3).SV(4),SE(2),SE(3),SE(4),VV,EE)
V (LAYERS)=VV CALL ERRVEL (T,V,@IM,DELH,NTR.ITERS,ERROR.LAYERS,NL,HZERO,TT,EE 1,XDIST,NCASE,AL@HA,D,SCRAT,VD,PAANG)
175 DUME=MINIF (EF,SE (3)) 180 IF (EF.EQ.DUME) 190.185 185 SE (2)=EE $ SV(2)=VV
GO TO 195 190 SE (2)=SE (3) $ SV (2)=SV (3) SE (3)=EE $ SV (3)=VV 195 VABS=ABSF
(SV (3)-SV (2)) IF (VABS.LE.5.0) 220.197 197 V (LAYERS)=SV (3) + (SV(3)-SV (2))
SV (4)=V (LAYERS) CALL FRRVEL (T.V,TIM,DELH,NTR,ITERS,ERROR,LAYERS,NL,HZERO,TT,EE
1,XDIST.NCASE.ALPHA,D,SCRAT,VD,RAANG) 198 CONTINUE SE (4)=EE 200 IF (EE.LT.SF (3))
205.210 205 VV=V (LAYERS) DUME=EE GO TO 190 210 CALL FITPARB (SV (2),SV (3),SV (4),SE
(2),SE (3),SE (4),VV,EE) V (LAYERS)=VV 215 CALL ERRVEL (T,V,TIM,CELH,NTR,ITERS,ERROR,LAYERS,NL,HZERO,TT,EE
1,XDIST.NCASE,ALPHA,D,SCRAT,VD,PAANG) 225 DIFF=ABSF (SV (3)-VV) IF (DIFF.LF.5.0)
220.175 220 IF (EF.LE.TOL) 639,163 639 IF (VV.LE.VMIN) 96,166 640 COMTINUE 750 RETURN
END
SUBRQUII@F @RRVEL (T,V,TI@@DEL@,NT@,ITERS,ERROR,LAYERS,NL,HZERO,TI
1,EE.XDIST,MCASE.AL@HA,D,SCRAT,VV,RAANG) DIMENSION ALPHA (1).D(1),SCRAT (1),VV(1),RAANG
(1) DIMENSION IT (1).T (1).TIM(1).V (1).XDIST (1) CALL TIMFCV(T.V,T)@,XDIST,NTR,ITFRS,FRROR.LAYERS.NL,HZERO,NCASE,@
IPHA,D,SCRAT,VV,RAANG) EE=0.0 DO 10 I=1,@TP 10 EE=EE + ABSF (TT (I) - TIM (I)) RETURN
END SUBROUTINF TIMECV (T,V,TIM,H,NTP,TIERS,ERROR,LAYERS,NL,HZERO,NCAS@, 1ALPHA,D,SCRAT,VV,RAANG)
DIMENSION ALPHA (1),D (1),SCRAT (1) DIMENSION T (1),V (1),TIM (1),H (1),VV (60),RAANG
(1) DIMENSION T20UM (100) KLI=LAYERS ZOT=0.0 DO 10 T=1.KLI ZOT=ZOT+T (I) 10 T2DU@
(I) = 2.0 * ZOT 25 CALL TWOTPHI (VV,D,ALPHA,H,NTR,ERROR,KLI,TIM,SCRAT,SCRAT (200),
SCRAT 1 (400),SCRAT (600), RAANG ,SCRAT (1000),SCRAT (1200),SCRAT (1400),SC 2RAT
(1600),T2DUM,SCRAT (1800),SCRAT (1860),SCRAT (1920)) 30 RETURN END
SU@R@U@@UF
FTTPAR@ (X1,X2,X3,Y1,Y2,Y3,ARK,A@P) 5 A1=X1*X1 B1=X1 C1=1.0 D1=Y1 A2=X2*X2 B2=X2
C2=1.0 D2=Y2 A3=X3*X3 B3=X3 C3=1.0 D3=Y3 10 A2@=A1*A2 B2F=A1*B2 C2E=A1*C2 D2E=A1*D2
AlF=A2*A1 BlF=A2*B1 ClE=A2*C1 DlF=A2*D1 A4=A1E-A2E B4=B1E-B2E C4=C1E-C2E D4=D1E-D2E
A2F=A3*A2 15 B2F=A3*B2 C2E=A3*C2 D2E=A3*D2 A3F=A2*A3 B3F=A2*B3 C3E=A2*C3 D3F=A2*D3
A5=A2E-A3E B5=B2E-B3E C5=C2E-C3E D5=D2E-D3E H4E=B5*B4 C4F=B5*C4 20 D4F=B5*D4 B5E=B4*B5
C5E=B4*C5 D5E=B4*D5 B6=B4E-B5E C6=C4E-C5E D6=D4E-D5E CN5=D6/C6 BN5=(D4-C'@S*C4)/B4
AN5=(D1-CN5*C1-BN5*B1)/A1 DERV=-2.0*ANS APK=BMS/DERV AMP=APK*AP@*ANS + APK*ENS+CNS
RETUR@ END
SU@POUTI@@ TWOT@HI@V,@@ANG,SPDIST,@TR,@RROR,KLI,TIME,X,Y,ALPRA.
-
1TWI,PM.TW,VVV.ZZZ,T@@@,REFT,XRAY,YPAY,ALRAY) DIMEMSIOM XRAY (1),YRAY
(1),ALPAY (1) DIMENSIOM V (1),Z (1),ANG (1),SRDIS@(1),TIME (1), X (1),Y (1),ALPHA
(1), 1T@@ (1),P@ (1),TM (1),VVV (1),ZZZ(1),TTTM (1) 2,REFT (1).TPE@ (60),@P@ (60)
IF (X(2000).EO.1.0) 3115,3116 3115 CO@TINUE X (2000)=0.0 DO 3114 T=1,60 3114 TDEL
(1)=PM (I) SPED=Z (1) Z (1)=0 3116 CONTINUE PI=3.14159265 $ PII=PI/2.
-
NNNL=KLI+1 3113 PPM(KLI)=ASIMF(TDEL (KLI)@.5/SPRD*V (2)) 3200 COMTINUE
DO 3210 I=1,NNNL C PRINT 10,(REFI(I),I=1.NNNL) C PRINT 10,(PRM(1),I=1,NNNL) C PRINT
10,(V(1),I=1.NNNL) C PRINT 10,(TTTM (I),1=1,NNNL) C PRINT 10,(ANG (I),I=1,NNNL)
C PRINT 10,(Z (I),I=1,NNNL) 3210 TTTM (I)=ANG (I) 10 FORMAT (6E16.8) Z (KLI + 1)=40000.
$ ANG (KLI + 1)=PII CALL TILT (REFT,PRM,V,TTTM,Z,X,Y,ANG,NNNL,KLI) ALRAY (KLI)=ANG
(NNML)-PII C IMAGE SYSTEM FOR PRIMARIES IS COMPUTED.XXXXXXXXXXXXXXXXXXXXXXXXXXXXX
CALL IMAGE (V.Z.ANG,KLI,VVV,ZZZ,TTTM,ALLPHA,DZ) NLI=2@KLI+1 NSE=0 SH=0.
-
H=0.
-
PM (1)=TTTM (KLI+1) -PII IKEY=0 DANG=PI/72.
-
DO 1000 K=1,NTR IF (NSE) 902,@02,903 902 CONTINUE IF (IKEY.GF.1)
908.911 908 DO 909 I=2.NLI DIP=PII-TTTM (I) ZZZ (I)=ZZZ (I)-TAMF (DIR)*TRASP 909
COMTINUE 903 COMTINUE 911 COMTINUE ITER=2 DELH=0.0 60 COMTINUE DO 30 L=1,ITER C
PAY PATHS APE TRACFD.XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX NR=1 DPHI=PM
(1) D@@=DELH
C PRINT 601@@@@(1) 601 FORMAT (4H @AY5[@]6.@) CALL
TILTTW1 (TW1,P@,VVV,TTTM,ZZZ,X,Y,ALPHA,NL1,NP,TM,ALL@@A,@Z, 1 DELH,INTRS) TIXL=X(1999)
TIXL=TIXL+1.@ IF (TIXL.GT.4000.@)602,603 602 X(1999)=-5.0 GO TO 2000 603 CONTINUE
X(1999)=TIXL IF(H.EQ.0.0)5000,599 5000 YRAY (KL.@)=X(KLI+1) YRAY (KLI)=Y(KLI+1)
599 CONTINUE HDH=ABS@ (DELH-H) PRINT 20,HD@,PM(1),DDH,@@@@,H,ERROR 20 FORMAT (1H
6F16.8) IF (HDH-ERROR)50,50,600 600 CONTINUE PM(1)=PM(1)+DANG 30 CONTIN@E IF (TIEP-2)80,70,80
80 CONTINUE C PRINT 20.@@(1),DX,DANG C GRADIENT IS
[email protected])
DX=H-DEL@ DANG=PI/180. $ ITER=2 DANG=0.0001 GO TO 100 70 CONTINUE ITER=1 GRADH=(DELH-DDH)/DANG
DX=H-DELH XDX=DX/GRADH IF (XDX.GT.PII)4041.4042 4041 XDX=DANG 4042 CONTINUE C PRINT
20,GRADH,XDX,P@(1),DANG,DX PM(1)=XDX+DPHI DANG=0 100 CONTINUE HDH=ABSF(DELH-H) IF
(HDH.LT.ERROR)50,60 50 CONTINUE DANG=0.0001 C @TERATION IS COMPLETE AND DESIRED
RAY IS AT HAND, TWO WAY TIME IS JJ.XXX TIME(K)=TWT(@) H=S@DIST(K+1)*2.0 TRASP=-(SRDIST(K+1)=SRDIST(K))
IKEY=1 1000 CONTINUE 3010 CONTINUE 2000 CONTINUE RET@RN END
SU@POUTI@E
IMAGE (V,Z,TMM,NL,VV,ZZ.TTM.ALPHA,DZ) C V=LAYER VFLOCITY. Z=IMTEOFACF DEPIH. TMM=INTFRFACE
ANG@E C@@ DOW C NL=NUMBER OF LAYERS SYSIFM, NLIS=NUMBER LAYERS IMAGE SYS.
-
C OUTPUT YTELDS PARAZET@RS OF TMAGE SYSTEM.
-
C AN N LAYER SYSTEM IS AN N+1 INIERFACE SYSTEM C INTEPFACE (1) IS
SUPFACE, V(1)=Z(1)=TMM(1)=0.
-
DIMEMSION V(1),Z (1),TMM(1),VV(1),ZZ(1),TTM (1) V (1)=Z (1)=0.0 PI=3.14159265
PII=P1/2.
-
NLI=NL + 1 DO 10 I=1,N@@ VV (I)=V (I) $ JJ=2*ML-T+3 $ VV (JJ)=V (I)
ZZ (I)=Z (I) $ TTM (I)=TTM (I) 10 CONTINUE DO 15 I=1.NL ALPHA=(PII-TMM (NLT)) $
GAMMA=2.*ALPHA + TMM (NLI - I) TTM(MLT+I)=PT-GAMMA DZ=Z(NLI)-Z(NLI-I) DZZ=DZ*STNF
(TMM (NLI - I))/SINF (GAMMA) ZZ (NLI + 1) =Z (NLI) +DZZ 15 CONTIMUE RETURN END
SUBROUTIMF
TTLTTWI@R,PM,V,TMM,Z,X,Y,ALPHA,NL,NR,TM,ALLPHA,DDZZ, 1 DISP,IMTRS) C (X,Y) RAY IMTERCBRT
WITH INTER@ACES C R(J) TTME TO INIFR@ACF (J) C ALPHA (J) ANGLE OF RAY FROM VERTICAL.
AT INTERFACE (J) C NR=NUMREP OF TRIALS TO ESTIMATE DELTA (H)/DELTA (P) C NL=NUMBER
OF LAYERS IM SYSTEM + IMAGE C V=LAYED VFLOCITY, IMM=LAYER DIR, Z=INTERFACE DEPTH
DIMENSION P(1),PM(1),V(1),TMM(1),ALPHA(1),Z(1),TM(1),X(1),Y(1) KREFRAC=0 PI=3,14159265
PII=PI/2.
-
INTRS=0 DO 700 J=1.NR TM (2)=PM (J) TOT=TM(1)=TMM(1)=@ (1)=Y (1)=Z
(1)=0.0 DO 500 I=2,M@ TAU=TMM (I)-TM (I) TAAU=PII=TAU ARG=(V(I+1)/V(I) *SINF (TAAU))
ARSARG=ARSF (ARG) IF (ARSARG-1.) 73.74.74 73 TAUU=ASINF (ARG) TM (I+1)=TAUU + TMM
(1) -PII GO TO 75 74 KREFRAC=1 75 ANGLF=TMM (I) -PII DZ=Z (I)-Y(I-1)-X (I-1)*TANF
(ANGLE) IZ=DZ IF (IZ) 300.300,400 C 300 INTRS=1 $ GO TO 800 C IF (INTRS.EQ.1) 2000.599
300 COMTINUE 400 CONTINUE DIST=DZ*SIMF (TMM (I))/SINE (TMM(I)-TM (I)) TIME=DIST/V
(I) DELX=DIST*SINF (TM (I)) $ DELY=DIST*COSF (TM(1)) X (I)=DELX+X (I-1) $ Y (I)=DELY+Y
(I-1) LAST=I DELH=-Z (LAST)*SINF (ALLPHA)/SINF (PI-ALLPHA-TMM(LAST)) TOT=TOT+TIME
R (J)=TOT ALPHA (J) = (TM (I)) 500 COMTINUE 700 COMTINUE Y LAST =Y (LAST)-Z (LAST)
DISP=DELR-SORTF (Y LAST *Y LAST +X (LAST) *X (LAST)) *SIGNF (1.,X (LAST) C PRINT20,DISP.TOT,DEL@
20 FORMAT (1@ 4HDISP/10E20.4) 800 CONTIMUE RETURN END
SUBROUTIME
TILT(R,@M,V,TMM,Z,HORZ.DEPTH,ALPHA,@L,NR) DIMENSIONR (1),PM (1),V (1),TMM (1),HORZ
(1),DEPTH (1),ALPHA (1),Z(1), 1 TM (100),X (100),Y (100) PI=3,14159265 PII=PI/2.
-
ALPHA (1)=PII C DO 700 J=1.NR J=NR TM (2)=PM (J) TOT=TM (1)=TMM (1)=X
(1)=Y (1)=Z (1)=0.0 DO 500 I=2,NL TAU=TMM (I)-TM (I) TAAU=PII-TAU 75 ARG=(V (I+
1)/V (I)) * SINF (TAAU) IF (ABSF (ARG).GE.1.0) 76.73 76 V (I+1)=V (I+1)-V (I+1)
*.1 GO TO 75 73 TAUU=ASTNF (ARG) TM (I+1)=TAUU+TMM (I) -PII ANGLF=TMM (I) -PII DZ=Z
(I) -Y (I-1)-X (I-1) * TANF (ANGLE) DIST=DZ*SINF (TMM (I) @/SINF (TMM (1)-TM (I))
TIME=2,0*DIST/V (I) DELX=DIST*SINF (TM(I)) $ DELY=DIST*COSF (TM (I)) X(I)=DELX+X
(I-1) $ Y (I)=DELY+Y(I-1) TOT=TOT+TIME DELT=R (J)-TOT IF (DELT) 100,100,500 100
DD=DFLT/TIME DX=DD*DELX DY=DD*DELY DEPTH (J)=Y (I) +DY $ HORZ(J)=X(I)+DX ALPHA (J+1)=TM
(I) + PII Z(J+1)=DF@TH(J)-HORZ (J)/TANF (ALPHA (J+1)) GO TO 700 500 CONTINUE 700
CONTINUE RETURN
Min besonderes Rechnersystem, das zur Verwendung
geeignet ist, wird von "Control Data Corporation unter der allgemeinen Modellbezeichnung
6600 geliefert und es umfasst die folgenden Bauteile: 6604 Zentral Computer, 65MK
Memory 6608 Disc System 6602 Console Display 6681 Data Channel Converter 3228 Magnetic
Tape Controller 607 Magnetic Tape Transport 3447 Card Rader Controller 405 Card
Rader 3256 Line Printer Controller, und 501 Line Printer.
-
Das besondere FORTRAN-Programm zum Ausführen des Verfahrens gemäss
der erfindung einschliesslich gewisser Modifikationen, wird nachstehend angegeben,
wonach eine kurze Beschreibung des Arbeitens des Programms erfolgt. Das Programm
ist in FORTRAN-Sprache angegeben, die für die meisten Digitalrechner geeignet ist.
Für besseres Verständnis der Verwendung von FORTRAN-Ausdrücken wird bezug genommen
auf "Introduction to FORTRAN" von S.O. Plumb, McGraw-Hill Bock Company, ew York,
New York.
-
Das vorgenannte Programm beginnt mit mehreren Dimensionsinstruktionen,
welche Speicherreihen für die zu erzeugenden Parameter reservieren.
-
Die Instruktion bzw. der Befehl 5 liest die Eingangsparameter, die
wie folgt sind: NTR = Nummer der Spuren IAYERS= Nummer der Schichten VZERO = Erste
Geschwindigkeitsschätzung DELH = Spurenabstand HZERO = Versetzung zur ersten Spur
(nicht verwendet) B0L = Fehler je Spur in Sec.
-
ERROR = Konvergenzfehler beim Modell nach dem Snell'schen Gesetz
NPW
= 1 = Absolutwert-tehlerschwelle 2 = Mittleres Quadrat - tehlerschwelle SPRD = Spurenabstand
unterirdisch Die Instruktion 15 stellt den Schwellenwert in Übereinstimmung mit
dem Fehler je spur in Sekunden, dem absoluten Wert der Fehlerschwelle oder die Fehlerschwelle
gemäss dem mittleren Quadrat und die Anzahl der Spuren ein Die Instruktionen 16,
17 und 18 lesen die Schichtenparameter von Bändern ab. Diese Schichtenparameter
sind DI(I), das die Zeit 20 der Schicht I ist, VVV(I), welches die mittlere Geschwindigkeit
zur Schicht I ist, und TDEL(I), welches die Zeitneigung der Schicht I ist. Diese
Parameter werden alle in dem zuvor erwähnten VIP-Verfahren gemessen und sind als
Eingangsparameter verfügbar.
-
Die Instruktionen 26 bis 30 ieiten das Programm ein.
-
Die Instruktionen 34 bis 39 leiten das Programm neu ein, wenn eine
Schicht fallengelassen wird0 Die Subroutine VELEST, die an der Instruktion 45 abgerufen
wird, findet das Minimum in der Fehlerkurve. FLAG wird auf 1,0 eingestellt, wenn
die Schicht nicht angenommen wird0 Die Instruktionen 53 bis 58 numerieren die schichten
neu, wenn eine Schicht fallengelassen wird0 In der Instruktion 65 wandelt DIP RAY(II)
den Neigungswinkel der Schicht II von bogenmass in Grad um.
-
Die Instruktionen 659 und 60 speichern XRAY(II) bzw.
-
YRAY (II), welches die verlagerten Koordinaten in Fuss des Reflektors
sind.
-
Die Intervallgeschwindigkeit der Schicht II ist in V(II) durch die
Subroutine VELEST gespeichert worden0 Bei der Instruktion 70 wird die erste Geschwindigkeitsschätzung
VZERO eingestellt auf einen neuen Wert von VZERO.
-
Die Instruktion 75 ist der Zähler 2 in dem Fliessdiagramm. Das heisst,
jedesmal, wenn eine Geschwindigkeitsschätzung gemacht wird, wird der Zähler achrittweise
fortgeschaltet
und die Intervallgeschwindigkeit für die nächste
Schicht wird bestimmt.
-
Die Schleife, die bei der Instruktion 31 beginnt und mit der Instruktion
80 endet, wird eine Anzahl von Malen durchlaufen, die gleich der anzahl von LhYERS
abzüglich der Anzahl der fallengelassenen Schichten ist. LAYERS abzüglich der anzahl
von fallengelassenen Scnichten ist die Anzahl von Schichten in den Ausgangsinformationen.
-
@s wird nunmehr die Subroutine VELEST betrachtet.
-
Die Instruktion 5 iii der ubroutine bestimmt VMIN, die minimale annehmbare
Intervallgeschwindigkeit, die in diesem Fall 4500 Fuss je sekunde beträgt.
-
Die Instruktionen 10 bis 67 umfassen do-Schleifen, welche die dubroutine
ERRVEL abrufen. Hierdurch werden Kurven für geschätzte Laufzeiten für gegebene V(lAYERS)
und Zeitneigung, die in RAANG (LAYERS) ist, erzeugt. Dies ist eine Newton'sche Iteration,
um das Minimum in der Fehler kurve zu finden.
-
Die Instruktionen 55 bis 96 bestimmen, auf welcher beite des minimums
der letzte gute iiert der Newton'schen Iteration liegt.
-
Die Instruktionen 97 bis 163 bestimmen drei Punkte in der Näne des
Minimums, zwei auf einer Seite und einen auf der anderen beite, Dies wird ausgeführt
durch Inkrementieren der letzten weschwindigkeitsschätzung durch Schritte von 500
Fuss je Sekunde oder weniger.
-
Die Instruktionen 165 bis 215 sind eine parabolische Kurvenanpassung
der drei Punkte zum Schätzen von drei neuen Punkten, Die Kurvenanpassung wird wiederholt,
bis V(LAYERS) innerhalb + 5,0 Fuss je sekunde des Minimums liegt.
-
Die Instrustionen 225 bis 640 prüfen das Minimum, um zu sehen, dass
die Schwellbedingungen oder Schwellwertbedingungen erfüllt sind und V(LAYERS) grösser
als VMIN ist.
-
Die verbleibenden Subroutinen sind ein eil des Programms zur Modellbildung
und Laufzeiterzeugung. Die Subroutine
TWOTPHI fordert die Subroutinen,
welche die Schritte 4 und 5 gemäss dem Fliessdiagramm ausführen. Insbesondere erzeugt
die Subroutine TWOUPHI in den Instruktionen 3115 bis 3210 die Geschwindigkeiten,
Neigungen und Winkel des Schichtmodells.
-
Danach wird die Subroutine Bild abgerufen. Diese Subroutine benutzt
das Bildverfahren für Strahlenspurbildung.
-
Die Laufzeitkurven zu der interessierenden Zwischenfläche werden
in einer do-Schleife berechnet, die mit der Instruktion DO 1000 K-l, NTR beginnt.
Die do-Schleife endigt an der Instruktion 1000. Diese Schleife umfasst die -Subroutine
TILTTWT, welche die Strahlen und die Strahlenzwischenflächen für irgendeinen bestimmten
anfänglichen Strahlenwinkel erzeugt. Die keststellungen 620 bis 50 zwingen den Strahl,
zu dem gewünschten Quelle-Empfänger-Abstand oder-Strecke zu konvergieren, und zwar
unter Anwendung einer Gradiententechnik.
-
Die Nützlichkeit der erfindung wird am besten anhand der kig. 6 bis
9 demonstriert0 Fig. 6 ist ein gestapelter Aufeeichnungsaschnitt variabler iläche,
auf welchem die Oberflächenschichten bis herab zu 2,0 Sekunden alle nach links geneigt
sind. bei 1,0 bis 1,8 Sekunden erscheint eine geringfügige Verdickung. Die Reflektion
"vent " hat eine Zeitneigung Null.
-
Felddaten aus dem gleichen bereich wurden in einem VIP-Verfahren
verarbeitet, um die in den Fig. 7 und 8 dargestellten Daten zu erhalten. Die mittleren
Geschwindigkeiten und die Reflektionszeiten T0 von diesen Daten wurden als Eingang
für das Verfahren gemäss der Erfindung verwendet. Die Srgebnisse sind in Fig. 9
dargestellt. Gemäss Fig. 9 wurde die Intervallgeschwindigkeit bestimmt in Übereinstimmung
mit der vorliegenden Erfindung, und zwar sowohl mit Neigung, wie vorstehend beschrieben,
als auch mit der Annahme9 dass keine Neigung vorhanden war. Beide diese Intervallgeschwindigkeiten
sind in Fig. 9 als Funktion der Aufzeichnungszeit aufgezeichnet.
weiterhin
sind in Fig. 9 Pfeile dargestellt, die in dieser Figur die verlagerten Erscheinungen
orten bzw. lokalisieren und ihr Neigungsverhalten angeben. bs ist zu bemerken, dass
die Erscheinung A, wenn sie.in Raumkoordinaten verlagert ist, nach rechts geneigt
ist, so dass eine Neigungsumkehr erzeugt wird, die auf dem Zeitabschnitt der Fig.
6 nicht sichtbar ist. Die durch die oben genannten Figuren dargebotenen Daten sind
nachstehend in Tabellenform wiedergegeben.
-
Zusammenfassung mit Neigung TWO-WAY TIME INT-VEL DIP XDIST YDIST .140
4902 0.0 0 343 .310 5515 0.0 0 812 .620 5893 7.9 193 1714 .945 6773 12.9 492 2781
1.240 7711 6.8 446 3937 1.530 8004 .4 220 5115 1.790 9574 -16.2 -681 6295 2.125
8329 - 9.3 -329 7743.
-
Es sind viele Modifikationen des Vorstehenden ersichtlich. Beispielsweise
wurde vorstehend das Verfahren beschrieben in Verbindung mit Schichtmodellierung
und Strahlenspurenbildung nach dem Snell'schen Gesetz. latsächlich stellt Schichtmodellierung
nach dem Snell'schen Gesetz eine Annäherung dar, d.h. einen besonderen Fall der
Strahlenspurbildung in einem sich kontinuierlich ändernden Medium. Im allgemeinen
bestimmen Eikonal-Gleichungen den Strahlenpfad in einem sich kontinuierlich ändernden
Medium. Diese sind beispielsweise in "Mechanical Radiation", Lindsay, R.B., MoGraw
Hill, New York 1960, Seiten 41-47 beschrieben. Bei Anwendung auf den allgemeinen
Fall ist die Erfindung anwendbar auf reflektierende Zonen an Stelle von Schichten
und die Strahlenpfade werden aus xiikonal-Gleichungen an Stelle nach dem Snell?schen;
Gesetz bestimmt.
-
Min dreidimensionales Verfahren wird verwirklicht durch
eine
einfache Ausdehnung des oben offenbarten Verfahrens Zu den Werten von TO, Va und
, wird der Wert von #, das Azimut, hinzugefügt. Diese Parameter bestimmen vollständig
die elastische Schichtung für den dreidimensionalen 1 all, wie es die drei Parameter
TO, Va und # für den zweidimensionalen Fall tuno Das Verfahren ist in der Verwirklichung
dem oben beschriebenen Verfahren identisch, wobei eine bubroutine für ein dreidimensionales
Modell verwendet wird, um die reihe von Laufzeiten TNij zu erzeugen, wenn N sich
auf die N-te Schicht bezieht, i sich auf die i-te Strecke oder auf den i-ten Abstand
in X-Richtung, j sich auf den j-ten Abstand in der Y-Richtung bezieht, wenn X und
Y ein willkürliches orthogonales Koordinatensystem auf der Erdoberfläche darstellen
und der Orientierung bzw. Ausrichtung-der Flächenreihe entsprechen.