Sarma D.D. Geostatistics with Applications in Earth Sciences

Подождите немного. Документ загружается.

186 Geostatistics with Applications in Earth Sciences

************************************************************

C READ CO=NUGGET

EFFECT;

SILL,

RANGE , HORIZONTAL (HAF) AND

C VERT

ICA

L (VAF)

ANIS

OTRO

PIC

FACTOR

************************************************************

1 =1 , 2

READ( 6

,*)

CO(

, C(

,RANGE(

,ANG(

,HAF(I),VAF(

************************************************************

C PRI NTS THE PA

RAME

TERS READ

************************************************************

WRI TE (2 , 2 5 ) CO ( I )

( I )

,RANGE(I)

, ANG (I ) , HAF (I ) , VAF (I )

2 5 F

ORMAT

(lX

, ' VARIOG

RAM

ETER

S',

/ 5X ,

, T2

FlO

/ ,

1 5X,

'C'

,T22

,F13

, /

,5X

, 'RANGEE' , T2 5 , F10 .

/ , 5X, ' ANGLE OF ROTATION

2 T2

5,F1

0 . 2

,/5X,

'HOR . ANI SOTRO

PY'

, T 2 5

,F1

3 5X , ' VER

ANI SOTR

OPY',

T2 5

,F1

.2)

=ANG(I)

*0 .

0174532

9 2

ST(I)

=CO

S(RA)

S I T ( I ) =SI N (RA)

CONTINUE

************************************************************

WX2 =W

WY2=W

Y/ 2 . 0

WX4 =W

WY4=WIDY/ 4 . 0

WX8 =W

IDX/

8 . 0

WY8=WIDY/

************************************************************

1 = 0

4 0 1

=1+1

************************************************************

C READS DATA WITH X, Y COOR

DINATES.

HERE

=ACCUMULATION,

C Z2 =GRADE, Z3

=THICKNESS

. THE ORDER CAN BE

CHANGED

NEEDED .

************************************************************

READ( 6

0)XE(

,YN(I)

,Zl

, Z2

3 0

FORMAT(5F10

. 0)

EL =O

D(I

,1)

=Z2

GD(I

,2)

=Zl

D(I

,3)

=Z3

ZL(I)

=EL

DPT ( I

IF(XE(I)

Q. O. O) GO TO 5 0

GO TO 4 0

5 0 NN

-1

IF(KEY2

0)G

O TO

************************************************************

C WRI TES DATA WI TH

COO

RDI NATES .

REM

OVE C BEFORE WRI TE

STATEMENT

************************************************************

Computer Software 187

************************************************************

BLOCK

THE DEPOSIT

CALCULATE

THE VARI ANCE OF A

54 I

=l,

WRI TE (2 , 52 ) DPT (I ) , XE (I ) , YN (I ) , ZL ( I )

,GD(I,1)

,GD(I

,3)

,GD(I

,2)

FORMAT

( '

',FlO

5X,

3F10.

3F1 0 . 2)

CONTINUE

IDENTIFICATION OF

BLOCKS

I N GRID E

NVIRONME

NT.

FOUR

SITUATIONS

ARISE .

DECREASING

WNWARDS

AND

INCREASING

FROM LEFT TO RIGHT, 2) Y

I NCREASING

WNWARDS

AND

DECREASING FROM LEFT TO RIGHT , 3) Y

DECREASING DO

WNWARDS

AND

DECREASING FROM LEFT TO RIGHT 4) Y

INCREASING

WNWARDS

AND

X I S INCREA

FROM

LEFT TO RIGHT .

'CONTROL' TAKES CARE OF THESE VARI OUS SITUATIONS .

CALL

BLOCK

(

BLKVAR

(

, 1 )

CALL

BLOCK

(

BLKVAR

(2) , 2 )

WRI TE (2

BLKVAR

FORMAT

(/5X,

' BLOCK VARIANCES

:',

T25,

F12

T40

F12

51)

3 3

(11,22

, 33 ,44)CONTROL

=YMAX+

0 . 5*WIDY

=XMIN- 0 . 5* WIDX

GO TO 5 5

YMIN

=YMAX

XMAX

=XMIN

YM=YMIN- 0 . 5* WIDY

=XMAX+

0 . 5*WIDX

TO 5 5

XMAX

=XMIN

=YMAX+0

.5*

WIDY

=XMAX+0

.5*

WIDX

TO 55

YMIN

=YMAX

=YMIN-

5* WIDY

=XMIN- 0 . 5* WIDX

GO TO 5 5

667

************************************************************

NOTE: I F B(3) VALUE IS

NOT

VAILABLE

THEN

SET

ZERO

. IT IS

W 2D.

THE

VARIABLES

USED

HERE

ARE G

RAD

E AND

ACCUMULATI

(ACCU)

. THESE

CAN

CHAN

GED.

************************************************************

B(3)

=0 . 0

WRI TE (2 , 62 )

FORMAT(lX

0X, , KRIGED RESULTS ' )

WRI TE (2 , 6 67 )

ORMAT

(4X, 'X-C

ORD

' ,4X, ' Y-

RD' , 6X, ' G

DE' ,5X, ' S. E' ,6 X, ' ACCU. ' 1,

, 'S . E')

188

222

666

777

888

999

5555

Geostatistics with Applications in Earth Sciences

************************************************************

ROW IS

FIXED

COLUMN

VARIES

************************************************************

DO 80 J

=NBR,NER

TO(66

,77,88

,99)CONTROL

B(2)=Y

M-

J *WIDY

222

B(2)=YM+J *WIDY

222

B(2)=Y

M-

J *WIDY

GO TO 222

B(2)

=YM+J*

WIDY

************************************************************

AND

OES

FOR

VARI

OUS

COLUMNS

************************************************************

80 K=

NBC

NEC

(666,777

,888

999)CONTROL

B(l )=

+K*WIDX

5555

B(l )=

XM-

K*WIDX

5555

B(l)

=XM

-K*

WIDX

5555

B(l )=

+K*WIDX

************************************************************

FIND

NEAREST

OLES

FOR

KRIGING

EXCLUDES

THE

BLOCK

THE

ORDINATES

THE

OCK

ESTIMATED

MATCH

WI TH

THE

ORDINATES

THE

'NEAREST'

NEIGHBOURHOOD

OCK/H

OLE

************************************************************

=RMAX2

60 L=l,NN

IF(

B(2) . EQ.YN(L) .

D.B (l ) .EQ.XE(L)

)GO

=B(2) -YN(L)

=B(l) -XE(L)

=B(3) -ZL(L)

************************************************************

OMPUTES

THE

DISTANCE

BET

WEEN

OLES

************************************************************

DIS

=(DX*C

ST(2)+DY*SIT(2))**2+(HAF(2)*(DY*CST(2) -DX*SIT(2)))**2

1 +(VAF(2)*DZ)**2

IF(DIS.

.R2)GO TO 60

(11.

EQ.0)

100

IF (I1.GE .NK.

D.DIS.GT.DI ST

(II)

)

100

=II+1

DIST

(II)

=DIS

NDATA

(II)

(11.

. 1)

1000

99 98

99 99

Computer Software 189

************************************************************

NDATA

STANDS

FOR

THE

NEIGHBOURHOOD

DATA

HERE

THE

PROVISION

FOR

51 POINTS

WHICH

REASONABLY

GOOD

ARRANGE

NEIGH

BOURHOOD

DATA

INCREASING

ORDER

DISTANCE

CALLING

SUBROUTINE

ARRANGE

************************************************************

CALL

ARRANGE

(DIST,

NDATA

II)

IF(II.GT

.NK)II

=NK

CONTINUE

LE . 2 ) GO TO 80

************************************************************

PRINTS

THE

NEIGHBOURHOOD

ACTUALLY

TAKEN

FOR

KRIGING

DEPENDING

THE

SPECIFIED

LIMIT,

************************************************************

WRI TE (2 , 1000 )

FORMAT

(3X, ' NEI GHBOURHOOD

DATA

ACTUALLY

USED FOR KRIGING

99 99 I

WRITE(2,9998)XE(NDATA(I)) ,

(

ATA(I ) ) ,GD(

ATA(I ) , 1 )

,GD(NDATA(I)

,2)

FORMAT(4F1

0.2)

CONTINUE

************************************************************

KRIG THE

BLOCK

THERE ARE MORE

THAN

TWO

DATA

POINTS

BY INVOKING KRIG. FOR

************************************************************

76 M=1

CALL

KRIG(B,II

,M,ERROR)

************************************************************

ERROR

OCCURS

IN THE

MATRIX

, THEN:

************************************************************

(ERROR. EQ. 0 ) GO TO 61

WRI TE (2 , 73 )J , K

73 FORMAT('ERROR IN

MATRIX

SOLUTION BLOCK'

,2I5)

DO 90 0

=NDATA(KN)

WRI TE (2 , 85 )XE (LL) , YN (LL )

85 FORMAT(2F1

900

CONTINUE

TO 80

************************************************************

CALCULATE

THE KRIGED ESTIMATE

AND

THE KRIGING

VARIANCE

C ************************************************************

GR(M) =O

(KEYl.

GE . 1 )

707

FORMAT

(lX

, ' HOLE' .rii , 'NORTH' ,T22 , ' EAS

, T32 , 'ELEV'

,T42,

'GRADE',

T52,

WEIGHT

70 L

=l,

LL=NDATA(L)

190 Geostatistics with Applications in Earth Sciences

=AG+GA(L)*R(L)

************************************************************

KRIGES

FOR

THE

FIRST

AND

SECOND

VARIABLES

************************************************************

GR(M)

=GR(M)+GD(LL

,M)*R(L)

NTINUE

VARKG(M)

=BLK

VAR(M)

-AG

-R(II+1)

VARKG(M)

=SQRT(VARKG(M))

CONTINUE

************************************************************

NOTE

FOR

KRIGING,

THE

ELE

VATI

ORE

THE

NEAREST

HOLE

ASSIGNED

THE

BLOCK.

************************************************************

=NDATA(l)

ELE

=ZL(KK)

WRI TE(2, Sl )B(1 ) , B(2) , (GR(I),VARKG(I),I=1,2)

FORMAT(2F10.

0,4F10

.2)

CONTINUE

707

STOP

END

************************************************************

SUBROUTINE

BLOCK

FOR FOR

COMPUTING

BLOCK

VARIANCE

************************************************************

SUBROUTINE

BLOCK

(BLK

VAR

,M)

************************************************************

THIS

PROGRAM

COMPUTES

THE

BLOCK

VARIANCE

THE

DEPOSIT

FOR

GIVEN

SIZE B

USING

EQUALLY

SPACED

DATA

POINT

APPROXIMATION.

************************************************************

COMMON/BLK/

WX2

WY2

WX4,WY4,

WXS

,WYS

COMM

ON/GAM/CO

(2) ,C(2)

,RANGE(2)

,CST(2),SIT(2)

,HAF(2) ,VAF(2)

DIMENSION

P1(3)

,P2(3)

GMOY

=O.

P1(1) =

WX2

P1(2) =

WY2

P1(3) =0.0

P2(3)

=0.0

10 1=1,4

P2(1) =

WX2+(I

3)*

WX4+

WXS

10 J =1,4

P2(2)

=WY2+(J

-3)*WY4+WYS

CALL

COVAR

(P1, P2,

GH,M)

GMOY

=GMOY+GH

CONTINUE

BLKVA

CO(M)+C(M

GMOY/ 16.

RETURN

END

Computer Software 191

************************************************************

SUBROUTINE

KRIG

.FOR

FOR

COMPUTING

KRIGED

ESTIMATES

FOR

GRIDDED

BLOCKS

************************************************************

SUBROUTINE

KRIG(B,NS,M,ERROR)

COMMON

/SYS/IPT,IOUT,IPR,KEY1

COMMON

/KRG/NDATA

(51)

,GA

(51) ,R (51)

COMMON

/DRILL/DPT(3000) ,YN(3000) ,XE(3000) ,ZL(3000),GD(3000,3)

COMMON

/GAM/CO(2),C(2),RANGE(2),CST(2),SIT(2) ,HAF(2) ,VAF(2)

DIMENSION

P1(3)

,P2(3)

,B(3)

,A(51,52)

15 K

=l,NS

=NDATA(K)

P1(1)

=XE(KK)

P1(2)

=YN(KK)

P1(3)

=ZL(KK)

=K+1

IF(K1.GT.NS)GO

10 L=K1,NS

=NDATA(L)

P2(1)

=XE(LL)

P2(2)

=YN(LL)

P2(3) =ZL(LL)

CALL

VGAM(P1,P2,GH

,M)

A(K,L)

=CO(M)+C(M)

-GH

A(L,K)

=CO(M)+C(M)

-GH

CONTINUE

11 A(K,K)

=CO(M)+C(M)

CALL

COVAR(B,

,GMOY

,M)

NS1

=NS+1

NS2

=NS+2

A(K,NS2)

=CO(M)+C(M)-GMOY

GA(K)

=A(K,NS2)

A(K,NS1) =

A(NS1,

CONTINUE

A(NS1

,NS2)

=1.0

A(NS1

,NS1) =0.0

IF(KEY1.NE.1)GO

C WRITE(2, 90) (GA(IX),IX=1,NS2)

WRI TE(2, 18)

FORMAT(lHO,'AUGMENTED

MATRIX')

19 IX

=l,

NS1

C WRITE(2 , 90) (A(IX,IY)

,IY

=1,NS2)

CONTINUE

C90

FORMAT(15F7

.2)

CALL

WEI GHTS (A,R,NS1,ERROR)

************************************************************

C R

CONTAINS

THE

SOLUTION

192 Geostatistics with Applications in Earth Sciences

************************************************************

RETURN

END

************************************************************

C SUBROUTINE

ARRANGE

.FOR FOR

ARRANGING

THE

NEIGHBOURHOOD

DATA

C IN ASCENDING

ORDER

************************************************************

SUBROUTINE

ARRANGE

(DS ,NDATA,N)

C THIS

A SORTING ROUTING IN INCREASING

ORDER

DIMENSION DS

(1)

, NDATA

(1)

DNE

W=DS(N)

NEW=NDATA(N)

N1=N

-1

20 I

=1,N1

IF(DNE

W.LT

.DS(I)

)GO TO

CONTINUE

RETURN

30 JK =O

40 I =K,N1

=N1-JK

JK =JK+1

DS(J+1)

=DS(J)

NDATA(J+1) =NDATA(J)

CONTINUE

DS (K) =DNEW

NDATA(K)=NEW

C WRI TE (2 , 51 ) (DS(K) ,K

=l,N)

51 FORMAT(10F8

.2)

RETURN

END

************************************************************

C SUBROUTINE WEI GHTS . FOR FOR

COMPUTING

THE

LAMDAS

AND

C THE

LAGRANGIN

PARAMETER

'MD'

************************************************************

SUBROUTINE WEI GHTS (A,

X,N

ERROR)

C THIS SUBROUTINE

COMPUTES

LAMDAS

AND

'MD'

DIMENSION

A(51,52)

,X(51)

ERROR

MP=N+1

666

=l,N

C WRI TE (2 , 555 )A (I , I )

555

FORMAT(lX

,F10.2)

666

CONTINUE

10 I

=l,N

=I+1

10 J

(LEQ.

888

IF (A(I ,

.EQ. O)

TO 30

F=(

-A(J,I))/A(I

,I)

9 K=IP,MP

A(J,K) =A(J,K)+F*A(I,K)

CONTINUE

20 I =l ,N

X(I)

=A(I,N+1)/A(I,I)

WRITE

*, X(I )

CONTINUE

ALAM

=O.O

DO 40 I =1,N-1

ALAM

=ALAM+X(I)

ONTINUE

NNN

-1

WRITE(2,888)NNN

,ALAM

FORMAT

(lX,

'NEIGHBOURHOOD

=' ,14

RETURN

ERROR

END

Computer Software 193

SUM

LAMDAS

=' , FS. 2 )

************************************************************

OUTINE

COVARIANCE

OMPUTES

THE

COVARIANCE

EEN A

NEIGHBOURHOOD

DATA

POINT/

HOLE

AND

A BLOCK

USING

POINT

APPROXIMATION.

************************************************************

SUBROUTINE

COVAR(P1

,P2,GH,M)

THE

REFERENCE

POINT

A BLOCK IS ITS

CENTRE.

COMMON

/GAM/CO(2)

,C(2),RANGE(2),CST(2)

,SIT(2)

,HAF(2) ,VAF(2)

COMMON

BLK/WX2

,WY2

,WX4

WY4

,WX8,WY8

DIMENSION

D(3)

,P1(3),P2(3)

GMOY

D(3) =P1(3)

D(1) =P1(1)

-2*

WX4

WX8

10 1=1,4

D(1) =D(1)+

WX4

D(2) =P1(2)

-2*WY4-WY8

10 J =1,4

D(2)

=D(2)+WY4

CALL

VGAM(D,P2

,GH

,M)

GMOY

=GMOY+GH

10 C

ONTINUE

=GM

OY/16

RETURN

END

************************************************************

SUBR

OUTINE

VGAM

FOR

COMPUTING

THE

VARIO

GRAM

VALUES

FOR

VARI

OUS

DISTANCES

. M=l

MEANS

G-T

PRODUCT;M

=2 M

EANS

DE;

************************************************************

SUBROUTINE

VGAM

(P1, P2,

GH,

THIS

ROUTINE

OMPUTES

THE

VARIOGRAM

VALUE

FOR

GIVEN

DISTANCE

194

Geostatistics with Applications in Earth Sciences

C USING

THE

SPHERICAL

MODEL.

MMO

N/G

/CO(2 ) ,C(2 )

,RANGE(2

) ,CST

(2),SIT(2

)

,HAF

(2 ) ,

VAF

(2 )

NSIO

N P

l(3

) ,P2 (3)

=Pl (1)

-P2(

1 )

=Pl (2)

-P2(2

)

=Pl (3)

-P2(3

)

=SQRT

((XD

*CST

(

M)+YD

*SIT (M))*

*2+HAF(

*HAF

(M)*(YD*CST( M

*SIT(

)**2+

VAF(

M)*

VAF

(

M)*ZD*ZD

)

IF(D.GT

.RANGE(

M) )GO

=CO

(

M)+C

(M)*(1 .5*D

/RANGE

(M)

-0

.5* (D

/RANGE

))**3)

IF(D.EQ .O

.O)GH

RET

=CO

(

M)+C

(M)

RET

Bibliography

Bibliography is divid ed into two sections:

I. Deals with bibliographical list in the field

mathematical statistics and

related fields - specially those relevant to geological sciences (this

includes the references cited in the text also); and

II. Deals with a bibliographical list relevant to geological data analysis,

and geostatistical methodology as propounded by the French School.

SECTION I

Mathematical

Statistical Methods

Aitchison, J. and Bro wn, J.A.C., 1957.

The

Lognormal Distribution.

Cambridge University Press, London, 176 p.

Bartlett, M.S., 1946. On the Theoretical Specification

Sampling Properties

Autocorrelated Time Series. Jr!. Royal Stat. Soc., B8, 27.

Box, G.E.P. and Jenkins, G.M., 1976. Time Serie s Analysis, Forecasting

and Control, Holden-Day, San Franc isco, California, 575 p., 2nd Edn.

Cooley, W. and Lohne s, P.R., 1971. Multivariate Data Analy sis. Wiley and

Sons, New York, 564 p.

Cooley , J. and Tukey , S.W., 1965. An Algorithm for the Machine Calculation

Complex Fourier Series. Math. Comp., 19, 297-301.

Cramer, H., 1945. Mathematical Methods

Statistics (First publi shed in

Sweden, Uppasala by Almqvist & Wiksells). Princeton University Press.

Princeton, 579 p.

Cressie, N.A.C., 1996. Statistics for Spatial Data. John Wiley and Sons ,

New York, 900 p.

Doob, J.L., 1953. Stochastic Processes. John Wiley and Sons, New York,

654 p.

Feller, W., 1957. An Introduction to Probability Theory and its Applications,

Vol. I, 2nd Edn. Asia Publishing House, New Delhi, 461 p.

Feller, W., 1966. An Introduction to Probability Theory and its Applications,

Vol. II. Wiley Eastern Pvt. Ltd., New Delhi, 626 p.