Wavelet representation of localized multiscale flow structures in high-resolution variational data assimilation ´ Fournier,∗a Mark Dixon,b Ross N. Bannistera & Sue P. Ballardb Aime a U. Reading Dept. Meteorology, b UKMO Joint Centre for Mesoscale Meteorology September 22, 2008
Abstract This work aims at improved computation of multiscale structures and covariances in the UK Met Office limited-area data assimilation system, in particular the Unified Model (UM) control-variable transform. Such structures are exhibited in samples of the covariance matrices computed from UM forecast increments. Three-dimensional fields of horizontal velocity are Helmholtz decomposed to obtain 3D velocity-potential χ and streamfunction ψ increment fields. The Wishart distribution implies that estimating a rank-Ns representation of the full covariance matrix requires ≥ Ns + 2 independent field samples, which would be prohibitive at full rank. Rank reduction is achieved in the vertical by standard structure functions and in the horizontal by special symmetric (linear phase response) biorthonormal wavelets that include boundary corrections at every scale and can be applied to non-dyadic data lengths. Estimates of the time-average and standard-deviation χ and ψ fields for the leading 4 vertical structures κ ≤ 4 appear to be reasonable despite too few samples being available. The wavelet compression of 9 spatial samples of ψ-covariance matrix for κ = 1 explains > 98% of the covariance with < 2% of the matrix entries. Compression is achieved using the simplest possible criteria, so the possibility of much greater compression is strongly suggested. Methods to construct inhomogeneous anisotropic covariance models are cited, and other potential technical improvements are discussed.
1
Introduction
1.1
Purpose: Joint scale & location dependence of covariance
This work aims at improved computation of multiscale structures and covariances in the UK Met Office (MO) limited-area data assimilation system (in particular the Unified Model control-variable transform), in line with MO plans [39]. As NWP codes and observational data both increase in spatiotemporal resolution and complexity, data assimilation (DA) must better represent and compute strongly nonlinear multiscale structures in the atmosphere, and their covariances. Although in this work the background covariance is considered to be time-independent, its scale dependence is itself location dependent, a feature that should greatly improve the overall model flow dependence.
1.2
Background: Multiresolution analysis of nonlinear multiscale structures
The present research draws upon a background in high-order dynamically [26, 27, 50] and statically [24, 28] adaptive refinement methods and multiresolution analysis (MRA) methods [20, 21, 22, 23, 25] for analysis and simulation of these kinds of structures. Several especially relevant meetings [e.g., 36, 44, 56], reports [e.g., 1, 3, 4, 13, 14, 15, 34, 38, 40, 41, 42, 43, 48] and projects [e.g., 2, 12] provide recent context for this project. Several results from MRA are immediately suggestive for improving DA.1 MRA is a special kind of scale-selection method [cf. 37] that provides appropriate customized projection operators. Moreover, it has long been known in the wavelet-analysis community, that covariance and similar 2-point functions that can be trusted to decay as a function of point-separation no more slowly than a certain form of estimate, are guaranteed to yield sparse wavelet representations with a factor O[N log N ] (assymptotically ∗ Corresponding 1 After
author:
[email protected] a 2007/2/28 proposal to MO by A. Fournier.
1
O[N ]) significant entries per space coordinate discretized to N points [5 pp. 155, 160; 6 p. 10246a]. This Qd Pd implies O[ a=1 Na ] cost for d coordinates, but recent methods may reduce this to a=1 O[Na ] [6, 7]. Also there is an established record for wavelet-based representation of covariance, e.g., for a 143-member ensemble of N1 N2 = 128 × 128-pixel face images the “best basis” wavelets [e.g., 20] represent 90% of the variance with 2% of the coefficients, and after appending a wavelet-approximated Karhunen-Lo`eve analysis, represent 88% with 0.06% [54].
2
Data: Seven hours of 1-km MO Unified Model output
We used a 3D data set of horizontal wind ~u[t, λn1 , ϕn2 , zl ] = ~un,l = un,l~en1 + vn,l~en2 predicted from the Nt :=#T = 6 hourly times t ∈ T:={2005/7/4 9:00, . . . 14:00} plus the verifying analysis at t = ta :=15:00, output from the MO Unified Model v6.1, on Nq := N1 N2 = 360 × 287 uniformly gridded horizontal spherical coordinates rotated [52, §2.1.3] to geographical origin {−2.5◦ , 52.5◦ }, {λn1 , ϕn2 } ∈ X = [λ1 , λN1 ] × [ϕ1 , ϕN2 ] = [−2.193◦ , 2.653◦ ] × [−2.623◦ , 1.238◦ ] S (with boundary normal ~ν n = ±~ena defined for n ∈ B := {1, . . . N1 } × {1, N2 } {1, N1 } × {1, . . . N2 }), and N3 = 76 non-uniform hybrid height coordinates [52, §B.4] zl ∈ [z1 , zN3 ] = [2.5m, 37.14km].
(1)
Each prediction and analysis was computed using the same inhomogeneous lateral boundary conditions, so the t − ta increments 4~un,l [t] := ~un,l [t] − ~un,l [ta ] (2) satisfy the homogeneous Dirichlet (impermeable) boundary condition ~ν n · 4~un,l = 0
3
∀n ∈ B.
(3)
Analytical methods
3.1
Velocity-potential and streamfunction computation
To perform the Helmholtz decomposition into the potential χn,l and streamfunction ψn,l , ~ n χl + ⊥ ∇ ~ n ψl 4~un,l = ∇
~ n · 4~ul χn,l = (∇−2 ∇) , −2⊥ ~ ∇)n · 4~ul ψn,l = (∇
(4)
and 4v l = (CN1 ⊗ SN2 )ˆ v l ∈ RNq ,
(5)
(3)
⇐⇒
it’s useful to expand 4ul = (SN1 ⊗ CN2 )ˆ ul ∈ RNq where N N Cm,n := cN m cn
q
2 N −1
cos π(m−1)(n−1) N −1
N and Sm,n :=
q
2 N −1
sin π(m−1)(n−1) N −1
(6)
denote the symmetric orthonormal cosine- and pseudo-orthonormal sine-Fourier bases: m = m, ´ 0, N m−1 ´ 2 π N sin (m+ m) ´ sin π CN CN = IN , SN SN = 0 ⊕ IN −2 ⊕ 0, †C N S = cm m m ´ 2 N −1 , m= 6 m ´ N −1 sin π m+m−2 ´ π m−m ´ sin 2
N −1
2 N −1
−1/2 and cN . Then (4) becomes m := (1 + δ1,m + δN,m )
ˆ l ∈ RNq , ψ l = (SN1 ⊗ SN2 )ψ m1 Ξ1 , m2 Ξ2 {χ ˆm,l , ψˆm,l } = −(|Ξ||m|)−2 {ˆ um,l , vˆm,l } m2 Ξ2 , −m1 Ξ1 ˆ l ∈ RNq , χl = (CN1 ⊗ CN2 )χ
(7)
where
and Ξ := {(λN1 − λ1 )−1 , (ϕN2 − ϕ1 )−1 }π
2
is a 2D fundamental wavenumber.
(8)
3.2
Covariance-matrix estimation
The unbiased estimate of the full covariance matrix Be ∈ RNx ×Nx is
where
B = (1 − Nt−1 )−1 (hψ †ψi − hψi †hψi) ∈ RNx ×Nx , X hAi := Nt−1 A[t]
(9) (10)
t∈T
is the sample average and Nx := N1 N2 N3 ≈ 8 × 106 is the total spatial data size. If independently for each t, ψ[t] were Nx -variate normally distributed with exact covariance matrix Be , it can be shown [∵ 32 pp. 90, 92, 98; 51 p. 24] that B−1/2 BB−1/2 is Wishart distributed with mean INx , e e Nt −Nx −2 Nx mode Nt −1 I 1 and variance Nt −1 (INx + 1);
(11)
however Nt ≥ Nx + 2 would be required, which is prohibitive. Also, these (1 + Nx )Nx /2 ≈ 3 × 1013 entries are far too numerous to compute. To follow good statistical practice we should transform from the full ψ ∈ RNx to a sparse representation Qψ ∈ RNs of rank Ns ≤ Nt − 2 Nx . Then (11) becomes [∵ 32, p. 96] (QBe †Q)−1/2 QB †Q(QBe †Q)−1/2 , Wishart distributed with mean INs , Nt −Ns −2 Ns mode Nt −1 I 1 and variance Nt −1 (INs + 1),
(12)
which still decreases slowly as Nt increases. Note that QB †Q is s.p.d.2 Part of our goal is to build Q so that (12) is useful.
3.3
Horizontal 2D symmetrized biorthonormal wavelet analysis
In place of the common global spectral analysis using sinusoids, a suite of biorthonormal wavelet analyses c [19]. An were applied. These were implemented by the free on-line MATLABr -library, WaveLab introduction to orthonormal wavelet analysis focusing on meteorological applications, including extensive c demonstrations of WaveLab, was provided by Fournier [20]. Orthonormal wavelets are inferior to biorthonormal ones in respect to necessarily possessing some asymmetries and hence nonlinear phase responses [11 pp. 492, 535; 54 pp. 164, 166]. The method used here possesses no asymmetry and also minimizes boundary artifacts by using particular extension operators E at every scale and spatial interval of the analysis and synthesis [9, 10]. The method used here also eliminates the usual severe restriction to a dyadic number N = 2J ∈ {2, 4, 8, . . .} of points along each coordinate. Biorthonormal wavelet analysis is a kind of frame analysis [e.g., 54, pp. 14, 156] that in 1D involves transforming a sequence rn of any number N of function values uniformly spatially sampled at coordinates λn , into a sequence of N scaling-function and wavelet coefficients r˜i :=
N X
?N Wn,i rn ,
(13)
n=1
and back again, rn =
N X
N Wn,i r˜i .
(14)
i=1
Whereas the index n only describes spatial location, the index i ∈ {ij−1 + 1, . . . ij } describes both wavenumber bandcenter and bandwidth ∝ ij and spatial location ∝ i/ij (both as precisely as the Heisenberg principle allows [cf. 54, p. 24]), where ij−1 := dij /2e =: ij − ˜ıj−1 ∈ {1, . . . iJ−1 } and J := dlog2 N e
defines a dyadic partition of
N =: iJ ,
is the maximum level.
(15) (16)
c The WaveLab commands that perform (14) and many other calculations in this section are schematized in Table 1. 2 Proof:
although Q is singular, B = †B is s.p.d. ⇔ B = A †A (∃A) ⇒ QB †Q = QA †(QA) ⇔ QB †Q is s.p.d.
3
c Table 1: WaveLab functions. The labelling of matching the literature. Eq. Expression † 13 r˜ := †rW?N 14 r = WN r˜ 15 {i0 , . . . iJ } 16 [N, J] ?N ?N ?N ?N † N 19,20 I n {V 1 , . . . V Jm , W 1 , . . . W N −Jm } ? 19 VN i 20? WN i ¯ 21 ↓j y † 22 yXj [f ] † 23 r¯ E+ n,0,2 ,3 24 h?
primal and dual is self-consistent but not necessarily Function call = FWT SBS(†r, Jm , l? , l) = IWT SBS(˜ r , Jm , l? , l) = dyadpartition(N ) = dyadlength(r) ? = FWT SBS(†I N n , Jm , l , l) ? N = IWT SBS(I i , Jm , l , l) ? = IWT SBS(I N im +i , Jm , l , l) ¯) = † UpSampleN( †y = symm aconv(f , †y) = symm iconv(fliplr(f ), †y) = extend(†r¯ , 1 + 2 , 1 + 3 ) = MirrorSymmFilt(fliplr(l)) lshift(· · · lshift( †y)), p ≥ 0 = rshift(· · · rshift( †y)), p ≤ 0 ←−−−−−−−−−−−−→
yTpj = †y(T1j )p
25
†
26 27 26,27 29 30 29,30
yXj [l? ]↓j 1 † yTj Xj [h? ]↓j [†r¯ L?j , †r¯ H?j ] † † ¯ ↓j X[l− ] y − † † ˜ ↓j T−1 y j X[h ] Lj β + Hj α [l? , l]
|p|
= DownDyadLo PBS( †y, l? ) = DownDyadHi PBS( †y, l) = DownDyad SBS(†r¯ , l? , l) ¯ , l) = UpDyadLo PBS( †y = UpDyadHi PBS(˜ y , l? ) = † UpDyad SBS( †β, †α, l? , l) = MakeBSFilter(Type, Par)
†
The primal and dual wavelet transforms WN and W?N are biorthonormal and complete: WN W?N = WN †W?N = IN ,
†
(17)
so reconstruction is accurate to machine precision and no condition number is increased. No full transform matrix need ever be stored, but †r˜ ≡ †rW?N and r = WN r˜ are computed by a recursive hierarchy on J − Jm levels; e.g., †r˜ by post-multiplying the row †r by W?N := {L?J {L?J−1 · · · {L?Jm +2 {L?Jm +1 , H?Jm +1 }, H?Jm +2 } · · · , H?J−1 }, H?J } ?N ?N ?N N ×N , = {V ?N 1 , . . . V im , W 1 , . . . W N −im } ∈ R
(18)
and similarly for WN , obtained by removing the ? everywhere. The minimum level Jm ∈ {0, . . . J} controls the trade-off between im := iJm ∈ {1, . . . N } scaling-functions and N − im wavelets in the analysis and synthesis, e.g., Jm = 0 implies N − 1 wavelets and a single scaling-function of uniform value ∀λn , and Jm = J implies W?N = IN (no wavelets). Equation 18 implies that the entry in row n, column i of {V
?N ?N 1 , . . . V im }
=
J−J m −1 Y
L?J−k ∈ RN ×im
(19)
k=0
approximates the value at λn of the ith dual scaling-function Vi? [λ] as a product of dual low-pass filter matrices L?j ∈ {0, l? }ij ×ij−1 ; similarly the entry in row n, column i − ij−1 of J−j−1 Y
?N {W ?N ij−1 +1 , . . . W ij } = (
L?J−k )H?j ∈ RN טıj−1
(20)
k=0
approximates the value at λn of the ith dual wavelet Wi? [λ] as a product of dual low-pass factors and one dual high-pass filter matrix H?j ∈ {0, ±l}ij טıj−1 . For example, Fig. 1 shows a certain primal basis 4
Table 2: New MATLABr commands. Also see header and other comment lines of individual M-files. Purpose Load raw variables Create vertical structure functions Fix b.c. (37) Compute χ & ψ (7) Spectrum of ψ (not shown) Apply F Plot certain F κ χ & F κ ψ (not shown) Load preliminary data Compute χ & ψ statistics Prepare 2D wavelet transform Load basis experiment Create 25 l & l? Create 25 WN1 & W?N1 and (17) errors Test reconstruction of norm ˜ and get FWT SBS timings Create 25 ψ Fig. 1 Fig. 2 Fig. 3 Fig. 4 Fig. 5 Figs. 6, 7 Sample true ψ covariance (33) Fig. 8 Fig. 9 Compute (33) ⇐ (35), Figs. 10 & 11 Figs. 12, 13, 14 & 15 Figs. 16 & 17 and get IWT SBS timings Plot sorted |ψ˜i | (not shown) Plot vertical means (not shown) Compute boundary errors Fast C transform (6) Do ’job’ if ’y’ to ’ ?’ Wavenumbers → spectrum (not shown) Declare global variables Wavenumber annuli (not shown) Spectral shells of ψ (not shown) Compute (5,7,8)
Command analysop(’load’) analysop(’modesL76’) modesL76 analysop(’fix bc’) analysop(’Helmholtz’) analysop(’psi anal’) analysop(’vtransform’) analysop(’vtransform plot’) prelstat(’load’) prelstat(’chpsvs’) prelstat(’WdWp’) basitest(’load’) basitest(’MakeBSFilter’) basitest(’WdWp’) basitest(’recon’) basitest(’wpsi’) basitest(’1plot’) basitest(’xplot’) basitest(’mplot’) basitest(’2Dplot’) analysop(’modes plot’) prelstat(’vmodeplot’) vmodeplot({λ1 , . . .}, {ϕ1 , . . .},chpsvs,{’data1’ ...}) prelstat(’covaanal’) covaanal(hψi,{’data1’ ...},n samp,v mod) prelstat(’corrplot’) corrplot(n samp,B t,st dψ ,i x,i y,{λ1 , . . .}, {ϕ1 , . . .}) prelstat(’covaplot’) covaplot(n samp,B t,i x,i y,{λ1 , . . .}, {ϕ1 , . . .}) prelstat(’psiw’) analysop(’u v exnerp’) basitest(’psir’) basitest(’sova’) prelstat(’vmeanplot’) vmeanplot({’data1’ ...}) bcerror({’data1’ ...},update period) ˆ = dct i(v) v evalinpu(’?’,’job’) Fwn2rps globaset(’u v w ...’) shelpass analysop(’u v exnerp shells’) [χ, ψ] = uv2vpsf(u, v, [π/Ξ1 , π/Ξ2 ])
5
CDF−3−9 wavelet
Wi[λ]
0.2 0 −0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
(λ−λ1)/(λN −λ1)
45
0.9
23
1
123
basitest.m
12
6
mode i
11−Aug−2008
Figure 1: Primal wavelets Wi [λ] from the spline-based l and l? with 3 and 9 vanishing moments respectively [11]. Each wavelet is offset ⊥ λ and colored ∝ i ∈ {1, . . . 45}, where ij ∈ {1, 2, 3, 6, 12, 23, 45, 90, 180, N1 = 360} is the dyadic partition. for Jm = 0. (The new MATLABr commands created for the figures and analyses in this work are summarized in Table 2.) No L?j or H?j need ever be stored except for a small number of entries from the primal and dual low-pass filter sequences l and l? of lengths 2P − 1 and 2P ? − 1. Each †r¯ L?j or †r¯ H?j is computed as a restriction n no + Rn := I0 ∈ {0, 1}ij ×n of a decimation (downsampling) +
+
+
↓j := Iij ⊗ {10 } ∈ {0, 1}2ij ×ij of a convolution Xj [f1 , . . . f2P −1 ] :=
P −1 X
(21) +
+
fP +p Tpj ∈ {0, f }2ij ×2ij
(22)
p=1−P
between a symmetrized extension [9, §I.A] {¯ r1 , . . . r¯n }Esn, := {¯ r1 , . . . r¯n , 0, . . . 0,s¯ rn+2 −1 , s¯ rn+2 −2 , . . . s¯ r3−3 , s¯ r2−3 } ←−−−→
(23)
1
∈ {0, r¯ , s¯ r }1×(2n+1 +2 +3 −2) of r¯ and a dual low- or high-pass filter sequence l? or h? := (−1)P +1 {l2P −1 , −l2P −2 , l2P −3 , . . . − l2 , l1 },
(24)
respectively, where ( Tpj
:=
0 + + I2ij −mod[p,2ij ]
+
Imod[p,2ij ] 0 6
) +
+
∈ {0, 1}2ij ×2ij
(25)
is a periodized translation, i+ j := ij + o − 1 describes the symmetry lp = l2P +o−p
about its centre P + o/2 ∈ {P, P + 21 }
of {l1+o , . . . l2P −1 }
and l? has the same symmetry type as l. (By convention o = 1 ⇒ l1 = l1? = 0.) Putting all this together, L?j = H?j
=
? ij−1 E+ , ij ,0,o,o Xj [l ]↓j R
(26)
1 ? ˜ ıj−1 . E+ ij ,0,o,o Tj Xj [h ]↓j R
(27)
These low- and high-pass filters are inverted using biorthogonal filters: {L?j , H?j } †{Lj , Hj } = †{Lj , Hj }{L?j , H?j } = Iij , where
(28) − ij † E+ ij−1 ,0,1−oj ,o ↓j Xj [l ]R
†
Lj =
(−1)
∈ {0, l}
ij−1 ×ij
o
, ? ˜ij−1 ×ij
− ij Hj = E˜ıj−1 ,ooj ,o+oj −ooj ,1 †↓j T−1 j Xj [h ]R ∈ {0, ±l }
†
(29) (30)
− (noting both these E factors have i+ j columns), {f1 , . . . fn } := {fn , . . . f1 } and oj := ij (mod 2). The ? 2b(P + P )/2c − 1 equations min[P ? −1,P −1−2q]
X
lP? ? +p lP +p+2q = δ0,q ,
q ∈ {d1 − (P + P ? )/2e, . . . b(P + P ? )/2 − 1c},
p=max[1−P ? ,1−P −2q]
combined with eq. 24 (and its primal version) provide necessary (to within 2 free parameters not shown here, and assuming finite P and P ? ) and sufficient conditions for (28) and thence (17), and admit many kinds of solutions e.g., using Fourier series [11, p. 495]. c For l and l? , WaveLab provides many options, including interpolating [17], average-interpolating [18], spline-based [11, Table 6.1] and image-compression “optimal” [53] filters. Figures 2 & 3 show one primal and one dual wavelet from each of the 25 wavelet transforms considered in this work, vs space and wavenumber coordinates. In d dimensions (13,14) are generalized to tensor-product transforms between columns of spatial values or wavelet coefficients r, r˜ ∈ RN1 ···Nd : r. r = (WN1 ⊗ · · · ⊗ WNd )˜
r˜ := †(W?N1 ⊗ · · · ⊗ W?Nd )r,
(31)
Figure 4 shows 2D primal wavelets with fixed i2 and successively halving i1 values.
3.4
3D orthogonal decomposition and covariance rank reduction
As described elsewhere [8, eq. 1.17], an orthogonal vertical-structure-function matrix F = {F 1 , . . . F N3 } ∈ RN3 ×N3 is provided by diagonalizing a weighted vertical-covariance matrix. The columns F κ ∈ RN3 are plotted in Fig. 5. The full 3D biorthonormal transforms between r, r˜ ∈ RNx are a further modification of (31): r˜ := †(W? ⊗F)r,
r = (W ⊗F)˜ r,
W? :=W?N1 ⊗W?N2 ∈ RNq ×Nq ,
W :=WN1 ⊗WN2 ∈ RNq ×Nq . (32)
A first step in reducing the rank of the covariance-matrix estimate (9) is to consider the vertically transformed B horizontally sampled at an arbitrary reference point λn := {λn1 , ϕn2 }, and decompose it using (32): ˜ n,κ B n,κ := †Fκ BFκ I n =: WB p p =: †I n B n,κ diag( †Fκ BFκ ) ρn,κ ∈ RNq ,
(33) (34)
N2 Nq 1 where Fκ := INq ⊗ F κ ∈ RNx ×Nq is the 3D vertical transform, I n := I N is the horizontal n1 ⊗ I n2 ∈ {0, 1} Nq sampler and ρn,κ ∈ [−1, 1] is the correlation between λn and every other point. Using the intrinsic decorrelating properties of wavelets [op. cit.] and F, (33) can be approximated using a sparse columnmatrix ( I i , Mi is true, s † Nq Nq ×Ns ˜ ˜ B n,κ := M MB n,κ ∈ R , where M := {M 1 , . . . M Ns } ∈ {0, 1} , M k := (35) 0, otherwise.
7
Ave−Interp−2
0
Deslauriers−7
0.5
Deslauriers−5
Deslauriers−3
Deslauriers−1
25 primal & dual wavelets at indexes 16 & 21 1
−0.5
−0.5
CDF−1−5
0
CDF−1−3
0.5
CDF−1−1
Ave−Interp−6
Ave−Interp−4
1
−1
CDF−3−1
0
CDF−2−8
0.5
CDF−2−6
CDF−2−4
CDF−2−2
1
−0.5
CDF−4−4
−0.5
CDF−3−9
0
CDF−3−7
0.5
CDF−3−5
CDF−3−3
1
−1
Villasenor−5
0
Villasenor−4
0.5
Villasenor−3
Villasenor−2
Villasenor−1
1
−0.5 0
0.5
0
0.5
0
0.5
0
0.5
0
0.5
(λ−λ1)/(λN −λ1) (λ−λ1)/(λN?−λ1) (λ−λ1)/(λN −λ1) −λ1) )/(λN −λ1) c N 1 Figure 2: W1 16 [λ] (blue) and W21 [λ] (green) vs λ for1 25 filters (λ−λ from1)/(λ WaveLab as (λ−λ labeled. 1 1 1
8
16 1024
2
524288
1024
524288
65536
65536
8192 1024
8192 1024 128
CDF−3−7
524288
CDF−3−5
16
524288 65536 8192 1024 128 16 2 0.25
524288 65536 8192 1024 128 16 2 0.25 0.0312 0.0039
16
2
524288
524288
524288
65536
65536
65536
8192
8192
1024 128
1024 128
Villasenor−3
16
8192
16 16 1 4 16 64 1+wavenumber m
8192 1024
Ave−Interp−2 CDF−1−5
65536
128
Villasenor−2
Villasenor−1
CDF−3−3
128 1024
1024
524288
2 1 4 16 64 1+wavenumber m
1024 128 16
524288 65536 8192 1024 128 16 2 0.25 0.0312 0.0039 0.0005
524288 65536 8192 1024 128 16 2 0.25 0.0312 0.0039 0.0005 0.0001
2 1 4 16 64 1+wavenumber m
8192
1024 128 2
524288 65536 8192 1024 128 16
524288 65536 8192 1024 128 16
524288
65536
8192
16
128
CDF−2−6
CDF−2−4
CDF−2−2
8192
8192
65536
524288
8192
524288 65536
65536 65536
65536
128
CDF−3−1
16
524288
CDF−4−4
128
524288
524288 65536 8192 1024 128 16 2 0.25 0.0312 0.0039 0.0005 0.0001
524288
Villasenor−5
1024
524288 65536 8192 1024 128 16 2 0.25 0.0312 0.0039
CDF−1−3
8192
Ave−Interp−6
Ave−Interp−4
65536
CDF−1−1
524288
CDF−2−8
128
CDF−3−9
1024
524288 65536 8192 1024 128 16 2 0.25 0.0312
Villasenor−4
8192
8192
Deslauriers−5
65536
Deslauriers−3
Deslauriers−1
65536
Deslauriers−7
25 primal (blue) & dual (green) wavelet spectra at indexes 72 & 72 524288 524288
1024 1 4 16 64 1+wavenumber m
65536 8192 1024 128 1 4 16 64 1+wavenumber m
Figure 3: As in Fig. 2 except spectra (scaled to max at 219 ) vs 1-plus-wavenumber. Each ordinate ? lower-limit is the lesser of the spectrum values at wavenumber 1 of W16 and W21 . Labelling the same ordinates in each plot illustrates the different bandwidth decay rates.
9
0
0.20
0.8
W64[x]W16 W[y] [x]W16[y] 64
0.6
0.4
0
0.2
0 0.20
0.4 0.2
0.8
0.4
0.6
0.8
W8[x]W16[y] 0.6 0.4 0.8 0.6
0.8
W16[x]W16 W[y] [x]W16[y] 16
0.6 0.8
0.8
0.6
0.4
0 0.2
0.2
0.4
0.6
0.8
0.8 0.6
0.8
0
0.8
0
0.2
0.6
0.8
0.8
0.8
0.4 0.2
0.6 0.4
0.8 0.6
0.8
0.2
W2[x]W16W [y]2[x]W16[y] 0.4 0.6 0.8
0
0.6 0.4
0.8 0.6
0.8 0.4 0.6 0.2 0.4
0
0.8
4
16 4
W [x]W W [y][x]W [y]
16
2
16 2
0 0.20
0 0.2
0.4 0.2
0.6 0.4
0.8 0.6
0.8
0.2 0
3.2
0.8
0.8
3 0.6
0.6
0.6
0.4
0.4
0.2
0.2
0.2 2.8 0.4
0.4
2.6
0
0 2.4 0.2
0.2
0.2
0.4
0.6
0.8
level (m) 0
11−Aug−2008
0 20.2 0
0.4 0.2
0.6 0.4
0.8 0.6
0
0.8
0 0.20
0
1.8
basitest.mbasitest.m 1.6
level (m)
2.2 0.4 0.2
3.4
3.2
3.2
3
3
2.8
2.8
2.6
2.6
2.4
2.4
2.2
2.2 0.8 0.6
0.8
1.8 1.6
11−Aug−2008 11−Aug−2008
2 1.8 1.6
1.4
1.4
1.4
1.2
1.2
1.2
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
−0.5 −0.25
−0.5 −0.25
0 0.25 0.5 60
70
80
90
100
−0.5 −0.25 0 0.25 0.5 50
60
70
80
90
100
P⋅eχ variance % P⋅eψ variance % Figure 5: Vertical structure functions −0.7 < Fl,κcorranal.m < 0.7 vs task zl (1) and cumulative percentage of variance 24−Apr−2008 plotmodes explained up to index κ, for χ (left) and ψ (right).
10
4
3.6
3.4
20.4 0.6
0
x 10
basitest.m
16
3.6
0
W4[x]W16[y]
0.4
0.4 0.4 0.2 0.2
0.4 0.2
0.6
0.6
W [x]W W [y][x]W [y]
0.8
0.6
e%
0 0.2
0.6
0.4
90
0.4 0 0.20
0.6
0.6 3.4
100
0.6
Figure 4: As in Fig. 1 except 2D primal wavelets W27−j [λ]W16 [ϕ] vs πΞ1 (λ − λ1 ) & πΞ2 (ϕ −11−Aug−2008 ϕ1 ) with 11−Aug−2008 W [x]W [y] basitest.mbasitest.m 0 0 0 0 2 16 0 0.2 0.4 0.6 0.8 0 . .0.4 0.4 0.8 0.6 0.8 0 0.20 0.4 0.2 0.6 0.4 0.8 0.6 0.8 ∈ 0.2 {1, .0.2 6}. 0.6 4 4 x 10 x 10 3.6
0
0.8 0.6
0
W16[x]W16[y]
0j
0.8
0.8
0.8
0.2
0
0
0.2
0.8
0 0.20 0.2
W8[x]W16W [y]8[x]W16[y] 0.4 0.6 0.8
0.8
0
W8[x]W16W [y]8[x]W16[y]
0.4 0
0.6 0.4
W40[x]W16W [y]4[x]W16[y] 0.20 0.4 0.2 0.6 0.4
0.6
0
0.2
0.4 0.2 0.2
0.2
0.4
0.8
0.2
0.4
0.6
0.8
0.4
0.4
0.6 0.2
0.6
0.2 0.4
0
0.8
0.8 0.6
0.6
0.6
0.8
0.4
0.4
0
0
0.8
0.6
0.8
0
0.8
0.6 0.4
level (m)
0
0 0.2
0.2
0.4 0.2
0.8
0.2 0 0.20 0.2
0.20
0
0.4
0.4 0
0
0.8
0.6 0.2
0.4
0.2
0 0.2
0.4
0.6
0
0.8
W32[x]W16 W[y] [x]W16[y] 32
0.8
0.2 0.4
0.8 0.6
0.6
0.8 0.4
0.8
0.8
0.6
0.6 0.8
0.6 0.4
W16[x]W16 W[y] [x]W16[y] 16
W32[x]W16[y]
0.8
0.4 0.2
0 0.25 0.
P⋅eχ
Table 3: Qualities (eq. 36 in %) of sparse representation (33) using (35), at 9 reference points. (£M would be n1 - as well as n2 -symmetric but for a small artifact of the λn1 generation.) £M εM ◦ ◦ ◦ ◦ λn1 = −1.40 0.22 1.84 −1.40 0.22◦ 1.84◦ ϕn2 = 0.60◦ 98.49 98.34 98.46 1.33 0.20 0.70 −0.69◦ 98.41 98.26 98.37 0.78 0.10 0.36 ◦ 98.49 98.34 98.46 1.67 0.15 0.67 −1.99
That is, M †M ∈ {0, 1}Nq ×Nq is a diagonal mask matrix of rank Ns Nq for a judicious criterion Mi . Referring to (12), since n and κ are arbitrary, we may identify Q := †(Fκ W? M) ∈ RNs ×Nx of rank Ns . Then Fκ WMQB †(Fκ WMQ) is the rank-Ns reconstruction of B. We use the criterion “WnN11,i1 WnN22,i2 6= 0” for each n in (33). That is, only 2D primal wavelets whose support contains λn are included. The error and memory savings can be measured by εM := 1 − corrcoef(B n,κ , B sn,κ )
and
£M := 1 − Ns /Nq ,
(36)
respectively, where corrcoef is the usual correlation coefficient. Many other kinds of mask matrices are discussed in the literature, especially thresholding [e.g., 35].
4
Results
As an initial assessment of the data, consider the time average (10) and standard-deviation fields of χ and ψ shown in the top two rows of Figs. 6 & 7, respectively, all projected on the 4 F κ that explain the most vertical variance. The lower Nt rows show the instananeous fields. Although Nt is small, the ratio of standard deviation to time average appears to be reasonable. For the most part, the instananeous fields are similar to each other but of decreasing amplitude and increasing smoothness as t increases. This trend is consistent with χ and ψ being computed from the increments defined by eq. 2. Exceptions to this trend may well be due to data assimilated inbetween certain forecast lead times, which is not the concern of the present investigation. Figure 8 shows the correlation function ρn,1 (34). The contour levels are chosen to accumulate near |ρn,1 | = 1 in order to accentuate the inhomogeneous anisotropic structures, of which a better understanding and model are the main concerns of this study. Figure 9 shows the covariance (33) for 9 different reference points λn . Four of the λn exhibit much greater covariance ranges than the other five. The shape of the Fig.-9 contours illustrates the inhomogeneity and anisotropy of the covariance. For sparse-representation experiments we used the CDF-3-9 wavelet shown in Figs. 1 (1D), 4 (2D) and 2 & 3 (row 4 column 3). For all experiments Jm = 0. Figure 10 shows the sparse-representation covariance, and Fig. 11 shows the difference between the fields in Figs. 9 & 10. The quality of the representation is summarized in Table 3.
5
Concluding remarks
We have presented an anisotropic and inhomogeneous covariance model based on severely truncated wavelet expansions. Some other methods to construct covariance models are discussed or reviewed elsewhere [e.g., 29 p. 725; 30 p. 2450; 33; 47]. In future work we consider the following. 1. Multiple dependent fields may be mutually orthogonally analyzed e.g., using multiwavelets [e.g., 26, and citations therein]. 2. Also, nesting a limited-area model (LAM) into a global model (GM) raises issues of lateral boundary conditions (BCs) [14] [cf. 31, 55]. These issues may be addressed by constructing the generalized low-pass projection ℘ from the LAM onto the GM subdomain it replaces, and its complement high-pass projection 1 − ℘, since in ran ℘ no BC issues arise. 3. DA requires many repeated convolutions with non-local and/or non-separable kernels. These could be accelerated using appropriate local, separable operators [cf. 6, 46, 49]. 11
† Figure p 6: Column κ shows surface plots vs λn1 & ϕn2 of: Fκ hχi (row 1); (1 − Nt−1 )−1 h †Fκ χ †Fκ χi − †Fκ hχi †Fκ hχi (row 2); †Fκ χ[t + kh] for t = 9:00 (row k + 3 > 2). The color map is fixed for each κ.
12
Figure 7: As in Fig. 6 except for ψ instead of χ.
13
−1
−1 2
2 1/2
2
−1 2 1/2
(〈ψψ[x ]〉−〈ψ〉〈ψ[x ]〉)/(〈ψ 〉−〈ψ〉 ) (〈ψ[x ] 〉−〈ψ[x ]〉 ) −2 −2 0 0 0 0 ° ° ° ° ={−1.40 ,−0.69 }0,.5,.75,...1−2−7 for x ={−1.40 , x0.60 }, ±contours 0 0 −2 0 2 −2 1 1
°°
°°
xx ={ ={ 0.22 0.22 ,−0.69 , 0.60 }} 00 0 2
−2 −2 1
0
0
0
−1
−1
−1
−2 −2 1
−2 x0={−1.40°,−0.69°} 0 2
−2 1
−2 1
0
0
−1
−1
−1
−2 1
°
−2
°
x0={−1.40 ,−0.69 ,−1.99 } 0 2
−2 1
°
°
x0={ 0.22 ,−0.69 ,−1.99 } 0 2
−2 −2 1
0
0
0
−1
−1
−1
−2
−2 °
°
−2 x0={ 0.22°,−0.69°} 0 2
0
−2
°
xx ={ ={ 1.84 1.84°,−0.69 , 0.60°}} 00 0 2
x0={ 1.84°,−0.69°} 0 2
° x0={ 1.84°,−1.99 ,−0.69 } 0 2
−2
°
°
x0={−1.40 ,−1.99 } x0={ 0.22 ,−1.99 } x0={ 1.84°,−1.99°} −2 0 2 −2 0 2 −2 0 2 1 1 1 Figure 8: Contours of the κ = 1 correlation function ρn,1 [λ, ϕ] (34) vs λ & ϕ for uniform reference points ◦ ◦ ◦ ◦ (indicated by white crosshairs) λn1 = −1.40 , 0.22◦ , 1.84◦ (west to east) and 0 0 0 ϕn2 = −1.99 , −0.69 , 0.60 (south to north). Contour levels are 0, ±1/2, ±3/4, . . . ± 127/128. −1
°
−1
prelstat
−1
26−May−2008
4. An error-covariance matrix informed by a very large range of scales may be significantly non−2 diagonal, but still amenable to −2 −27] [cf. 16]. fast MRA-based operator calculus [6, x0={−1.40°,−1.99°} x0={ 0.22°,−1.99°} x0={ 1.84°,−1.99°} −2 0 2 −2 0 2 −2 0 2 Acknowledgements 1 1 1 We The MATLABr package MEXNC 0 thank Zhihong Li for useful discussions. 0 0 [45] is gratefully acknowledged. The first author thanks Mark A. Taylor for obtaining references [9, 10]. −1
prelstat
A
−1
−1
Preliminary data
26−May−2008
−2 −2 −2 For preliminary experiments we used a similar data set as described in §2, except with N1 N2 N3 = −2 0 2 −2 0 2 288−2× 360 × 38 0points in 2
{λn1 , ϕn2 } ∈ X = [λ1 , λN1 ] × [ϕ1 , ϕN2 ] = [−5.682◦ , 4.650◦ ] × [−4.452◦ , 8.472◦ ] and
[z1 , zN3 ] = [10m, 36.08km]
at a single time t = 2007/7/20 10:00 from UM v6.4. To satisfy (3) we replace (2) by a simple projection
prelstat
D {uD l , vl }
where
Πn¯N,n
{ul − (Π N1 ⊗ IN2 )ul , v l − (IN1 ⊗ Π N2 )v l }, 2 π n n ¯ −1 ¯ −1 := δn,1 cos2 π2 N −1 + δn,N sin 2 N −1 =
26−May−2008
(37) (38)
n ¯ −1 n ¯ −1 is an interpolating projection. In (38) the implicit sequence cos π N −1 could be replaced by 1 − 2 N −1 or any sequence interpolating {1, N }, {1, −1}; the present choice affects every n, but only affects the
14
1 −2
0
2
1 −2
0
2
0
0 〈ψψ[x ]〉−〈ψ〉〈ψ[x ]〉 for
−1
0
0
x ={−1.40°, 0.60°}
−1
0
1 −2
0
2
0 °
−1
°
x0={ 0.22 , 0.60 }
°
°
x0={ 1.84 , 0.60 }
−1 −2 contours −5.00e−01:2 :3.50e+00
2 −2 contours 0.00e+00:2 :1.60e+01
−2 contours −2.00e+00:2−1:5.00e−01
° ,−0.69°} 1−2 x0={−1.40 0 2 0 contours 0.00e+00:2 :5.00e+00 0 1
° ,−0.69°} 1−2 x0={ 0.22 0 2 2 contours 0.00e+00:2 :2.00e+01 0 1
° ,−0.69°} 1−2 x0={ 1.84 0 2 2 contours 0.00e+00:2 :2.00e+01 0 1
−1 0
−1 0
−1 0
−2 −1 x0={−1.40°,−0.69°} −2 0 2 −2 contours 0.00e+00:20:5.00e+00
−2 −1 x0={ 0.22°,−0.69°} −2 0 2 −2 contours 0.00e+00:22:2.00e+01
−2 −1 x0={ 1.84°,−0.69°} −2 0 2 −2 contours 0.00e+00:22:2.00e+01
1 −2
0
2
0
1 −2
0
2
0
−1
−1
x0={−1.40°,−0.69°}
1 −2
0
2
0 −1
x0={ 0.22°,−0.69°}
x0={ 1.84°,−0.69°}
−2 contours 0.00e+00:20:5.00e+00
−2 contours 0.00e+00:22:2.00e+01
−2 contours 0.00e+00:22:2.00e+01
° ,−1.99°} 1−2 x0={−1.40 0 2 contours −3.00e+00:2−1:1.00e+00 0 1
° ,−1.99°} 1−2 x0={ 0.22 0 2 contours 0.00e+00:22:1.60e+01 0 1
° ,−1.99°} 1−2 x0={ 1.84 0 2 contours −2.00e+00:21:8.00e+00 0 1
−1 0
−1 0
−1 0
−2 −2 −1 −1 ° ° x ={−1.40 ,−1.99 } x ={ 0.22°,−1.99°} 0 0 −2 0 2 −2 0 2 −2 contours −3.00e+00:2−1:1.00e+00 −2 contours 0.00e+00:22:1.60e+01 1 −2
0
2
1 −2
0
2
−2 −1 x ={ 1.84°,−1.99°} 0 −2 0 2 −2 contours −2.00e+00:21:8.00e+00 1 −2
0
2
0 Figure 9: As in Fig. 8 except samples0(33) of the covariance matrix. The0lower and upper contour limits and interval are indicated for each plot by “l:i:u”, and the color map is fixed. −1 −1 −1 x ={−1.40°,−1.99°} x ={ 0.22°,−1.99°} x ={ 1.84°,−1.99°} 0 0 0 −1 0 and Ξ because 2}, 1 prelstat 26−May−2008 fundamental wavevectors ∀m ∈ {1, . . . N −2 −2 −2 contours −3.00e+00:2 :1.00e+00 contours 0.00e+00:2 :1.60e+01 contours −2.00e+00:2 :8.00e+00 n o 1−2 1−2m + (1 + (−1) 1−2 Π N {cos π (m−1)(·−1) Π N {sin π (m−1)(·−1) 0 N −1 } = 2 12 1 − (−1) 0 N −1 } =2 0. 0 m ) cos π2N·−1 −1 , 0
B
0
0
Preliminary results
−1 −1 −1 Figure 12 shows surface (z = z1 ) plots of un,1 and vn,1 , Figs. 13 & 14 show their projections (37) to prelstat satisfy (3), and Fig. 15 shows the potential χn,1 and streamfunction ψ−2 single-precision n,1 . Due to the26−May−2008 −2 −2 data, (3) is verified only to 9 digits or so. Figure 16 shows a ψ l=12field and its−2reconstruction −2 0 0 matrix of2 1s for 0 WMψ˜l=1 2 (where M−2is a diagonal the 16 ψ˜i,1 of largest magnitude and 0s for the other ψ˜i,1 ) for all 25 bases, as a function of λn1 and ϕn2 . The error fields are shown in Fig. 17.
References [1]
R.N. Bannister, 2007: “Can wavelets improve the representation of forecast error covariances in variational data assimilation?”, Mon. Wea. Rev. 135, 387–408.
[2]
R.N. Bannister, The waveband summation transform Research Project.
[3]
G.M. Baxter, S.L. Dance, A.S. Lawless & N.K. Nichols, 2006: “Scale analysis of observational information in data assimilation”, Numerical Analysis Report 11/2006, Reading.
[4]
G.M. Baxter, S.L. Dance, A.S. Lawless, N.K. Nichols & S.P. Ballard, 2007: “Towards data assimilation for high resolution nested models”, Numerical Analysis Report 2/2007, Reading.
[5]
G. Beylkin, R. Coifman & V. Rokhlin, 1991: “Fast wavelet transforms and numerical algorithms I”, Comm. Pure Appl. Math. 44, 141–183.
[6]
G. Beylkin & M.J. Mohlenkamp, 2002: “Numerical operator calculus in higher dimensions”, Proc. Nat. Acad. Sci. 99, 10246–10251.
prelstat
26−May−2008
15
1 −2
0
2
1 −2
0
2
0
0 〈ψψ[x ]〉−〈ψ〉〈ψ[x ]〉 for
−1
0
0
x ={−1.40°, 0.60°}
−1
0
1 −2
0
2
0 °
−1
°
x0={ 0.22 , 0.60 }
°
°
x0={ 1.84 , 0.60 }
−1 −2 contours 0.00e+00:2 :3.50e+00
2 −2 contours 0.00e+00:2 :1.60e+01
−2 contours −2.00e+00:2−1:5.00e−01
° ,−0.69°} 1−2 x0={−1.40 0 2 0 contours 0.00e+00:2 :5.00e+00 0 1
° ,−0.69°} 1−2 x0={ 0.22 0 2 2 contours 0.00e+00:2 :2.00e+01 0 1
° ,−0.69°} 1−2 x0={ 1.84 0 2 2 contours 0.00e+00:2 :2.00e+01 0 1
−1 0
−1 0
−1 0
−2 −1 x0={−1.40°,−0.69°} −2 0 2 −2 contours 0.00e+00:20:5.00e+00
−2 −1 x0={ 0.22°,−0.69°} −2 0 2 −2 contours 0.00e+00:22:2.00e+01
−2 −1 x0={ 1.84°,−0.69°} −2 0 2 −2 contours 0.00e+00:22:2.00e+01
1 −2
0
2
0
1 −2
0
2
0
−1
−1
x0={−1.40°,−0.69°}
1 −2
0
2
0 −1
x0={ 0.22°,−0.69°}
x0={ 1.84°,−0.69°}
−2 contours 0.00e+00:20:5.00e+00
−2 contours 0.00e+00:22:2.00e+01
−2 contours 0.00e+00:22:2.00e+01
° ,−1.99°} 1−2 x0={−1.40 0 2 contours −2.50e+00:2−1:1.00e+00 0 1
° ,−1.99°} 1−2 x0={ 0.22 0 2 contours 0.00e+00:22:1.60e+01 0 1
° ,−1.99°} 1−2 x0={ 1.84 0 2 contours −2.00e+00:21:6.00e+00 0 1
−1 0
−1 0
−1 0
−2 −2 −1 −1 ° ° x ={−1.40 ,−1.99 } x ={ 0.22°,−1.99°} 0 0 −2 0 2 −2 0 2 −2 contours −2.50e+00:2−1:1.00e+00 −2 contours 0.00e+00:22:1.60e+01 1 −2
0
2
1 −2
0
2
−2 −1 x ={ 1.84°,−1.99°} 0 −2 0 2 −2 contours −2.00e+00:21:6.00e+00 1 −2
0
2
0 0 substituted in (33). For each λ0 , the same contour interval is Figure 10: As in Fig. 9 except with (35) n used as in Fig. 9. The error 0.1% < εM < 1.7% (36) and memory savings 98.2% < £M < 98.5% (36) −1 all λ . −1 −1 over x n={−1.40°,−1.99°} x ={ 0.22°,−1.99°} x ={ 1.84°,−1.99°} 0 0 0 −1 covaplot −2 contours −2.50e+00:2 :1.00e+00 [7]
2 −2 contours 0.00e+00:2 :1.60e+01
1 21−Aug−2008 −2 contours −2.00e+00:2 :6.00e+00
G. Beylkin & M.J. Mohlenkamp, 2005 “Algorithms for numerical analysis in high dimensions”, SIAM J. Sci. Comp. 26, 2133–2159.
1−2 1−2 0 A.C. Lorenc2& M. Thurlow,11999: 0 2 −2 “VAR scientific 0 documentation 2 paper 13: Control J. Bray, N.B. Ingleby, variable transforms”, version 1999/11/16.
[8] [9]
0 C.M. Brislawn, 1995:
[10]
“Preservation of subband symmetry in multirate signal coding”, IEEE Trans. Signal Process. 43, 3046–3050.
0
0
C.M. Brislawn, 1996: “Classification of nonexpansive symmetric extension transforms for multirate filter banks”, Appl. Comp. Harm. Anal. 3, 337–357.
−1 A. Cohen, I. Daubechies & J.-C. Feauveau, 1992: −1“Biorthogonal bases of compactly supported wavelets”, −1 Comm. Pure Appl. Math. 45, 485–560.
[11]
[12]
Concerted Research Action 731 – Propagation of uncertainty in advanced meteo-hydrological forecast systems, ongoing from 2005/6/27.
covaplot −2 [13]
[14]
−2
−2
21−Aug−2008
M.J.P. Cullen, 2003: “Four-dimensional variational data assimilation: A new formulation of the background-error covariance matrix based on a potential-vorticity representation”, Q.J.R. Meteorol. Soc. 129, 2777–2796.
−2 −2 forecasting”, Physica 0 D 196, 1–27. 2 −2area data assimilation 0 2 S.L. Dance, 2004:0“Issues in high2resolution limited for quantitative precipitation
[15]
I. Davison, 2006: “Scale analysis of short term forecast errors”, MSc thesis in Mathematics, Reading.
[16]
A. Deckmyn & L. Berre, 2004: “Using wavelets for the modelling of LAM background error covariances” [36].
[17]
G. Deslauriers & S. Dubuc, 1989: “Symmetric iterative interpolation processes”, Constr. Approx. 5 49–68.
[18]
D. Donoho, 1993: “Smooth wavelet decompositions with blocky coefficient kernels”, Recent Advances in Wavelet Analysis, Academic Press, pp. 259–308.
covaplot
21−Aug−2008
[19]
D. Donoho, A. Maleki, M. Shahram, J. Buckheit, M. Clerc, J. Kalifa, S. Mallat, T. Yu: “WaveLab . . . library of MATLAB routines for wavelet analysis”, (http://www- stat.stanford.edu/~wavelab/).
[20]
A. Fournier, 2000: “Introduction to orthonormal wavelet analysis with shift invariance: Application to observed atmospheric blocking spatial structure”, J. Atmos. Sci. 57, 3856–3880 (http://dx.doi.org/10.1175/1520- 0469(2000)057< 3856:ITOWAW> 2.0.CO;2).
[21]
A. Fournier, 2002: “Atmospheric energetics in the wavelet domain. Part I: Governing equations and interpretation for idealized flows”, J. Atmos. Sci. 59, 1182–1197 (http://dx.doi.org/10.1175/1520- 0469(2002)059< 1182:AEITWD> 2.0.CO;2).
[22]
A. Fournier, 2003: “Atmospheric energetics in the wavelet domain II: Time-averaged observed atmospheric blocking”, J. Atmos. Sci. 60, 319–338 (http://dx.doi.org/10.1175/1520- 0469(2003)060< 0319:AEITWD> 2.0.CO;2).
16
−1
0
0
−1
x ={−1.40°, 0.60°} 0
−2 −2 contours −7.50e−01:2 :7.50e−01
1 −2
0
2
−1
°
〈ψψ[x ]〉−〈ψ〉〈ψ[x ]〉 for 0
0
x ={−1.40°, 0.60°}
1 −2
−1
0
−2 −2 contours −7.50e−01:2 :7.50e−01
°
°
x0={ 1.84 , 0.60 }
−1 −2 contours −1.50e+00:2 :1.00e+00 −2 contours −2.50e−01:2−3:2.50e−01
0
2
0
0 −1
°
x0={ 0.22 , 0.60 }
1 −2
0
2
0 °
−1
°
x0={ 0.22 , 0.60 }
°
°
x0={ 1.84 , 0.60 }
−1 −2 contours −1.50e+00:2 :1.00e+00 −2 contours −2.50e−01:2−3:2.50e−01
° ,−0.69°} 1−2 x0={−1.40 0 2 contours −7.50e−01:2−2:1.25e+00 0 1
° ,−0.69°} 1−2 x0={ 0.22 0 2 contours −1.50e+00:2−1:1.50e+00 0 1
° ,−0.69°} 1−2 x0={ 1.84 0 2 contours −2.00e+00:2−1:2.00e+00 0 1
−1 0
−1 0
−1 0
−2 −2 −2 −1 −1 −1 ° ° ° ° x0={−1.40 ,−0.69 } x0={ 0.22 ,−0.69 } x0={ 1.84°,−0.69°} −2 0 2 −2 0 2 −2 0 2 −2 contours −7.50e−01:2−2:1.25e+00 −2 contours −1.50e+00:2−1:1.50e+00 −2 contours −2.00e+00:2−1:2.00e+00 1 −2
0
2
0 −1
1 −2
0
2
0 −1
x0={−1.40°,−0.69°}
1 −2
0
2
0 −1
x0={ 0.22°,−0.69°}
x0={ 1.84°,−0.69°}
−2 contours −7.50e−01:2−2:1.25e+00 −2 contours −1.50e+00:2−1:1.50e+00 −2 contours −2.00e+00:2−1:2.00e+00 ° ,−1.99°} 1−2 x0={−1.40 0 2 contours −7.50e−01:2−2:7.50e−01 0 1
° ,−1.99°} 1−2 x0={ 0.22 0 2 contours −2.00e+00:2−1:1.50e+00 0 1
° ,−1.99°} 1−2 x0={ 1.84 0 2 contours −1.50e+00:2−1:1.00e+00 0 1
−1 0
−1 0
−1 0
−2 −2 −2 −1 −1 −1 ° ° ° ° x ={−1.40 ,−1.99 } x ={ 0.22 ,−1.99 } x ={ 1.84°,−1.99°} 0 0 0 −2 0 2 −2 0 2 −2 0 2 −2 contours −7.50e−01:2−2:7.50e−01 −2 contours −2.00e+00:2−1:1.50e+00 −2 contours −1.50e+00:2−1:1.00e+00 1 −2
0
2
1 −2
0
2
1 −2
0
2
0Figure 11: Difference between Figs.0 9 & 10. Much smaller contour intervals 0 are used, as indicated. −1
−1
x0={−1.40°,−1.99°}
covaplot −2 contours −7.50e−01:2−2:7.50e−01
1−2
0
2
−1
x0={ 0.22°,−1.99°}
−1 21−Aug−2008 −2 contours −2.00e+00:2−1:1.50e+00 −2 contours −1.50e+00:2 :1.00e+00
1−2
0
2
1−2
0
0
0
−1
−1
−1
covaplot −2
−2
−2
−2
covaplot
0
2
x0={ 1.84°,−1.99°}
−2
0
2
−2
0
2
21−Aug−2008
0
2
Figure 12: Surface wind components un,1 (left) & vn,1 (right) vs λn1 & ϕ21−Aug−2008 n2 .
17
D Figure 13: As in Fig. 12 but un,1 − uD n,1 & vn,1 − vn,1 .
D Figure 14: As in Fig. 12 but uD n,1 & vn,1 (eq. 37, impermeable boundary as black line).
[23]
A. Fournier, 2005: “Instantaneous wavelet energetic transfers between atmospheric blocking and local eddies”, Journal of Climate 18, 2151–2171 (http://dx.doi.org/10.1175/JCLI3381.1).
[24]
A. Fournier, 2006: “Exact calculation of Fourier series in nonconforming spectral-element methods”, J. Comp. Phys. 215, 1–5 (http://dx.doi.org/10. 1016/j.jcp.2005.11.023).
[25]
A. Fournier, 2006: “Multiresolution analysis and simulation using spectral elements”, NCAR Computational and Information Systems Laboratory Annual Report (http://www.cisl.ucar.edu/nar/2006/catalog/Fournier1.pdf).
[26]
A. Fournier, G. Beylkin & V. Cheruvu, 2005: “Multiresolution adaptive space refinement in geophysical fluid dynamics simulation”, Lecture Notes Comp. Sci. Eng. 41, 161–170.
[27]
A. Fournier, D. Rosenberg & A. Pouquet, 2005: “Dynamically adaptive computation of multiscale coherent-structure interactions using GASpAR”, International Conference on High Reynolds Number Vortex Interactions, Toulouse (http://www-vortex.mcs.st- and.ac.uk/~jean/toulouse/abstracts/Fournier/ Fournier_abs.pdf).
[28]
A. Fournier, M.A. Taylor & J.J. Tribbia, 2004: “The spectral element atmosphere model (SEAM): High-resolution parallel computation and localized resolution of regional dynamics”, Monthly Weather Review 132, 726–748 (http://dx.doi.org/10.1175/1520- 0493(2004)132< 0726:TSEAMS> 2.0.CO;2).
[29]
G. Gaspari & S.E. Cohn, 1999: “Construction of correlation functions in two and three dimensions”, Q.J.R. Meteorol. Soc. 125, 723–757.
[30]
T. Gneiting, 1999: “Correlation functions for atmospheric data analysis”, Q.J.R. Meteorol. Soc. 125, 2449–2464.
[31]
V. Guidard & C. Fischer, 2004: “About the use of a global analysis in a LAM assimilation”, [36].
[32]
A.K. Gupta & D.K. Nagar, 1999: Matrix Variate Distributions, CRC Press, 367 pp.
[33]
D.T. Hristopulos, 2002: “New anisotropic covariance models and estimation of anisotropic parameters based on the covariance tensor identity”, Stochastic Environmental Research and Risk Assessment 16, 43–62.
Figure 15: As in Fig. 12 but χn,1 & ψn,1 (7). 18
econstruction of ψ from the 16 ( 0.02%) largest wavelets 8 6 4 2 0 −2 −4 Deslauriers−1 8 6 4 2 0 −2 −4
Deslauriers−3 8 6 4 2 0 −2 −4
Ave−Interp−4 8 6 4 2 0 −2 −4
Ave−Interp−6 8 6 4 2 0 −2 −4
CDF−2−2 8 6 4 2 0 −2 −4
Villasenor−1
356 360 364
CDF−3−9
Villasenor−3
Villasenor−4
Villasenor−5 8 6 4 2 0 −2 −4
356 360 364
Figure 16: Contours of ψ and its reconstructions.
19
CDF−4−4 8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4 356 360 364
CDF−3−1 8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4 356 360 364
CDF−2−8
CDF−3−7
Villasenor−2
CDF−1−5 8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
CDF−1−3
CDF−2−6
CDF−3−5
Ave−Interp−2 8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
CDF−1−1
CDF−2−4
CDF−3−3
Deslauriers−7 8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
8 6 4 2 0 −2 −4
Deslauriers−5 8 6 4 2 0 −2 −4
356 360 364
onstruction error of ψ from the 16 ( 0.02%) largest wavelets 8 6 4 2 0 −2 −4 Deslauriers−1
Deslauriers−3
Deslauriers−5
Deslauriers−7
Ave−Interp−2
8
8
8
8
8
6
6
6
6
6
4
4
4
4
4
2
2
2
2
2
0
0
0
0
0
−2 −4
−2 −4
−2 −4
−2 −4
−2 −4
Ave−Interp−4
Ave−Interp−6
CDF−1−1
CDF−1−3
CDF−1−5
8
8
8
8
8
6
6
6
6
6
4
4
4
4
4
2 0
2 0
2 0
2 0
2 0
−2 −4
−2 −4
−2 −4
−2 −4
−2 −4
CDF−2−2
CDF−2−4
CDF−2−6
CDF−2−8
CDF−3−1
8
8
8
8
8
6 4
6 4
6 4
6 4
6 4
2
2
2
2
2
0 −2
0 −2
0 −2
0 −2
0 −2
−4
−4
−4
−4
−4
CDF−3−3
CDF−3−5
CDF−3−7
CDF−3−9
CDF−4−4
8
8
8
8
8
6 4
6 4
6 4
6 4
6 4
2 0
2 0
2 0
2 0
2 0
−2 −4
−2 −4
−2 −4
−2 −4
−2 −4
Villasenor−1
Villasenor−2
Villasenor−3
Villasenor−4
Villasenor−5
8
8
8
8
8
6 4
6 4
6 4
6 4
6 4
2 0
2 0
2 0
2 0
2 0
−2 −4
−2 −4
−2 −4
−2 −4
−2 −4
356
360
364
356
360
364
356
360
364
356
360
364
356
Figure 17: Reconstruction error as in Fig. 16 but with half the contour interval.
20
360
364
[34]
N.B. Ingleby, 2001: “The statistical structure of forecast errors and its representation in The Met. Office Global 3-D Variational Data Assimilation Scheme”, Q.J.R. Meteorol. Soc. 127, 209–231.
[35]
I.M. Johnstone & B.W. Silverman, 1997: “Wavelet threshold estimators for data with correlated noise”, J. Roy. Stat. Soc. B 59, 319–351.
[36]
Joint SRNWP/Met Office/HIRLAM Workshop on Variational Assimilation: Towards 1-4km Resolution, 2004/11.
[37]
Z. Li, S.P. Ballard & A.C. Lorenc, 2004: “Development of scale-selective data assimilation for a 4km version of the Unified Model” [36].
[38]
A.C. Lorenc, 2003: “Modelling of error covariances by 4D-Var data assimilation”, Q.J.R. Meteorol. Soc. 129, 3167–3182.
[39]
A.C. Lorenc, 2007: “Ideas for adding flow-dependence to the Met Office VAR system” [56].
[40]
A.C. Lorenc & T. Payne, 2007: “4D-Var and the butterfly effect: Statistical four-dimensional data assimilation for a wide range of scales”, Q.J.R. Meteorol. Soc. 133, 607–614.
[41]
A.C. Lorenc & F. Rawlins, 2005: “Why does 4D-Var beat 3D-Var?”, Q.J.R. Meteorol. Soc. 131, 3247–3257.
[42]
A.C. Lorenc, S.P. Ballard, R.S. Bell, N.B. Ingleby, P.L.F. Andrews, D.M. Barker, J.R. Bray, A.M. Clayton, T. Dalby, D. Li, T.J. Payne & F.W. Saunders, 2000: “The Met Office global three-dimensional variational data assimilation scheme”, Q.J.R. Meteorol. Soc. 126, 2991–3012.
[43]
B. Macpherson, P. Andrews, S. Harcourt, B. Ingleby, B. Chalcraft, A. Maycock, R. Renshaw, C. Parrett, S. Anderson, M. Sharpe, D. Harrison & M. Gibson, 2001: “The operational mesoscale data assimilation system 1999-2001: Implementation of 3DVar and later upgrades”, Forecasting Research Technical Report 374.
[44]
Mesoscale Variational Data Assimilation Workshop, 2000/5.
[45]
“MEXNC . . . mex-file interface to NetCDF files for MATLAB . . . roughly a one-to-one equivalence with the C API for NetCDF” http://mexcdf. sourceforge.net/
[46]
R.J. Purser, 2004: “A geometrical approach to the synthesis of smooth anisotropic covariance operators for adaptive meteorological data assimilation” [36].
[47]
R.J. Purser, W.-S. Wu, D.F. Parrish & N.M. Roberts, 2003: “Numerical aspects of the application of recursive filters to variational statistical analysis. Part II: Spatially inhomogeneous and anisotropic general covariances”, Mon. Wea. Rev. 131, 1536–1548.
[48]
F. Rawlins, S.P. Ballard, K.J. Bovis, A.M. Clayton, D. Li, G.W. Inverarity, A.C. Lorenc & T. Payne “The Met Office global four-dimensional variational data assimilation scheme”, Q.J.R. Meteorol. Soc. 133, 347–362.
[49]
A. Rhodin & H. Anlauf, 2007: “Representation of inhomogeneous, non-separable covariances by sparse wavelet-transformed matrices” [56].
[50]
D. Rosenberg, A. Fournier, P. Fischer & Annick Pouquet, 2006: “Geophysical-astrophysical spectral-element adaptive refinement (GASpAR): Object-oriented h-adaptive fluid dynamics simulation”, J. Comp. Phys. 215, 59–80 (http://dx.doi.org/10.1016/j.jcp.2005.10.031).
[51]
D.B. Rowe, 2003: Multivariate Bayesian Statistics CRC Press, 329 pp.
[52]
A. Staniforth, A. White, N. Wood, J. Thuburn, M. Zerroukat, E. Cordero, T. Davies & M. Diamantakis, 2006: “Joy of U.M. 6.3 — Model formulation”, Unified Model Documentation Paper 15 (http://www.metoffice.gov.uk/research/nwp/publications/papers/unified_model/).
[53]
J.D. Villasenor, B. Belzer & J. Liao, 1995: “Wavelet filter evaluation for image compression”, IEEE Trans. Image Proc. 4, 1053–1060.
[54]
M.V. Wickerhauser, 1994: Adapted wavelet analysis from theory to software, A.K. Peters, 486 pp.
[55]
M.A. Wlasak, S.P. Ballard & M.J. Cullen, 2004: “Limited area 4DVar over a North Atlantic European domain” [36].
[56]
Workshop on flow-dependent aspects of data assimilation, 2007/6.
21
