The invention discloses a seismic data multiple elimination method. The method comprises the following steps: S1, inputting three-dimensional data d(offx,t), and according to three-dimensional space coordinate information, extracting channel sets as three-dimensional surface element channel sets dbin(y,x,t); S2, calculating nonlinear trend fitting coefficients p(x) and p(y), writing out conversionoperators Lx, Ly and Lxy, and calculating Radon domain data M; S3, cutting the Radon domain data to obtain a multiple model mra, and performing calculation and conversion to obtain a time-space domain multiple model; S4, performing subtraction to obtain three-dimensional surface element channel sets p after multiple elimination; and S5, sequentially arranging the channel sets according to scalaroffsets to obtain three-dimensional data after multiple elimination. The high-frequency calculation is restrained through a low-frequency result; the low-frequency result does not generate an alias frequency in operation, so that a relatively good high-resolution effect still can be achieved under the condition of sampling sparseness; and meanwhile, the process of improving the resolution by a conventional iterative method is avoided, so that the calculation efficiency is greatly improved.