The Richardson extrapolation is widely utilized in numerically solving the differential equations since it enjoys both high accuracy and ease of implementation. For the Riesz space fractional diffusion equation, the fractional centered difference operator is employed to approximate the fractional derivative, then the Richardson extrapolation methods for two difference schemes are constructed with the help of the asymptotic expansions of the discrete solutions. Specifically, for the Crank–Nicolson difference scheme, the extrapolation method contains two extrapolation formulae that achieve the fourth order and the sixth order both in the temporal and spatial directions, respectively. The extrapolation method for the compact difference scheme involves one extrapolation formula by which the sixth order can be obtained when the time step size is proportional to the squares of the space step size. The maximum norm error estimates of the extrapolation solutions are proved by the discrete fractional Sobolev embedding inequalities. The extension to the high dimensional and nonlinear cases is also demonstrated. Numerical results verify the theoretical convergence orders and efficiency of our methods.