Multiple removal is one of the important preprocessing steps in a seismic data processing sequence. For an accurate multiple prediction, a full 3D method is required. However, the application of such methods is often limited by practical and economic constraints. Therefore, in practical situations, 2D prediction methods are used, such as 2D surface-related multiple elimination. However, the resulting predicted multiples may have significant temporal shift, spatial mismatch, and amplitude inconsistency, compared with true 3D multiples. Adaptive multiple subtraction based on 2D blind separation of convolved mixtures (BSCM) has been proposed to estimate the 2D matching filter in a single gather. To improve the flexibility of the adaptive multiple subtraction for the inconsistencies between the 2D predicted multiples and true 3D multiples, we evaluated the adaptive multiple subtraction as a problem of 3D BSCM. In the proposed method, the predicted multiples were modeled as the convolution of the true multiples with a 3D kernel, whose third dimension is in the gather direction. By maximizing the non-Gaussianity of the estimated primaries, the iterative reweighted least-squares algorithm was exploited to obtain the 3D matching filter, which is the inverse of the 3D kernel. To avoid the possible overfitting problem introduced by the 3D matching filter, the proposed method fit several seismic gathers using one 3D matching filter. In addition, by using the non-Gaussian maximization criterion, the proposed method alleviated the orthogonality assumption used by the least-squares subtraction method. Furthermore, the proposed method can eliminate the temporal and spatial mismatches between the 2D predicted multiples and true 3D multiples better than the 2D BSCM subtraction method. Tests on synthetic and field data sets demonstrated the effectiveness of the 3D BSCM subtraction method.