I want to calculate the velocity gradient using a finite difference scheme. The geometry looks like this:

1 2 3
(wall) o..|..o....o (fluid)

Now I am interested in the velocity gradient at the last lattice node (number 2) in front of the wall (fullway bounceback boundary condition). Is is too naive to approximate it by

du/dz = [u_2 / (h / 2) + (u_3 - u_2) / h] / 2

?
h is the lattice constant. So basically, I calculate the gradient between nodes 2 and 3 and the gradient between the wall and node 2 and then just use the mean value.
This approach is only used for post-processing, i. e. comparing with another data set. It is not used within the LB code.

I could imagine that what you want is an estimate of du/dz(z=2) which is second-order accurate with respect to the grid spacing, and thus compatible with the accuracy of the LB scheme. You can obtain this by Taylor-expanding the velocity around the spot of interest (z=2), using

this is a very nice answer to my question. Thanks a lot! I will immediately implement this approach.
By the way, can you recommend a review of finite difference methods used in LBM?