Backus average

[mspringer64]

There seems to be an error in the Backus average. For some choices of the interval length lb the averaged logs do not follow the original logs. This seems to be a normalization problem. The resulting acoustic impedance, i.e. the product VP*RHO of the averaged quantities is correct, i.e. there is a cancellation of errors which suggest that the error is multiplicative. Also the error becomes small if the averaging interval becomes large. The attached figure illustrates the problem, where the averaging length is once equal to 6 log increments and once for 7 log increments.

[mspringer64]

The implementation of the anisotropic Backus average in the Bruges package has two bugs
which are related to the moving average implementation

1.) For even number of increments the moving average values are shifted by one increment.
This moving average is only centered at the correct value if the number of increments
is odd. This can be easily seen by taking the moving average of a function y = x.
2.) The normalization in the moving average can easily become wrong by rounding errors. If the
log increment is dz = 0.1524 then an averaging intervall of lb = 1.0668 would correspond to exactly
7 samples. However, in the moving average there is a boxcar function which like
boxcar = np.ones(int(length))/length
Here length = lb/dz. The problem is that the int(1.0668/0.1524) evaluates to 6
while 1.0668/0.1524 evaluates to 6.9999999999999...

To solve the problem one can:
1.) Restrict the average to odd numbers. For instance
n = int(dz/lb)
length = n + n%2 -1
Pass the length to the moving_average and not int(dz/lb) in backus.
Alternatively one can calculate = 0.5 * (moving_average(f, n+1) + f)
2.) Change the boxcar in the to
boxcar = np.ones(int(length))
boxcar = boxcar/len(boxcar)
This will always guarantee that boxcar is normalizesd to unity.

Here is a modified version of the moving_average
def moving_average(a, length, mode='same'):
"""
Computes the mean in a moving window using convolution. This is
the implementation from the bruges package but corrected for the
bug when the length is even. Also the normalisation error which might
occur by rounding errors is removed.

Example:
    >>> test = np.array([1,1,9,9,9,9,9,2,3,9,2,2,np.nan,1,1,1,1])
    >>> moving_average(test, 5, mode='same')
    array([ 2.2,  4. ,  5.8,  7.4,  9. ,  7.6,  6.4,  6.4,  5. ,  3.6,  nan,
            nan,  nan,  nan,  nan,  0.8,  0.6])
"""
mm   = round(length)
nn     = mm + 1 - mm%2 # make sure nn is odd
n2     = nn//2
padded = np.pad(a, n2, mode='edge')
boxcar = np.ones(nn)/nn
smooth = np.convolve(padded, boxcar, mode='same')
smooth = smooth[n2:-n2]
if mm%2 > 0:
    return smooth
else:
    return 0.5*(smooth+a)