stochastische modellering van de waterkwaliteit
van de Zuid-Nederlandse Noordzee
J.R. Moll
verslag onderzoek
R 2176 deel
R 2360
2. Relatie inputs-HISTOS chloride gegevens 4 2.1 Inleid ing , 4 2.2 Fysische processen en tijdschalen 5 2.3 ARMAX - modellering 7 2.4 Recursieve parameter schattingstechniek 10 2.5 Resultaten 13 2.6 Interpretaties en conclusies 13 3. Relatie meetnetten 16 3.1 Inleiding , 16 3.2 Lineaire regressie , 17 3.3 Resultaten 19 3.4 Interpretaties en conclusies. 23 4. Samenhang waterkwaliteltparameters 24 4.1 Inleiding 24 4.2 Hoofdkomponentenanalyse , 26 4.2.1 Inleiding , 26 4.2.2 Meetkundige voorstelling 27 4.2.3 Fartitie van de variantie 28 4.2.4 Dimensiereduktle 28 4.3 Resultaten , 29 4.4 Interpretaties en conclusies,. 33
5. Ruimtelijke correlatie , 35 5.1 Inleiding 35 5.2 Beschrijving van de Kriging-techniek ,..* 37 5.2.1 Mathematische formulering van het probleem. 37 5.2.2 Krigen bij afwezigheid van trend 38 5.2.3 Aanwezigheid van een ruimtelijke trend , 40 5.2.4 De gegeneraliseerde intrinsieke hypothese.*... 41 5.2.5 Schatting van de gegeneraliseerde covariantiefunktie 43 5.3 Resultaten 45
5.4.1 Keuze van modelorde en s t ruk tuur. , 47 5.4.2 Interpretatie van de modelcoëfficiënten..» • 47 5.4.3 Het effect van outliers onder de meetgegevens op de berekeningen. 48 5.4.4 Interpretatie van de figuren..., , «. 49 5.5 Representativiteit van Hlstos-lokaties... •, «, 50 5.6 Interpretatie en conclusies » 51 5.7 Het gebruik van een verdunningslijn * 52 5.8 Interpretatie en conclusies ,., , 53
6. Integratie van deelstudies * 54 6.1 Veranderingen in de tijd...,, » 54 6.2 Veranderingen in de plaats 54 6.3 Differentiatie per waterkwaliteitsparameter 55 6.4 Veranderingen in meetnet , , 55
7. Samenvatting? conclusies en aanbevelingen ., 56
LITERATUUR
2.1 Inputsignaal ui als funktie van de tijd 2.2 Inputsignaal u2 als funktie van de tijd 2.3 Parameter 91 als funktie van de tijd 2.4 Parameter 92 als funktie van de tijd 2.5 Parameter 63 als funktie van de tijd 2.6 Parameter 94 als funktie van de tijd
2.7 Chloride-concentratie, gemeten te HA-10, dec 1981 2.8 Chloride-concentratie, berekend te HA-10, dec 1981 2.9 Autocorrelatie - funktie van de residuen
4.1 Scatterplot van de scores van de variabelen op de geroteerde factoren nr 1 en 2
4.2 Scatterplot van de scores van de variabelen op de geroteerde factoren nr 1 en 3
5.1 Lokaties van het Helikopter-meetnet periode jan - mei 1981 5.2 Lokaties van het Helikopter-meetnet periode juni 1981 - dec 1982 5.3 Gelnterpoleerde waarden van de concentratie NH4, 17-02-'81
5.4 Standaardafwijking interpolaties van de concentratie NH4, 17-02-'81 5.5 Gelnterpoleerde waarden van de concentratie NO2/NO3, 17-02-'81
5.6 Standaardafwijking interpolaties van de concentratie NQ2/NO3, 17-02-'81 5.7 Geinterpoleerde waarden van de concentratie Ntot, 17-02-'81
5.8 Standaardafwijking interpolaties van de concentratie Ntot, 17-02-'81 5.9 Geinterpoleerde waarden van de concentratie P04, 17-02-'81
5.10 Standaardafwijking interpolaties van de concentratie P04, 17-02-'81 5.11 Geinterpoleerde waarden van de concentratie Ptot, 17-02-'81
5.12 Standaardafwijking interpolaties van de concentratie Ptot, 17-02-'81 5.13 Geinterpoleerde waarden van de concentratie Si, 17-02-'81
5.14 Standaardafwijking interpolaties van de concentratie Si, 17-02-'81 5.15 Geinterpoleerde waarden van de concentratie Cl, 17-02-'81
5.16 Standaardafwijking interpolaties van de concentratie Cl, 17-02-*81 5.17 Gelnterpoleerde waarden van de concentratie Cpart, 17-02-'81
5.20 Standaardafwijking interpolaties van de concentratie DOC, 17-02-*81 5.21 Geinterpoleerde waarden van de concentratie NH4,08-12-'81
5.22 Standaardafwijking Interpolaties van de concentratie NH4, 08-12-'81 5.23 Geinterpoleerde waarden van de concentratie NO2/NO3, 08-12-'81
5.24 Standaardafwijking interpolaties van concentratie NO2/NO3, 08-12-'81 5.25 Geinterpoleerde waarden van de concentratie Ntot, 08-12-'81
5.26 Standaardafwijking interpolaties van de concentratie Ntot, 08—12—*81 5.27 Geinterpoleerde waarden van de concentratie P04, 08-12-'81
5.28 Standaardafwijking interpolaties van de concentratie P04, 08-12-'81 5.29 Geinterpoleerde waarden van de concentratie Ptot, 08-12-'81
5.30 Standaardafwijking interpolaties van de concentratie Ptot, 08-12-'81 5.31 Geinterpoleerde waarden van de concentratie Si, 08-12-81
5.32 Standaardafwijking Interpolaties van de concentratie Si, 08-12-'81 5.33 Geinterpoleerde waarden van de concentratie Cl, 08-12—*81
5.34 Standaardafwijking interpolaties van de concentratie Cl 08-12-'81 5.35 Gelnterpoleerde waarden van de concentratie Cpart, 08-12-'81
5.36 Standaardafwijking interpolaties van de concentratie Cpart, 08-12-'81 5.37 Geinterpoleerde waarden van de concentratie DOC, 08-12-'81
5.38 Standaardafwijking interpolaties van de concentratie DOC, 08-12-'81 5.39 Geinterpoleerde waarden van de concentratie NH4, 09-02-'82
5.40 Standaardafwijking Interpolaties van de concentratie NH4, 09-02-'82 5.41 Geinterpoleerde waarden van de concentratie NO2/NO3, 09-02-'82
5.42 Standaardafwijking Interpolaties van de concentratie NO2/NO3, 09-02-'82 5.43 Geinterpoleerde waarden van de concentratie Ntot, 09-02-'82
5.44 Standaardafwijking interpolaties van de concentratie Ntot, 09-02-'82 5.45 Geinterpoleerde waarden van de concentratie P04, 09-02-'82
5.46 Standaardafwijking Interpolaties van de concentratie P04, 09-02-'82 5.47 Geinterpoleerde waarden van de concentratie Ptot, 09-02-'82
5.48 Standaardafwijking interpolaties van de concentraties Ptot, 09-02-'82 5.49 Geinterpoleerde waarden van de concentratie Si 09-02-'82
5.50 Standaardafwijking interpolaties van de concentratie Si, 09-02-'82 5.51 Geinterpoleerde waarden van de concentratie Cl, 090-02-'82
5.52 Standaardafwijking interpolaties van de concentratie Cl, 09-02-'82 5.53 Geinterpoleerde waarden van de concentratie Cpart, 09-02-'82
5.55 Geinterpoleerde waarden van de concentratie DOC, 09-02-'82
5.56 Standaardafwijking interpolaties van de concentratie DOC, O9-O2-'82 5.57 Geinterpoleerde waarden van de concentratie NH4, 20-07-'82
5.58 Standaardafwijking interpolaties van de concentratie NH4,20-07-'82 5.59 Geinterpoleerde waarden van de concentratie N02/NO3, 20-07-'82
5.60 Standaardafwijking interpolaties van de concentratie NO2/NO3, 20-07-'82 5.61 Geintepoleerde waarden van de concentratie Ntot, 20-07-'82
5.62 Standaardafwijking interpolaties van de concentratie Ntot, 20-07-'82 5.63 Gelnterpoleerde waarden van de concentratie P04, 20-07-'82
5.64 Standaardafwijking Interpolaties van de concentratie P04, 20-07-'82 5.65 Geinterpoleerde waarden van de concentratie Ptot, 20-07-f82
5.66 Standaardafwijking interpolaties van de concentratie Ptot, 20-07-'82 5.67 Geinterpoleerde waarden van de concentratie Si, 20-07-'82
5.68 Standaardafwijking Interpolaties van de concentratie Si, 20-07-'82 5.69 Gelnterpoleerde waarden van de concentratie Cl, 20-07-'82
5.70 Standaardafwijking interpolaties van de concentratie Cl, 20-07-'82 5.71 Geinterpoleerde waarden van de concentratie Cpart, 20-07-'82
5.72 Standaardafwijking Interpolaties van de concentratie Cpart, 2O-O7-'82 5.73 Geinterpoleerde waarden van de concentratie DOC, 20-07-'82
5.74 Standaardafwijking interpolaties van de concentratie DOC, 20-07-'82 5.75 Geinterpoleerde waarden van de concentratie Cl, 17-02-'81, op basis
Histos lokaties
5.76 Standaardafwijking interpolaties van de concentratie Cl, 17-02-'81, op basis Histos lokaties
5.77 Geinterpoleerde waarden van de concentratie Cl, 08-12-'82, op basis Histos lokaties
5.78 Standaardafwijking interpolaties van de concentratie Cl, 08-12-'82, op basis Histos lokaties
5.79 Geinterpoleerde waarden van de concentratie Cl, 09-02-f82, op basis
Histos lokaties
5.80 Standaardafwijking interpolaties van de concentratie Cl, O9-O2-'82, op basis Histos lokaties
5.81 Gelnterpoleerde waarden van de concentratie Cl, 20-07-'82, op basis Histos lokaties
basis Histos lokatles
5.83 Geinterpoleerde waarden van de residu-concentraties Ntot, 08—12—'81 5.84 Standaardafwijking Interpolaties van de residu-concentraties Ntot,
08-12-'81
5.85 Geinterpoleerde waarden van de residu-concentraties P04, 08-12-'81 5.86 Standaardafwijking interpolaties van de residu-concentraties P04,
08-12-f81
5.87 Geinterpoleerde waarden van de residu-concentraties Si, 08-12-'81 5.88 Standaardafwijking interpolaties van de realdu-concentraties Si,
08-12-'81
5.89 Geinterpoleerde waarden van de residu-concentraties Ntot, 09-02-'82 5.90 Standaardafwijking interpolaties van de residu-concentraties Ntot,
09-02-'82
5.91 Geinterpoleerde waarden van de residu-concentraties P04, 09-02-'82 5.92 Standaardafwijking interpolaties van de residu-concentraties P04,
09-02-'82
5.93 Geinterpoleerde waarden van de residu-concentraties Si, 09-02-'82 5.94 Standaardafwijking interpolaties van de residu-concentraties Si,
b X-coëfficlënt nr j bij input-variabele nr i (i - 1,..3; j - 1
t..q )
c^ MA-coëfflclënt nr at (m - l,...r)
d. getransformeerde AR-coëfficlënt nr 1 (i » l,..«p')
g(.) getijstand
p aantal AR-coëfflciënten
p. modelparameter (i » 1,.».4)
p' maximum van p en r
p'' minimum van p en r
q
±aantal X-coëfflciënten bij intput variabele nr i (i • 1,...3)
r aantal MA-coëfficiënten
t tijdindex
u.(.) input variabele nr 1 (i
aI,...3)
v(.) stroomsnelheid langs de kust in ZW-rlchtlng
y(.) gemeten chloride concentratie
C1„
Q(.) berekende chloride concentratie ten NO van Histos-lokatie
Cl
z w(.) berekende chloride concentratie ten ZW van Histos-lokatie
M
9parameter in schattingsalgorlcme
P(.) schatting van covariantiematrix van 9
Q (.) Haringvliet-afvoer
H
Q
T(.) Rijn-afvoer te Loblth
jjQ (.) Schelde-afvoer
T
. getransponeerde
cx(.) hoek tussen windrichting en landinwaarts gerichte normaal op de kust
S^ looptijd van inputvariabele nr i
e parameter in schattingsalgoritme
e(.) witte ruis proces
X(.) vergeetfactor
9(.) parametervector
9 schatting van de parametervector
T lntegratievariabele
<|) (.) toestandsvector
b regressie-coëfficiënten x verklarende variabele y_ te verklaren variabele
y(.) meting van de te verklaren variabele y(.) schatting van de te verklaren variabele y gemiddelde van y
var variantie J kostenfunktie e(.) residu
U verwachting
E covariantieraatrix u vector met dimensie p
x, z stochastische vector met dimensie p \ Lagrange-multiplicator
T
getransponeerde h hoofdkomponent X eigenwaarde nr i tr spoor v/e matrix U ortogonale pxp matrix
I identiteitsraatrlx van dimensie p
fi factor nr i
B matrix van factor ladingen
B' gereduceerde matrix van factor ladingen (factormatrix) y willekeurig meetpunt
c coëfficiënt in gegeneraliseerde covariantiefunktie
cd karakteristieke correlatie-afstand
e
1,e„ exponenten
f (x) polynoora in x met index m
h afstandvektor
i,J Indices voor meetpunten
k orde
l(k) aantal polynotnen van orde kleiner of gelijk k
ra index voor polynoora
n aantal meetpunten
u Lagrange multiplicator
x punt in twee-dimensionale ruimte
x„ willekeurig, niet bemeten punt
x meetpunt
z(x) geregionaliseerde stochastische funktie
z schatter voor z
z schatting voor z
E verwachtingsoperator
G gegeneraliseerde aangroeiing
K(.) gegeneraliseerde covariantiefunktie
N(h) aantal puntenparen in afstandklasse h
P, (x) polynoom van orde k in x
8 coëfficiënt van een polynoom
y semivariogram
Y experimenteel semivariogram
S Dirac-delta funktie
A weegfaktor
W verwachting
£. ,§„ coördinaten in twee-dimensionale ruimte
p(.) correlatie-coëfficiënt
1. Inleiding en vraagstelling
In 1981 en 1982 werd door de hoofdafdeling Milieu en Inrichting van de Delta-dienst een Inventariserend onderzoek verricht naar de kwaliteit van het kust™ water voor het Deltagebied. Dit onderzoek was er vooral opgericht relaties te vinden tussen enerzijds de waterkwaliteit voor de kust en anderzijds de afvoer van Rijn, Maas en Schelde alsmede de waterbeweging langs de kust.
Daarnaast werden de seizoeninvloeden onderzocht*
De vraagstelling van dit onderzoek is mede ingegeven door de ingrijpende ver-anderingen die zich in dit gebied zullen voltrekken. Hierbij kan op dit moment onder andere gedacht worden aan het tot stand komen van de Oosterscheldedam en de plannen voor werken langs de Noordzeekust voor o.m. energieopwekking, aan-landing van goederen en berging van baggerspecie.
In het kader van bovengenoemde onderzoek hebben met behulp van een helicopter vierwekelljkse bemonsteringen plaatsgevonden van 66 3 68 meetpunten, rond hoogwater kentering* Een groot aantal waterkwaliteitsparameters waaronder nutriënten zijn daarbij onderzocht.
Parallel hieraan hebben waarnemingen plaatsgevonden met HISTOS meetpalen, waar hoogfrequente metingen van de waterstand en soms ook de stroomsnelheid, wind-snelheid en -richting, golfhoogte, sallniteit en temperatuur gedurende meer-dere jaren hebben plaatsgevonden.
De uit het onderzoek verkregen gegevens kunnen op verschillende wijzen worden verwerkt ora de relaties tussen de waterkwaliteit en de externe faktoren te on-derbouwen. In dit verband kunnen worden genoemd:
- het onderzoeken van direkte (statistische) verbanden tussen verkregen meet-reeksen, leidend tot (stochastische) black-box modellen.
- het ontwikkelen van verklarende modellen, waarin de invloed van de water-beweging en processen expliciet tot uitdrukking komen.
Beide benaderingen verschillen in aanpak en in karakter van het resultaat. In een aantal opzichten kunnen zij elkaar aanvullen.
Bij brief d.d. 26 oktober 1984 (knr 22292) is het Waterloopkundig Laboratorium door de Deltadienst gemachtigd beide benaderingen parallel uit te werken. Bij brief d.d. 9 september 1985 (knr 21236) is tevens machtiging verleend een aan-tal nadere vragen uit te werken, betrekking hebbende op de stochastische mo-dellering. In dit rapport wordt verslag uitgebracht van de stochastische ana-lyse, alsmede van de uitwerking der nadere vragen.
Bij deze stochastische analyse moeten er, gezien de aard en de hoeveelheid van de beschikbare gegevens, vier aspekten worden onderscheiden.
1. De temporele korrelaties tussen wind, getij- en afvoercijfers (de input-grootheden) en de waterkwaliteitsparameters zoals gemeten met de HISTOS meetpalen.
2. De relatie tussen de waterkwallteitsgegevens van het HISTOS meetnet, en de overeenkomstige gegevens van het helikoptermeetnet op, of nabij, de lokatie van de HISTOS-palen.
3. De samenhang van diverse waterkwaliteitsparameters van het hellkopterraeet-net op de onder 2. beschouwde lokaties. Hierbij kunnen extra lnputgroothe-den in de samenhang een rol spelen.
4. De ruimtelijke korrelatie en ruimtelijke patronen voor deze waterkwali-teitsparameters uit het helikoptermeetnet.
De struktuur van het stochastisch model is weer te geven als in figuur 1.1, waarbij de gebruikte symbolen de volgende betekenis hebben:
v (x ,t) je waterkwaliteltsparameter uit het HISTOS-tneetnet op lokatie
x en tijd t
je waterkwaliteitsparameter uit het helikopter-meetnet op lokatie
x en tijd t
aantal lokaties in HISTOS meetnet aantal lokaties in helikopter meetnet
MHIS aantal waterkwaliteitsparameters In HISTOS meetnet
aantal kwaliteitsparameters In helikopter meetnet
Figuur l.l geeft het kader aan voor dit verkennend onderzoek. Indien er sta-tistische relaties zijn te identificeren voor elk blok in figuur 1.1 zijn er in feite rechtstreekse verbanden te leggen tussen externe inputs en de water-kwaüteitsparameters op elk punt in het helikoptermeetnet. Deze verbanden kun-nen diekun-nen als basis voor voorspellingen in vergelijkbare situaties.
De genoemde vier aspekten bepalen de opzet van dit verslag.
In hoofstuk 2 t/m 5 worden achtereenvolgens de afzonderlijke aspekten nader aan de orde gesteld. Hoofstuk 6 bevat de integratie van de resultaten van deze vier deelstudies en hoofdstuk 7 tenslotte de conclusies en aanbevelingen.
Het hier beschreven onderzoek is uitgevoerd en gerapporteerd door ir. J.R. Moll i.s.ra. ir. T. Schilperoort.
2. Relatie inputa-HISTOS chloride gegevens
2.1 Inleiding
Voor het leggen van de relatie hydro-meteorologische inputs - Histos-chloride gegevens is gebruik gemaakt van een adaptieve ARMAX-modellering; dit is een tijdreeksmodellering, waarbij de parameters kontinu worden aagepast aan de ak-tuele meetgegevens.
De belangrijkste elementen van deze ARMAX-modellering zijn:
- een voorbewerking van de meetgegevens uitgaande van bekende fysische rela-ties
- het modelleren van deze bewerkte gegevens met een simpele ARMAX-structuur - een adaptieve parameterschattingstechniek
Deze drie elementen komen resp. in paragraaf 2.2, 2.3 en 2.4 aan de orde. Re-sultaten, verkregen met deze modellering zijn opgenomen in paragraaf 2.5. In-terpretaties en conslusies volgen in paragraaf 2.6. Doel van de modellering is om variaties in het op de Hlstos-lokaties gemeten chloridegehalte te verklaren uit variaties In de meteorologische omstandigheden en de grootte van de ri-vierafvoeren.
De voor dit doel beschikbare meetgegevens zijn opgenomen in Tabel 2.1.
soort meetgegeven lokatie bemonsteringsinterval
chloride, -2.50 NAP HISTOS 10 min chloride, -8.00 NAP HISTOS 10 min getijwaterstand HISTOS 10 min windsnelheid B6-2, OS-12 10 min windrichting BG-2, OS-12 10 min waterstand ' Loblth 24 uur afvoer Schelde 10 dagen
Tabel 2.1 Beschikbare meetgegevens
De maand december 1981 is geselekteerd als geschikte periode om het model op te zetten. In die periode is er sprake van hoogwatergolven op de grote rivie-ren, waardoor de identificatie van de Invloed van de grootte van de
rivieraf-In Tabel 2.1 staan twee series Histos-chlorlde-metlngen vermeld. Alleen het chloride op -2.50 m NAP is in dit onderzoek betrokken.
ARMAX-tnodellen zijn de laatste jaren frekwent toegepast in zowel waterkwall-teits-, als waterkwantiteitsstudies, zie bijv. 0'Connell en Clarke (1981), Beek en van Straten (1983), Young (1984).
2.2 Fysische processen en tijdschalen
Voor het toepassen van een ARMAX-modeliering moeten zinvolle inputsignalen worden geselekteerd. Deze inputsignalen zijn al dan niet bewerkte meetgegevens van fysische grootheden. Voor het selekteren van de inputsignalen moeten de relevante fysische processen worden beschouwd, en de karakteristieke
tijdschalen hiervan.
Er is onderscheid gemaakt tussen drie soorten processen, gelet op de richting ervan, zie Tabel 2.2.
aoort Proces
A Dwars op de kust
1. Lozing Haringvlietsluizen bij laagwater
2. Afvoer Schelde
3. Stroming t.g.v. lokale wind 4. Opzet t.g.v. wind op
Noordzee-schaal
B Langs de kust
1. twaalf uurs getijstroom 2. maandelijks getij
3. stroming, geïnduceerd door lokale wind 4. opzet t.g.v. wind op Noordzee-schaal kar. tijdschaal 6 uur enkele dagen enkele dagen enkele dagen enkele uren week enkele uren enkele dagen In de vertikaal 1. stratificatie/destratlficatie 2. neerslag/verdamping enkele dagen enkele uren
Tabel 2.2 Fysische processen en karakteristieke tijdschalen
Niet alle in Tabel 2.2 opgenomen processen zijn expliciet meegenomen bij de constructie van de inputsignalen voor het ARMAX-raodel. Zo ontbraken gegevens om het stratificatie/destratificatie proces (Cl) redelijk te kunnen beschrij-ven. De invloed van neerslag en verdamping (C2) op het chloridegehalte is ver-waarloosd gelet op de ordegrootte ervan.
Bij de stroming t.g.v. lokale wind (A3) wordt gedoeld op stroming waardoor re-latief zout water uit zee wordt aangevoerd. Er is hier dan sprake van stroming geïnduceerd door aanlandige wind. Gelet op de geringe strijklengte van aflan-dige wind wordt stroming zeewaarts t.g.v. lokale aflanaflan-dige wind niet meegeno-men in het model.
De opzet t.g.v. wind op Noordzee-schaal (A4, B4) is niet rechtstreeks, maar indirect vla de getijstand gemodelleerd. Dit zelfde geldt voor stroming langs
de kust t.g.v. lokale wind (B3).
Het maandelijks getij is gelet op de relatief grote karakteristieke tijdschaal
ervan niet expliciet opgenomen in het model. Via de tijdsvarierende parameters
van het ARMAX-model is deze maandelijkse cyclus echter wel in te brengen.
Bovengenoemde vereenvoudigingen brengen het aantal te construeren
inputsigna-len op drie:
1) het effekt op het chloridegehalte t.g.v. stroming langs de kust uit het
Zuidwesten (incl. evt. invloed Schelde)
2) het effect op het chloride-gehalte t.g.v. stroming langs de kust uit het
Noordoosten (incl. evt. invloed Haringvliet)
3) het effekt op het chloride-gehalte t.g.v. stroming richting kust t.g.v.
lo-kale aanlandige wind.
2-3 ARMAX-modellering
Een ARMAX-model is een tijdreeksmodel met de volgende algemene structuur:
p q q
p q
7(t) -
1£
1a
iy(t-i) +jïj bj uj (t-6
rj)
met:
y(t) = chloride-concentratie op tijdstip t op Histos-lokatie
s(t) = witte-ruis proces
u.(t) « Input-variabele nr i (1 « 1,...3)
a = AR-coëfficlenten nr i
b* =• X-coëfficienten nr i bij input variabele 1 (i » 1....3)
c » MA-coëfficienten nr m
m
p = aantal AR-coëfficienten
r = aantal MA-coëfficienten
<H = aantal X-coëfficienten bij input variabele i (i = 1,...3)
Bij gegeven tijdreeksen y(t), u ^ t ) (i=l,...3) bestaat het toepassen van een ARMAX-model uit twee stappen:
i) Het indentificeren van de modelstructuur of modelorde d.w.z. de waarden van de parameters p, q., q„, q«, r, 6., 6„ en 5»
il) Het schatten van de waarden van de parametervector 0 ,
1 1 2 2 3 3 T
met 9 - [a ,,,,a , bn,...b , bn,...b , bn,...b , c1,...c ]
i p u q^ v ^2 u 4g l r
Het schattingsalgoritme voor de parameter vector 6 (stap li) wordt in de vol-gende paragraaf behandeld. Voor het vaststellen van de modelstructuur (stap i) zijn er in principe twee mogelijkheden: ofwel er wordt op grond van a priori aanwezige fysische kennis een structuur vastgesteld, ofwel er wordt op grond van statistische kriterla, zoals het Akaike Information Crlterion (Akaike,
1974), een keuze gemaakt.
Een combinatie van beide mogelijkheden is tenslotte ook nog mogelijk: op grond van de identificeerbaarheid van de parametervector 9 (stap ii) kan de a priori gekozen modelstructuur worden aangepast.
In dit onderzoek is op deze wijze te werk gegaan.
De a priori gekozen modelstructuur is
y(t) - ax y(t-l) + a2y(t-2) + bju^t) + bju^t) + b ^ U ) + e(t) + c^Ct-1)
(2.2)
met:
y(t) • te modelleren chloride-gehalte
u^(t) • inputsignaal, representerend de invloed op y(t) t.g.v. stroming langs de kust uit NO
U2(t) a inputsignaal, representerend de invloed op y(t) t.g.v. stroming langs de kust uit ZW
Ug(t) a inputsignaal, representerend de invloed op y(t) t.g.v. stroming richting de kust uit NW
e(t) a ruisproces
P = 2, q
x» O, q
2= O, q
3= O, r = 1, 6f Sf &f 0.
1 2 3 T
Hiermee wordt de parametervector 6 = [a.,a
2»b
o,b
o,b
o,c,] .
Bij deze structuur behoort een tijdstap van 3 uur.
De gemaakte keuze is als volgt toe te lichten. Op grond van de in Tabel 2.2
vermelde karakteristieke tijdschalen lijkt de keuze van 3 uur voor de tijdstap
goed mogelijk. Gelet op de alternerende aanwezigheid van u
]en u„ en de
tijd-schaal van die aanwezigheid (6 uur, oftewel een halve getijde-cyclus) ligt het
opnemen van meer b-parameters niet voor de hand. Voor het beschouwde probleem
en de hierin een rol spelende fysische processen is het introduceren van
loop-tijden niet noodzakelijk; bij de hierna te behandelen constructie van de
inputsignalen wordt dit nader uitgewerkt. Het aantal AR- en MA-parameters (cq.
2 en 1) lijkt in eerste instantie voldoende.
Definitie van inputsignalen
•
Uj(t) =» [ci
N0(t-l) - y(t-l)] . f ( ƒ V ( T ) dx) (2.3)
met:
C-J-NCI " chloride concentratie ten NO van de Histos-paal
y « chloride concentratie op de Histos-paal
t =• tijdindex
v =» stroomsnelheid langs de kust in ZW-richting
f(x) = x als x > 0
= 0 als x < 0
De stroomsnelheid v is niet als meting beschikbaar. Gekozen is voor gebruik
van de wel beschikbare getijstanden g:
• Uj(t) = [ci
N0(t-l) - y(t-l)] . f(g(t) - g ( t - O ) (2.4)
• u
2(t) - [Cl
zw(t-1) - y(t-l)] • f(g(t-l) - g(t)) (2.5)
met C l
z w= chloride-concentratie ten ZW van de Histos-paal
met:
Ü)(T) =• windsnelheid
a(-r) » hoek tussen de windrichting en landinwaarts gerichte normaal op de kust
De in (2,4) en (2.5) voorkomende variabelen C 1N O en C lz w kunnen op twee
manie-ren een waarde krijgen;
- ofwel de chloride-concentratie in het NO resp ZW is als meting beschikbaar aangezien daar een Histos-paal staat
- ofwel deze concentratie in het NO resp ZW wordt berekend uit: t
C 1N 0( t ) " ^ N O ^ "1* " Pl * J QH ( T ) d T + P2 * U3( t ) < 2'7 )
t
Clzw(t) - Clzw(t-1) - p3. ƒ Qg (T)dT + p4. u3(t) (2.8)
met:
QH » afvoer uit Haringvlietsluizen
Qe *• afvoer van de Schelde
o
p, a modelparameter i • 1,...4
De afvoer uit de Haringsvlietsluizen is als volgt geschat:
QH(t) - 2 . Qjt-24) . S(t) (2.9)
met:
Q L » afvoer te Lobith
S(t) - 0 bij vloed (sluis dicht) » 1 bij eb (sluis open)
De modelparameters P l i n (2.7) en (2.8) worden gecalibreerd aan de hand van
karakteristieke situaties.
2.4 Recursieve parameterschattingstechniek
In de literatuur worden vele ldentifikatie-algoritmen aangeboden.
Ljung en SöderstrSm (1983) geeft een overzicht en bevat een algemeen raamwerk voor de verschillende methoden. Ook het in deze studie toegepast
gemodiflceer-de adaptieve kleinste kwadraten algoritme wordt door genoemgemodiflceer-de auteurs
behan-deld. De geïnteresseerde lezer die meer details over het hieronder behandelde
algoritme wil weten wordt dan ook naar dit boek verwezen.
Het algoritme is gebaseerd op het kleinste kwadraten krlterium: de te
minima-liseren doelfunktie bestaat uit de kwadraten van gemeten waarden en door het
model berekende waarden. Het algoritme is recursief, hetgeen betekent dat
nieuw binnenkomende meetgegevens direkt kunnen worden gebruikt voor het
even-tueel bijstellen van de schattingen voor de modelparameters zonder dat het
ge-hele meetverleden opnieuw moet worden doorgerekend. De adaptiviteit van het
algoritme biedt de mogelijkheid in-de-tijd-varlerende parameters te volgen.
Onder verwijzing naar de notaties uit vergelijking (2.L) stel
p
1• =3 minimum van p en r
p' = maximum van p en r (2.10)
en definieer
^i
= ci ~
al 1 < p''
-
ci p " < i < r
- ~
ai p
( 1< 1 < p
1 =• l,...p' (2.11)
Defineer de parametervector 9 door:
0
T- [ d ^ . - . d , bj,...bj ,bj-,...b ,b
t...b
0
T- [d^.-.dp, bj,...bj ,bj-,...b2 ,b3
t...b3
>c
1 >...c
r] (2.12)
Defineer voorts:
*
T(t) - [y(t)..y(t-p' + 1 ) , u
1(c-fi
1)...u
1(t-6
1-q),
u
3(t-$
3)..u
3(t-6
3-q
3), -$
T(t-l) ê(t), -*
T(t-2) ê(t-l),.,.-$
T(t-r) 9(
(2.13)
waarbij de in bovenstaande vergelijking voorkomende schatting Q(
t) gegeven
wordt door het volgende algoritme:
9(t) - 8(t-l) + a(t-l) P(t-l) *(t-l) [y(t) - $
T(t-l) 8(t-l)] (2.14)
a ( t - l ) - [1 + $
T( t - l ) P ( t - l ) K t - n T
1(2.16)
P(0) - YIJ 8(0) - 0; j - P ' + q + r + s + m (2.17)
T
(2.18)
De parameters s en M_ moeten worden opgegeven; zij bepalen de karakteristieken van de 'vergeetfactor' X(t).
De matrix P is te interpreteren als een schatting voor de covariantiematrix van de schattingsfouten in 9.
In bovenstaand algoritme komt een numeriek gezien ongunstige stap voor: het berekenen van de matrix P vla (2.15). Dit kan leiden tot singulariteit van P. Dit is opgevangen door het algoritme te herschrijven. De resulterende
verge-lijkingen zijn gebaseerd op een UDUT -factorisatie van de matrix P. Voor
2.5 Resultaten
Be berekeningsresultaten met de a priori gekozen raodelstructuur (verg. (2.2)), wezen uit dat modificaties van deze structuur gewenst waren. Er bleek dat de parameter b« slecht Identificeerbaar was. Dit betekent dat het effect van ver-zilting t.g.v. stroming opgewekt door lokale aanlandlge wind niet expliciet identificeerbaar'was. Aangezien dit effect langs indirecte weg, nl. vla de in-putsignalen u.(t) en u^(t), ook al in het model gebracht was, en aangezien in de beschouwde periode (december 1981) niet sprake was van een overheersende aanlandige wind, is dit resultaat te verklaren. Hierop Is gekozen voor het schrappen van de terra b„u3(t) uit de modelstructuur (2.2)
Een tweede aanpassing bestond uit het schrappen van de term c..e(t-l) in (2.2). Het resulterend model is daarmee robuuster:
y(t) = aiy(t-l) + a2y(t-2) +" bju^t) + bju2(t) + e(t) (2.19)
Een laatste aanpassing bestond uit het inbouwen van een onder- en bovengrens voor berekende (niet gemeten) chloride concentratie ten NO van de beschouwde Histos-paal HA-10. Deze grenzen zijn gesteld op 18 gr/l en 5 gr/l.
Het gebruik van verg. (2.7) zou zonder deze grenzen kunnen lelden tot onrea-listisch hoge of lage concentraties.
Resultaten met het model (2.19) zijn weergegeven In de figuren 2.1 t/m 2.9. Figuur 2.1 en 2.2 bevatten de inputsignalen, figuur 2.3 t/m 2.6 de parametere-voluties, figuur 2.7 en 2.8 resp. gemeten en berekende chloride concentraties
(voorspelling 3 uur vooruit) en figuur 2.9 de autocorrelatlefunktle van de re-siduen. (
2.6 Interpretaties en conclusies
Het adaptieve parameterschattingsalgorltrae kent een inregelperiode van onge-veer 50-70 rekenstappen. Dit is goed zichtbaar in figuur 2.1 t/m 2.6 die resp. de Inputsignalen en de parameterevoluties bevatten. De parametervector
Flguur 2.1 bevat uj(t) gedefinieerd in (2.3). De sterk negatieve 'golf' in uj(t) rond tijdstip 104-120 (13-15 dec *81) correspondeert met een afvoergolf op de Rijn waarvan de top Loblth passeert op 12 december 1981. De alternerende aanwezigheid van U|(t) en U2(t) hangt samen met de richting van de getijde-stroom voor de kust.
Het verloop van de parameters in de figuren 2.3 t/m 2.6 geeft aanleiding tot twee opmerkingen:
1) na het inschakelverschijnsel variëren de parameters langzaam rond een zeker nivo
ii) de inputsignalen hebben de eigenschap de parameters continu aan te sturen, m.a.w., de inputsignalen zijn 'rijk' genoeg
Indien de parameters sterk fluctueren wijst dit op een instabiel model. Indien de parameters volledig konvergeren naar vaste nivo's betekent dit dat het ge-kozen model de realiteit perfekt beschljft, of dat het inputsignaal geen in-formatie bevat (bijv. nulslgnaal als input).
Het parameterverloop in figuur 2.3 t/m 2.6 is op het eerste gezicht bevredi-gend te noemen.
De identificeerbaarheid van de parameters is vastgesteld aan de hand van de berekende covariantlematrix P uit (2.15)
Het verloop van gemeten en berekende chloride concentratie bij HA-10 is uitge-zet in figuur 2.7 en 2.8. Wat opvalt is dat de berekeningsresultaten minder 'slingeren* dan de gemeten waarden. Dit is een gevolg van het feit dat de waarden van de modelparameters per tijdstap slechts in beperkte mate kunnen worden aangepast. Dit wordt gecontroleerd door de ingestelde
vergeet-factor M t ) in (2.18). Door deze vergeetvergeet-factor anders in te stellen is het mo-gelijk de meetwaarden uit figuur 2.7 aanzienlijk nauwkeuriger te reproduceren door modelberekeningen dan die in figuur 2.8 weergegeven. Het parameterverloop wordt dan echter veel grilliger.
Figuur 2.9 bevat het autocorrelogram van de residuen, de verschillen tussen gemeten en berekende waarden. He blijkt dat er geen autocorrelatie in de resi-du-reeks aanwezig is. Was dit wel het geval dan impliceerde dit dat het adap-tief schattingsalgoritme onvoldoende informatie onttrok aan de laatst binnen-gekomen meetwaarden.
Conclusie Is dat het model In staat is een kwantitatieve relatie te leggen tussen, hydrometeorologische Inputs en het chloridegehalte voor de Zuid-Neder-landse Noordzee kust. In concreto is deze relatie berekend voor het chloride-gehalte zoals gemeten op de Histos-meetpaal HA-10. De periode december 1981 is doorgerekend. Voor andere perloden en andere lokaties zijn soortgelijke bere-keningen uit te voeren. Gelet op de geringe variabiliteit in de op dit moment ter beschikking staande (Tabel 2.1) cijfers van de Schelde-afvoeren wordt ver-wacht dat het identificeren van de Invloed van deze input op het chloride ge-halte niet eenvoudig zal zijn.
Tot slot enige opmerkingen over voorspellingen.
De voorspelhorizon van het model (2.19) is beperkt tot 3 uur, aangezien ge-bruik gemaakt wordt van meetwaarden per 3 uur bij de constructie ven het In-putsignaal U2« Voor het voorspellen over een langere termijn is een
voorspel-ling van Inputsignaal u (t) nodig. Hiervoor kan op analoge wijze weer gebruik worden gemaakt van een ARMAX-model.
Door een aantal ABMAX-modellen op deze wijze te koppelen Is een geïntegreerd model op te zetten, waarbij de Inputsignalen uitsluitend externe invloeden
(wind, afvoeren) representeren. Indien voorspellingen van deze externe invloe-den te leveren zijn voor een zekere voorspelhorizon, zijn met dit geïntegreer-de mogeïntegreer-del ook chlorigeïntegreer-de gehaltes te voorspellen voor geïntegreer-deze termijn.
3. Relatie meetnetten
3.1 Inleiding
In dit hoofdstuk wordt verslag gedaan van het onderzoek naar het koppelen van meetgegevens van het Histos-meetnet aan die van het Helikopter-meetnet. Met het Histos-meetnet wordt hoog frekwent (elke 10 min.) op een gering aantal lo-katies een aantal waterkwaliteitsparameters bemonsterd; met het Helikopter-meetnet zijn gemiddeld eens per maand waterkwaliteitsparaineters bemonsterd op een groot aantal lokaties. Een belangrijke vraag is op welke wijze de twee meetnetten complementair zijn.
Hiertoe is onderzocht of de verschillende meetmethoden belangrijke verschillen opleverden voor op dezelfde plaats en op hetzelfde tijdstip bemonsterde water-kwaliteitsparameters.
Op de Histos-lokaties zijn chloride en temperatuur bemonsterd. Aangezien de temperatuurgegevens van de Helikopter-metingen niet betrouwbaar geacht werden is alleen getracht op basis van de chloride-metingen een vergelijking tussen de meetgegevens van de twee meetnetten te maken. Elke Hlstos-meetpaal waar chloride gemeten wordt heeft twee sensoren, één op -2.50 m en ëën op ± -8.00 m. Bij het relateren van de chloride gegevens uit het Helikopter-meetnet aan die van het Histos-meetnet is gebruik gemaakt van de gegevens van de sensor op -2.50 ra.
Een overzicht van de gebruikte Hlstos-, en Helikopter lokaties is opgenomen in Tabel 3.1 en 3.2. Het geringe aantal (één) stations in Tabel 3.1 is geen ge-volg van een doelbewust gemaakte keuze maar vloeit voort uit de beschikbaar-heid van meetgegevens.
Hlstos-paal R.D.M.-coördinaten Helikopter lokatie R.D.M.-coördinaten
6.OS-4 37870; 408740 65 40420; 409337
Histos-paal R.D.M.-coördinaten Helikopter lokatie R.D.M, -coördinaten
1.
2.
3. 4. 5. 6. 7. HA-10 BG-2 BG-8 OS-10 OS-15 OS-4 0S-9N 49827; 32815; 46220; 21388; 22202; 37870; 35338: 431569 421394 419185 403545 401070 408740 405235 1329
22 43 42 36 39 49122; 30151; 45835; 20151; 22862; 37124; 33711; 433281 420882 419717 403516 400839 409514 406492Tabel 3.2 Lokaties Histos en Helikopter periode juni '81 - nov '82
De Helikopter gegevens zijn vergaard tijdens vluchten rond tijkentering (hoog-water bij Hoek van Holland). Uitzonderingen hierop zijn de vluchten op
14-4-'81, l4-l-'82, 9-ll-f82 en ll-ll-'82.
Voor het vergelijken van de chloride gegevens van de beide meetnetten is ge-bruik gemaakt van lineaire regressierekening. Deze methode wordt in paragraaf 3.2 kort toegelicht. Voor een uitgebreidere beschrijving wordt verwezen naar handboeken zoals, Kreiszig (1970).
In paragraaf 3.3 zijn de resultaten opgenomen en in paragraaf 3.4 volgen de interpretatie hiervan en de conclusies.
3.2 Lineaire regressie
Het enkelvoudig lineair regressiemodel heeft de volgende vorm:
y •> a + Hierin bx is + s *•*
y
Xe
a,b =« te verklaren variabele = verklarende variabele =» residu = parameters (3.1)Voor de 'te verklaren' variabele y is de Helikopter chloride-meting gekozen, voor de 'verklarende' variabele x is de Histos-chloride-metlng op -2.50 m ge-kozen.
Bestaat de gegevensverzameling uit n metingen, dan is meting y(j) (Uj<n) te beschrijven met:
y(j) = a + b.x(j) + E(j) (3.2)
Beschikken we over de parameters a en b en de verklarende variabele x(j), dan is hiermee een voorspelling te maken van y(j):
y(j) - a + b.x(j) (3.3)
De parameters a en b zijn uit de metingen {(y(j), x(j)), Uj<n} te bepalen mef een kleinste kwadraten schattingsmethode:
min J(a,b) (3.4) a, b
met J - 2 £2(j) = t [y(j)- a - b.x(j)]2 (3.5)
J-l 3-1
Partieel differentieren naar de parameters a en b levert op:
f£°°
en
f " 0 (3.6)
een stelsel van twee lineaire vergelijkingen met twee onbekenden.
Van interesse is welk deel van de variantie van de variabele y verklaard wordt door het model. Een schatting voor de variantie van y is:
* In
var (y) - - Z [y(j) - y ]
2(3.7)
n j»l 1 n
waarin y - - Z y( j) (3.8)
n j-l '
Formule (3.7) is te herschrijven als:
n j-l n j«l
i n
.
i
n
. 2 n . <
3'
9>
- s [y(j) - y(j)]
2+ - E [y(j) - y ]
2+ -
zDe eerste term in het rechterlid is te interpreteren als een schatting voor de
residuele variantie, het niet door het regressiemodel te verklaren deel van de
variantie van v_ . De tweede term in het rechterlid is te interpreteren als een
schatting van het door het regressiemodel verklaarde deel van de variantie van
I •
De derde term in het rechterlid is nul bij gebruik van de kleinste-kwadraten
schattingsmethode voor de parameters a en b.
Tenslotte is van belang te weten welke betrouwbaarheid men toe mag kennen aan
de met (3.6) gevonden parameterwaarden.
Uitwerken van (3.6) levert:
b - ? y(j) [x(j) - x] / ï [x(j) - x]
2(3.10)
a * y - b x (3.11)
Uit (3.10) en (3.11) zijn rechtstreeks uitdrukkingen voor var (b) en var (a)
af te leiden:
var(b) - - S [y(j) - a - b x(j)]
2/ Z (x(j) - x )
2(3.12)
n j-1 j-1
var (a) * {- Z [y(j) - a - b x(j)]
2} . {- + x
2/ ï [x(j) - x]
2} (3.13)
n j-1 n j-1
De grootheden a . (var(a))^ en b. (var(b))"* hebben een t-verdeling met n-2
vrijheidsgraden als het regressiemodel correct is. Hiermee zijn
betrouwbaar-heidsintervallen voor de parameters te berekenen. Ook valt hiermee te toetsen
of de parameters significant van nul verschillen.
3.3 Resultaten
Er bleken 111 getallenparen (Helikopter-chloride, Histos-chloride) beschikbaar
te zijn. Zie voor een overzicht van deze gegevens Tabel 3.4. De Helikopter-,
en Hlstos-gegevens zijn niet in dezelfde eenheden weergegeven: mg/l resp, 0,1
g/l.
Visuele inspektie van Tabel 3.4 levert op dat de correspondentie tussen
chlo-ridemetingen te HA-10 (stationscode 1) en Helikopter lokatie 13 soms vrij
slecht is. Dit is te verklaren uit het feit dat de chloride-gradiënt in het gebied rond HA-10 dicht bij de Haringvlietsluizen vrij groot is. Een betrekke-lijk gering verschil in meetlokatie of beraonsteringstijdstip kan al tot aan-zienlijke verschillen in gemeten chloride-concentratie leiden. Op grond hier-van is besloten het station HA-10 niet in de berekeningen mee te nemen.
Gemiddelden en standaardafwijkingen van de twee reeksen chloride gegevens zijn opgenomen in Tabel 3.3.
gemiddelde standaard afwijking eenheid
Helikopter 17623.45 877.54 mg/l Histos 177.09 8.15 0,1 g/l
Tabel 3.3 Statistiek van Helikopter en Histoa chloride metingen.
De correlatiecoëfficiënt tussen de twee variabelen bedroeg ,999, zodat m.b.v. de Hlstos-metingen 99% van de variantie van de Helikopter-chloride cijfers verklaard kan worden.
De gevonden regressievergelijking is:
Helikopter (j) - 99,5 Histos (j) + e(j) mg/l ( U j < U l ) (3.U)
met als standaardafwijking voor de coëfficiënt voor Histos .15.
Bij het schatten van de parameters van deze regressievergelijking is a proiri gesteld dat de constante (parameter a uit (3.1)) een waarde nul krijgt, op fysische gronden.
3 1 O'S 1 7 6 1 G 4 1 4 8 1 V 6 1 « ' 8 i n f, i-b 8 1 H 61 £• 8 1 C 7 1 fe 8 1 f. i 7 2 1 8111714 810?14 81Cfcl3 81 Of. 11 810811 8 1 0 Ê 1 1 810908 8199£8 813908 810996 811013 811013 8111' 13 811110 811 U u 81111? 61113 0 811110 8 1 1 2 i) fl 8 1 1 2 s J 8 81120K 811203 811208 8201 17 8 2 0 1 0 7 8 2 $ l ij 7 8213107 820114 8 2 0114 820114 82011* 8 2 0114 820209 8202-37 820209 8 Z 0 2 D 9 82 r2?9 8 2 0 3 39 8 2 ;• 3 3 9 820309 8 2 :j 3 Ü 9 S 2 1 3 3 a 82 ?4 ^ 8 ? 1' 4 -j 5 S 2 C 4 3 5 32J4 23 8 2 3 4 '.13 3 2 f) r. ) <=> 8 2 0 5 U 4 1 " v ? 7 3 3 . 4 : IE.21 * e « 3 ? 1 2 , 5 "5 1 * . 4 8 1 * . 2 4 1 4 « 3 £ 1 1 . 1 9 13.04 12.39 1 2 . 5 D '5i,5 6 -1 -1 . 4 9 •-11.23 11,34 13.30 1 4 • 4.4 14.56 11.21 12.25 14.15 13.45 13.55 19.4S 1 0 . 4 8 • -12.40 1 2 , 1 7 1 2 . 2 7 1 1 . l r. 13.23 •12.57 13.09 10.36 11.43 1 4 . 1£ 13.5* 14.05 13.03 14,08 16.07 15.71 15.52 12.19 13,29 15.14 14.48 14.39 11.4^ 12.55 1 2 . 2 "'• 1 4 . 5 1 1 4 , 3 2 1 1 . 1 2 1 2 . 1 0 1 1.42 f. h ? 6 7 2 4 C 1 2 4 6 7 2 • 4 6 7-£ f 7 1 2 4 1 •• £ 4 6 7 2 4 f 7 1 2 4 6 7 1 2 4 f. 7 1 2 4 f, 7 1 X 1 4 7 1 1 17'=.'"!?. 1 f. ? 6 Q 1 8 1 ? 8 1 7 c' 7 7 1 7 7 £ "*. iaa53 1 7Q6 C I7«e4 17^41 17)60 179 83 16 " ë 7 17^46 1771 ? 1743? • 17135 17191 19345 1 ? r) 9 ? 18887 i6<;9? 18472 1 fi ? 9 ? 17393 17463 1U52 16&0P 16368 15534 1 4 7 6 8 1 6q 7? 17199 15812 16553 14514 15239 157 35 15043 15375 1511E 17731 17414 17 4 0 4 17*72 16525 17711 17659 17475 1 7 4 4 4 1 8 ü 19 17 6 ? S 1 f; q ? ° 17616 16845 1 * f. 9 o 173 6 C 1 f. *•(• 2 171 176 1P? 3 8 3 17? Ifi3
is r
181 178 17?. 186 172 176 18 2 182 176 175 191 191 189 158 184 183 181 T75 1*5 172 163 153 15C 166 174 157 163 127 157 1.64 1 5 f. 161 142 17Ê ! 76 175 175 151 1 7a ? 7 5 3 7 5 1 7T If ? 1 ? ': ! 7? ?76 169 1 ft f 1,76 1 6Q8 2 (: 3 C 4 B i l *"j y" "^ t 8 2 ö & ] 1 3 2 0 S J l 3 2 3 5 01
sais-n
3 2 0 6)1
S20S22
8 2 0622
3 2 062 2
82 1622
8 2 3 5 2 2
8 2 3 7 2 0
&a
8 3-3 7 2 fr
82)317
3-2-.m-7'
8-a-iat?
820817
820817
8 2 d 3 1 7
820 914
8 2^914
82Ü914
821312
821-112
821 '712
821
J12-821 J 1.2
821-312
8211.19
321139
32111
1*
321109-a21.1U9
821111
821111
321111
321111
8 2 1 U 1
321126
821125
821126
321126
321124-831126
0 9,49
K . 5 2
13,22
12.4^
12.24
12.34
14,52-14,22
16.34
16.1 1
16.20
12.55
15,52
13,24
15.14
11.41
12.44
12,14
.14-.44
14,17
14.28
12.03
11.27-13,41
fia.2i 9.2711,12
10.46
10.57
.18.33
1 1 . 4 510.59
11.54
11.27
13.44
13.17-12.17
11,41
11,50
2 1797^
* 17 5&1
5 17?37
f 17PS7
7 177C7
1 U737
3
sr6
7
1
2
•*.5
&
7
1
6
1
174P1
17727
17767
17748
17129
1312 5
17754
18.175
17851
17772
18024
13 213
17914
173 74
18114
17391
18593
18162
lÖS'Jl
17147
13939
5 135
6 i ° r
7
1
2
3
5
6
7
l
18.) 47
185 8.T
18666
1845-5
18155
1 8 7 ) 3
18.194
13914
18724
3 13 9 ? 9
185 1413599
13? 38
1« 121
X9L65
18 J 4 2
1 Sc.
1 ?nter,
177
177
.173
18?
177
183
177
133
178
179
188
173
J.S.1
1 9 1
13 8
183
186
189
177
185
185
182
18-7
13 3
135
132
184
186
.1.82
134
179
185
18 C
17 7
J. 7 7
179
Tabel 3.4 De voor de regressie-analyse beschikbare gegevens. De code
identificeert de Histos-paal, zie Tabel 3.2
3.4 Interpretaties en conclusies
Uit Tabel 3.3 blijkt dat gemiddelden en standaardafwijkingen van de twee soor-ten chloride-metingen vrijwel gelijk zijn.
De correlatiecoëfficiënt van .99 geeft aan dat de metingen goed vergelijkbaar zijn. De regressievergelijking (3.14) wijkt weliswaar iets af van de lijn die bij zeer goede overeenstemming te verwachten is:
Helikopter (j) - 0 + 100 Histos (j) + e(j) mg/l (Kj<;lll) <3.15)
maar het verschil is te verwaarlozen gelet op de waarde van de standaardafwij-king voor de coëfficiënt.
Conclusie is dat de Helikopter-metingen en Histos-metingen op -2.50 m NAP, die op dezelfde plaats en tijd zijn uitgevoerd zeer goed vergelijkbaar zijn. De twee meetnetten kunnen t.a.v de bemonstering van de chloride-concentraties als complementair worden beschouwd.
4. Samenhang waterkwallteitparameters
4.1 Inleiding
Tijdens de Helikopter-vluchten zijn 26 waterkwaliteitsparameters op 66 (tot 1-6-81:68) lokaties bemonsterd. Voor het leggen van relaties tussen waterkwali-teitsparameters en hydro-meteorologische Inputa (hoofdstuk 5) is het wenselijk het aantal te beschouwen parameters te reduceren door gebruik te maken van hun onderlinge samenhang. In dit hoofdstuk wordt verslag gedaan van het onderzoek naar deze onderlinge samenhang.
Het onderzoek heeft zich beperkt tot 9 van deze 26 waterkwaliteitspararaeters, nl. - chloride - silicium - aramoniura - nitraat + nitriet - particulair stikstof - orto-fosfaat - particulair fosfor
- particulair, organisch koolstof (Cpart) - opgelost organisch koolstof (DOC)
Uit het Helikopter-bestand zijn die lokaties geselekteerd die corresponderen met de posities van de Histos-palen. Door deze selektie wordt voorkomen dat
t.g.v. aanwezige ruimtelijke correlatie van de waterkwaliteitsparameters veel doublures in de te analyseren gegevens voorkomen. Door de gemaakte keuze is de koppeling naar de hydro-meteorologische inputs (hoofdstuk 2) gewaarborgd. Een overzicht van de gebruikte lokaties wordt gegeven in Tabel 4.1 en Tabel 4.2,
Histos-meetpaal Helikopter-lokatie HA-10 BG-2 BG-8 OS-10 OS-4 OS-9 46 30 38 25 65 66
Tabel 4.1 Keuze van. lokatiea Helikopter-meetnet gerelateerd aan positie Histos-meetpalen periode jan - mei'81
Histos-meetpaal Helikopter-lokatie HA-10 BG-2 BG-8 OS-10/OS-15 OS-4 OS-9/OS-9N 19 29 22 43 36 39
Tabel 4.2 Keuze van lokatles Helikopter-meetnet gerelateerd aan positie His-tos-meetpalen periode juni '81 - nov '82
Hiermee zijn 6 (lokaties) x 27 (helikoptervluchten)» 162 metingen geselekteerd van de genoemde 9 waterkwaliteitsparameters.
Voor het onderzoeken van de samenhang tussen deze 9 waterkwaliteitsparameters is gebruik gemaakt van hoofdkomponentenanalyse. Deze multivariate statistische techniek wordt kort beschreven in paragraaf 4.2. Resultaten van het onderzoek staan vermeld in paragraaf 4.3. Interpretaties en conclusies in paragraaf 4.4 sluiten dit hoofdstuk af.
4.2 Hoofdkomponentenanalyse
4.2.1 Inleiding
Hoofdkomponentenanalyse Is een 'klassieke' techniek voor het structureren van samenhang In groepen van variabelen* Voorbeelden van toepassing op waterkwali-teitsstudies zijn bijv. Mahloch (1974), Steele & Matalas (1974), en Groot en Malwald (1976). Het ook op dit terrein regelmatig toegepaste factoranalyse is een verwante techniek, voor hetzelfde doel inzetbaar, maar gefundeerd op een geheel ander uitgangspunt. Zie voor een gedetailleerde uiteenzetting over bei-de technieken Barman (1976).
De onderlinge lineaire samenhang van p waterkwaliteitsparameters (In het ver-volg aan te duiden als de 'variabelen') wordt gekarakteriseerd door de corre-latleraatrix ervan. Uit deze p x p matrix is op het oog moeilijk meer te halen dan de samenhang van afzonderlijke paren van variabelen. Het is wenselijk een eenvoudige structuur te vinden die de samenhang van alle variabelen
be-schrijft. Deze structuur kan dan gebruikt worden om de (vele) variabelen samen te vatten tot slechts enkele nieuwe variabelen zonder al te veel informatie te verliezen. Eén van de mogelijkheden om deze structuur te vinden is gebruik te maken van hoofdkomponentenanalyse.
Hoofdkomponentenanalyse construeert uitgaande van de p waargenomen variabelen een p-tal nieuwe variabelen, de hoofdkomponenten, op zo'n wijze datï
i) de eerste hoofdkomponent verklaart zoveel mogelijk van de varlantle In de p variabelen; de tweede hoofdkomponent draagt daarna zoveel mogelijk bij aan de verklaring, etc. Uiteindelijk is iedere variabele exact een li-neaire combinatie van de hoofdkoraponenten.
il) de eerste hoofdkoraponent vat de variabelen zo goed mogelijk samen; de tweede hoofdkomponent geeft een (numeriek gezien) wat minder belangrijk aspect weer; etc.
iii) de hoofdkomponenten zijn ongecorreleerd.
Een hoofdkomponenten analyse is geslaagd te noemen als
i) de belangrijkste hoofdkomponenten goed interpreteerbaar zijn
4.2.2 Meetkundige voorstelling
De verzameling meetpunten van de p variabelen is voor te stellen als een pun-tenwolk in een p-dimensionale ruimte. Indien men de punpun-tenwolk projekteert op een lijn, dan is de spreiding van de projekties In de regel afhankelijk van de richting van de lijn: de richting van een lange as van de puntenwolk levert een grote variantie, de richting van een korte as een kleine variantie. Het zoeken naar samenhang tussen de p variabelen kan op deze wijze meetkundig wor-den vertaald in het zoeken naar projektle-assen in een p-dimensionale ruimte.
Stel x is p-dimensionaal verdeeld met verwachting u en covariantiematrix £. Stel u e Rp met ju| 3 1. De projektie van x op u is (u x) u . De variantie
T
-hiervan is u 2 u . Elke u correspondeert met een willekeurige asrlchtlng. We
m
zoeken nu die u, waarvoor g (u) = u Z u maximaal, minimaal, dan wel stationair is.
Met behulp van de Lagrange-multiplikator X is die u te berekenen uit:
(g (u) - XuTu)« 0 (4.1)
uTu - 1 (4.2)
Hieruit volgt: E u - Xu (4.3)
Dit betekent: de variantie in de richting van u is extreem d.e.s.d. als u een (genormeerde) eigenvector van £ is; de variantie is dan gelijk aan de bijbe-horende eigenwaarde van Z.
De projecties op de assen zijn de hoofdkomponenten van x. De le hoofdkomponent
is de projektie op de langste as, de 2e die op de êên na langste as, etc. Bij
een zeer langgerekte puntenwolk heeft men het gevoel dat de le hoofdkomponent
reeds een aardige lineaire samenvatting geeft van x....x . Bij geringe excen-triciteiten zijn méér hoofdkomponenten nodig om een adequate samenvatting van X...X te geven.
4.2.3 Partitie van de variantie
Veronderstel dat de covariantiematrix £ nlet-singulier is. Dan zijn alle ei-genwaarden X positief (Ki<p) en te ordenen in dalende volgorde:
X„>XA> ... >X >0 (4.4)
1 2 p
© T
De i -hoofdkomponent van x is de stochastische varlable h • u x , met u de ge-normeerde eigenvector van E bij de i -eigenwaarde.
De variantie van de hoofdkotnponent is gelijk aan de corresponderende eigen-waarde, var h - var (u x) • u E u » u Xu =• X . Hoofdkomponenten bij verschil-lende eigenwaarden zijn altijd ongecorreleerd*
De totale variantie van het stel x,...x is var x,+ var x„ ...+ var x -1 -p -1 -2 -p. (ook te schrijven als tr (E)). Eenvoudig af te lelden is dat deze totale va-riantie Invariant is onder ortogonale transformaties. (Als u ortogonaal,
T
dan II II • I„ ) • Met X^ (Ki<p) de eigenwaarden voor E, U ortogonaal pxp, en [z, ...s
var Zj < X.
var z, + var z„ < X. + X_ (4.5) etc.
Door over te gaan op hoofdkomponenten wordt dus een nieuwe partitie van de to-tale variantie bewerkstelligd, waarbij met elke nieuwe hoofdkomponent een zo groot mogelijk deel van de variantie geëxtraheerd wordt.
4.2.4 Dimensiereduktie
We noteren gestandaardiseerde hoofdkomponenten h als factoren f.
±
j
^
J i - l,...p <4.6)
Zoals we hoofdkomponenten uit kunnen drukken in de oorspronkelijke variabelen, zo kunnen variabelen worden uitgedrukt in hoofdkomponenten, of in factoren:
De getallen In matrix B kunnen worden opgevat als regressiecoeffIcienten van de variabelen x op de factoren f ( H K p , l<j<p)
-i "3 Stel nu k < p.
Door regressie te plegen op f. t/m f, I.p.v. op f. t/m f wordt een zo goed mogelijke benadering van x gegeven door z met dimensie k:
[ z
r. z
p]
T= B* [ f
r. f
k]
T(4.8)
De matrix Bf ontstaat uit B door weglaten van de laatste p - k kolommen van B.
Als globale kwaliteltsmaat voor deze benadering dient de fractie verklaarde varlantie van de totale variantie:
k p Z \ f l X 1=1 i I=-l i
Deze fractie raoent men wel de communaliteit van x met de andere variabelen vla de k factoren. Het restant vormt de specificiteit van x,. Met behulp van
(4.8) kunnen van een willekeurig meetpunt y -(y ,...y ) de factorscores 1 P
worden berekend, noem deze v =• (v , ...v. ):
v = (B1) y . (4.9)
4.3 Resultaten
Voor het uitvoeren van de hoofdkomponentenanalyse Is gebruik gemaakt van sub-routine FACTOR/TYPE » PA1 uit het SPSS-pakket. De analyse ia uitgevoerd, In de correlatieschaal. Zie voor documentatie van SPSS: Nie, Huil, e.a. (1975). Ta-bel 4.1. bevat gemiddelde en standaardafwijking van de beschouwde waterkwali-teltparameters. Het aantal meetpunten is in de berekening tot 138 gereduceerd door het ontbreken van enkele data in het oorspronkelijke bestand.
standaard afw. .0509 .4413 .1314 .0286 .0412 .3968 1730.3707 1.0260 .3121
Tabel 4.1. Gemiddelde en atandaardafwijkingen waterkwallteitsparameters
Basisinformatie over de aanwezige samenhang tussen de parameters wordt gegeven in de correlatiematrix, Tabel 4.2. variabele
NH4
NO2/3 Npart P04 Ppart SICL
CpartDOC
gemiddelde .0597 .4324 .1743 .0535 .0493 .3160 17140.3188 1.6634 1.5232NH4
NO2/3 Npart P04 PpartNH4
N02/3 Npart P04 Ppart SI CL Cpart DOC 1.00000 .65594 -.11632 .76521 -.03804 .74666 -.69049-.11141
.20045
.65594
1.00000 -.00683 .74525 .26863 .92305 -.81931 .18997 .04283 -.11632 -.00683 1.00000 -.11541 .56279 -.07904 -.18267 .53276 -.04336 ,76521 .74525 -.11541 1.00000 .09106 .81143 -.62991 .00810 .20099 -.03804 .26863 .56279 .09106 1.00000 .15912 -.20072 .75883 .03919 SI CL CpartDOC
NH4
N02/3 NpartP04
Ppart SI CL Cpart DOC .74666 .92305 -.07904 .81143 .15912 1.00000 -.72753 .11769 .05727 -.69049 -.81931 -.18267 -.62991 -.20072 -.72753 1.00000 -.04904 -.28395-.11141
.18997 .53276 .00810 .75883 .11769 -.04904 1.00000 -.24328 .20045 .04283 -.04336 .20099 .03919 .05727 -.28395 -.24328 1.00000Het overzicht van de eigenwaarden in Tabel 4.3. laat zien dat met 3 hoofdkom-ponenten 83,0% van de variantie verklaard wordt.
Factor eigenwaarde 1 4.09234 2 2.30751 3 1.06971 4 .57871 5 .37203 6 .22948 7 .18777 8 .12725 9 .03521
Tabel 4.3. Eigenwaarden van de correlatiematrix
Vervolgens is de factormatrix (Bf ult(4.8)) berekend, Tabel 4.4,
PCT v.d. VAR. 45.5 25.6 11.9 6.4
4.1
2.5 2.1 1.4 .4 GUM PCT 45.5 71.1 83.0 89.4 93.6 96.1 98.2 99.6 100.0 FACTOR 1 .84021 .93066 .02298 .87564 .23838 .93579 -.87121 .12420 .21452 FACTOR 2 -.25983 .07394 .78659 -.15363 -.86292 -.04101 -.04027 .89205 -.22027 FACTOR 3 .00271 -.17319 .26024 -.04455 .15846 -.20492 -.19577 -.16771 .91457 NH4 NO2/3 Npart P04 Ppart SI CL Cpart DOCTabel 4.4. Ongeroteerde factormatrix
Schattingen voor de communallteit van de waterkwaliteifcsparameters zijn opge-nomen in Tabel 4.5.
VARIABELE
NH4
NO2/3 Npart F04 Ppart SI CL Cpart DOC COMMUNALITEIT .77347 .90160 .68698 .79234 .82656 .91937 .79896 .83930 .93098Tabel 4.5 Geschatte communaliteit van de waterkwaliteitsparameters
De factorruimte wordt opgespannen door drie factoren. Boor deze factoren te roteren, hetgeen neerkomt op het kiezen van een ander ortogonaal assenstelsel is de interpreteerbaarheid van de factoren te bevorderen. Het resultaat van een z.g. VARIMAX-rotatie leverde de factormatrix in Tabel 4.6 op.
NH4 NO2/3 Npart P04 Ppart SI CL Cpart DOC
Tabel 4.6 Geroteerde factormatrix FACTOR 1 .85434 .93323 -.09012 .88487 .12804 .95394 -.82960 .05527 .11374 FACTOR 2 -.15290 .15704 .81495 -.05210 .89950 .03997 -.17563 .86235 -.04438 FACTOR 3 .14209 -.07754 .12129 .08149 .03274 -.08814 -.28262 -.30430 .95712
Voor de volledigheid volgt de transformatiematrix in Tabel 4.7. FACTOR 1 FACTOR 2 FACTOR 3
FACTOR L FACTOR 2 FACTOR 3 .98631 -.10002 -.13108 .L2052 .97986 .15921 ,11251 -.17283 .97850 Tabel 4.7 Tranaformatiematrix
De scores van de variabelen op de geroteerde factoren zijn als scatterplots opgenomen in Fig. 4.1 en 4.2.
4.4 Interpretaties en conclusies
Uit Tabel 4.3 blijkt dat drie hoofdkomponenten tezamen 83,0% van de totale va-riantle verklaren terwijl de overige hoofdkomponenten elk afzonderlijk minder dan 10% verklaren. Tevens ia in Tabel 4.3 te zien dat de eigenwaarde van fac-tor K kleiner dan 1 is. Vaak kiest men een. waarde 1 als drempel voor het al dan niet opnemen van factoren in het model. Op grond hiervan lijkt een keuze voor drie factoren redelijk.
Tabel 4.5 laat zien dat bij deze keuze ten minste 68% van de variantie van el-ke afzonderlijel-ke waterkwaliteitsparameter verklaard wordt (Npart), terwijl in een enkel geval zelfs sprake is van 93% (DOC).
Met behulp van Tabel 4.6 zijn waterkwaliteitsparameters met factoren te asso-ciëren. Op factor 1 scoren hoog: NH4, NO2/N03, P04 en Si (positief), alsmede Cl (negatief). Op factor 2 scoren hoog: Npart, Ppart en Cpart, terwijl DOC op factor 3 hoog 3cort.
Een praktische interpretatie is de volgende: bevat een waterkwaliteitamonster veel 'factor 1', dan betreft het relatief zoet water met veel nutriënten» Be-vat een monster veel 'factor 2', dan la het aandeel van de particulaire para-meters (Npart, Ppart, Cpart) hoog. De waterkwaliteitspararaeter DOC blijkt niet gecorreleerd te zijn met de overige beschouwde parameters (zie ook Tabel 4.2) en wordt als aparte factor onderschelden ('factor 3*).
Het voorkomen van factor 1 hangt sterk samen met het afvoerregime van de grote rivieren. De factoren 2 en 3 vertonen geen periodiciteit.
Conclusie ia dat het 3-factorraodel als beschrijving van de samenhang van de 9 waterkwaliteltsparameters zowel vanuit statistisch standpunt gezien goed vol-doet, als vanuit praktisch gezichtspunt goed interpreteerbaar is.
5* Ruimtelijke correlatie
5.1 Inleiding
In dit hoofdstuk wordt verslag gedaan van het onderzoek naar de ruimtelijk
correlatlepatronen van de waterkwaliteitsparameters. Deze correlatiepatronen
zijn van belang voor het evalueren van de representativiteit van meetlokaties.
De correlatiepatronen kunnen worden gebruikt bij het optimaal inrichten van
het meetnet, d.w.z. bij het kiezen van het optimaal aantal meetstations en de
lokatie hiervan. In plaats van op ruimtelijke correlatie kan men zich ook
richten op de hiermee samenhangende ruimtelijke variabiliteit. De resultaten
van dit onderzoek geven tevens inzicht op welke afstand-schaal belangrijke
va-riaties in de waarde van waterkwaliteitsparameters kunnen optreden.
Het onderzoek ia uitgevoerd met de interpolatiemethode. De
Kriging-interpolator is de beste lineaire zuivere Kriging-interpolator In de zin der kleinste
kwadraten.
Kriging levert over het gehele gebied naast de geinterpoleerde waarden ook een
schatting van 4e hierbij behorende standaardafwijkingen. Hiermee zijn vragen
naar representativiteit van meetlokaties te beantwoorden. In paragraaf 5.2
wordt nader ingegaan op de theorie en het praktijkgebruik van Kriging.
Voor het beschrijven vande ruimtelijke correlatie en variabiliteit van 9
waterkwallteitparametars: NH4, NO2/NO3, NTGT, P04, PTOT, Si, Cl, CPART, DOC,
ia gekozen voor een viertal karakteristieke tijdstippen, zie tabel 5.1.
nr tijdstip datum wind
rivier-afvoeren
2 17 feb 81 NO 3 laag
12 8 dec 81 N W S Haringvliet+Schelde hoog
15 9 feb 82 ZW 4 Haringvliet hoog
21 20 juli 82 NO 5 laag
Bij het maken van de keuze hebben meegespeeld:
- variatie in windrichting en -snelheid - variatie in niveau van de afvoeren - spreiding in tijdstippen over het jaar
- volledigheid van het bestand van meetgegevens - bemonsteringen rand tijkentering uitgevoerd
Paragraaf 5.3 bevat de beschrijving van de ruimtelijke correlatie en variabi-liteit en paragraaf 5.4 bevat Interpretaties en conclusies.
De representativiteit van de HISTOS-lokaties komt aan de orde in paragraaf 5.5. Het onderzoek heeft zich beperkt tot de op de HISTOS-meetpalen gemeten waterkwaliteltparatneter, i.e. Cl.
Interpretatie van de resultaten en conclusies volgen in par. 5.6. Het gebruik van een verdunningslijn wordt onderzocht in par. 5.7. Interpretatie van de re-sultaten en conclusies sluiten dit hoofdstuk af.
5.2. Beschrijving van de Kriging-techniek
5.2.1 Mathematische formulering van het probleem
We beschouwen de grootheid waar we ons voor interesseren, de waterkwaliteits-parameter z, als een geregionaliseerde stochastische funktie (g.s.f.)
£ C^i»C2)• Hierbij zijn £. en §„ de coördinaten in een twee-dimensionaal ge-bied* z in het punt x « (5,,^) noteren we verkort als z(x).
De g.a.f. z(x) heeft een verwachting E z(x) =» u(x).
Voorts interesseren we ons voor de correlatiestruktuur van de g.s.f. Deze kan gekarakteriseerd worden door de covariantie (mits deze bestaat):
cov [ssCXjhaCXj)] E f z ^ ) » ^ )][jj<x
Ook kunnen we gebruik maken van het semivariogrant gedefinieerd door:
Het probleem ia nu het definiëren van een schatter Z(XQ) voor de g.s.f. in een willekeurig punt X Q in termen van de beschikbare waarnemingen van a(x) in de
punten {x.,x„,...x }. We beperken ons tot schatters die een lineaire korabina-tie van de funktlewaarden in de meetpunten aijn:
dus tot schattingen:
l
0 (5.2)
Tevens hebben we behoefte aan de nauwkeurigheid van deze schatter, oftewel aan var z(xQ).
We stellen de volgende eisen aan de schatter:
var[z(x
Q) -
Z ( X Q ) ]is minimaal (5.4)
Het oplossen van dit probleem verloopt in twee stappen;
1 De correlatiestruktuur van de g.s.f. wordt geschat en gekarakteriseerd met
behulp van cov(x
1,x
2), Y<Xj,x
2) of een nog nader te introduceren funktle.
2 De optimale schatter z(x
Q) wordt uitgedrukt in de meetpunten { x ^ . . . x^}
met behulp van deze correlatiefunktie. Hierbij wordt, var s(x
Q) rechtstreeks
uit deze correlatiefunktie berekend.
Een praktische vraag is hoe men op basis van êên realisatie van de g.s.f.
(me-tingen z(x) op de stations {x ,... x } ten tijde t) kan vaststellen of aan de
zuiverheidskonditie (5.3) voldaan is. Hiertoe veronderstelt men dat de g.s.f.
ergodisch is. Dit betekent dat men het onderzoek naar het gemiddelde gedrag
van schatters over meerdere realisaties mag vervangen door onderzoek naar het
gemiddeld gedrag over êên (de beschikbare) realisatie.
Het zal blijken dat de met bovenstaande aanpak gekonstrueerde schatters exakte
interpolatoren zijn, dat wil zeggen z(x
p) - z(x
p) als x
peen meetpunt is.
5.2.2 Krlgen bij afwezigheid van trend
Voor het beschrijven van de oplossingstechniek, het "Krigen", bekijken we
eer-ste het geval dat de g.s.f. tweede-orde stationair is. dat wil zeggen:
1) u(x) hangt niet van x af
2) cov(xj
(X2) bestaat, en is een funktle van alleen de verschilvektor
h " xj-x2*
De tweede eis zwakken we af tot de zogenaamde intrinsieke hypothese;
2') Y ( X ,X ) is eindig en een funktle van alleen de verschilvektor h •» Xi-x,.
Stap 1 van het bepalen van de oplossing van het probleem bestaat uit het
schatten van de correlatiestruktuur van de g.s.f. Dit kan door het
experimen-teel semivariogram Y(h) te berekenen met behulp van
waarbij z(x^+h) en z(x^) meetpunten op afstand h zijn en N(h) het aantal pun-tenparen in de afstandsklasse h is.
Dit experimenteel semivarlogram wordt vervolgens gesmoothed en geparametrl-seerd, met als resultaat y(h).
In het vervolg van deze paragraaf veronderstellen we y(h) bekend.
Stap 2 in de oplossingsprocedure bestaat uit het berekenen van y... y bij gegeven y(h) zö dat
£(
XQ) " ï \ f^*.») (5.6)
voldoet aan (5.3) en (5.4).
Als de intrinsieke hypothese van kracht is geldt onder de voorwaarde dat {A^i-1,.,,11} een kontrast is, dat wil zeggen
l \ ~ 0 (5.7)
dat n n nvar l X z(x ) « - l l y(x -x ) (5.8)
i»l X 1-1 j-1 X 2 We onderzoeken de schattingsfout z(x0) - z(xQ);var[z(x ) - z(x )] - var[ \ X z(x ) - z(x )j - var J X z(x ) (5.9)
i-1 u i-O 1 x waarbij X • -1 Uit (5.3) volgt; n
E[ l X z(x ) - z(x
n)] » 0
i-1
x i°
zodat n E l x. z(x,) « 0 i-O 1 ±en dus
Y X, - O (5.10) 1=0 1
Hieruit volgt dat [x |i = 0,...n} een kontrast is, en daarmee is (5.9) uit te werken tot;
var[|(x0) - z(x0)] = - \ \ Y ^ - X j ) + 2 I X± Y C ^ - X Q ) (5.11)
Het optlraaliseringsprobleem is nu te formuleren als: zoek die {X ji => l,..,n} waarbij (5.11) minimaal is onder voorwaarde (5.10).
Deze {X |i =• l,...n} zijn te vinden door gebruik te maken van de Lagrange-multiplicator u, en het volgende stelsel op te lossen:
n
l X y(x -x.) + u - Y ( X . - X Q ) (i=l,.,,n) (5.12)
IX, =1 (5.13)
1-1
Dit stelsel {(5.12),(5.13)} vormt het zogenaamde Kriging-stelsel.
5,2.3 Aanwezigheid van een ruimtelijke trend
De in de vorige paragraaf ontwikkelde oplossingmethode is niet toepasbaar als er een ruimtelijke trend aanwezig is in de verwachting van de g.s.f. z(x). In dit geval geldt niet meer dat y(x) - v, onafhankelijk van x, en dit is een
es-sentiële voorwaarden voor het rekenen onder de intrinsieke hypothese.
Een voor de hand liggende aanpak is het uit de data schatten van de ruimtelij-ke trend, het reduceren van de data met deze trend en het hierna toepassen van de procedure van paragraaf 5.2.2. Bij de schatting 2(XQ) wordt uiteindelijk de voor XQ gevonden trend-waarde opgeteld.
In de praktijk werkt deae aanpak niet. Het is bekend dat het op deze wijze verkregen experimenteel variogram van residuen een onzuivere schatting Is voor het "echte", onderliggende varlogram (Delhomme, 1978).
