An accurate and effective seismic wavelet estimation technique has extreme significance in the seismic data processing for analyzing the earth’s subsurface layer information. The seismic wavelet to be determined is modeled as a moving average (MA) process and assumed to be driven by a zero mean, non-Gaussian, statistically independent, and identically distributed (IID) process. In order to estimate the MA model parameter from the observed noisy seismic signal, we pose this as a blind system identification (BSI) problem. In the BSI, a set of multivariate polynomial equations is obtained by matching higher order cumulant of observed noisy data with a higher order moment of blind system’s impulse response. The Groebner bases that form the solution to this set of equations are obtained using the proposed algorithm. Numerical results demonstrate that the proposed method has a lower estimation error as compared to the previously reported methods.