The viscoelasticity of the subsurface media varies spatially, and such viscoelasticity can be represented concisely by a wave equation in the form of fractional temporal derivative (FTD). We have developed a strategy for simulating seismic waves propagating through a heterogeneous viscoelastic model. The FTD is transferred to fractional spatial derivatives (FSDs), and the FSDs are implemented through the fast Fourier transform (FFT), for improving the computational efficiency. However, the FFT implementation is not rigorously applicable to the heterogeneous model. In this paper, we have reformulated the FSD wave equation by introducing a spatial-position dependent filter. This spatial filter corrects the errors that are caused by the assumption of non-heterogeneity in the FFT implementation. This formulation appropriately represents the viscoelastic effect in seismic wave propagation, leading to the improvement on the accuracy of numerical simulation.