Hammerstein systems are composed by the cascading of a static nonlinearity and a linear system. In this paper, a methodology for identifying such systems using a combination of least squares support vector machines (LS-SVM) and best linear approximation (BLA) techniques is proposed. To do this, a novel method for estimating the intermediate variable is presented allowing a clear separation of the identification steps. First, an approximation to the linear block is obtained through the BLA of the system. Then, an approximation to the intermediate variable is obtained using the inversion of the estimated linear block and the known output. Afterwards, a nonlinear model is calculated through LS-SVM using the estimated intermediate variable and the known input. To do this, the regularisation capabilities of LS-SVM play a crucial role. Finally, a parametric re-estimation of the linear block is made. The method was tested in three examples, two of them with hard nonlinearities, and was compared with four other methods showing very good performance in all cases. The obtained results demonstrate that also in the presence of noise, the method can effectively identify Hammerstein systems. The relevance of these findings lies in the fact that it is shown how the regularisation allows to bypass the usual problems associated with the noise backpropagation when the inversion of the estimated linear block is used to compute the intermediate variable.