A sparse squared envelope is crucial for efficient and accurate diagnosis of bearing faults. Blind deconvolution (BD) is a well-established sparse feature enhancement method for the diagnostics of rolling bearings. Traditional BD methods, such as minimum entropy deconvolution (MED), are susceptible to random transients, making it difficult to enhance fault features of rolling bearings subject to strong random shocks. Deconvolution methods that take the fault characteristic frequency (or fault impulse period) of interest as an algorithm input parameter, such as maximum second-order cyclostationarity blind deconvolution (CYCBD), can alleviate this deficiency. However, bearing fault features are difficult to be enhanced by these methods when the specified characteristic frequency deviates from the actual value greatly. To overcome these problems, the modified smoothness index (MSI) of the squared envelope is proposed as the objective function of the deconvolution method, and a new BD method is developed to achieve a sparse squared envelope for fault diagnosis of rolling bearings. Furthermore, the methodology is extended to the frequency domain, and another new BD method that utilizes the MSI of the squared envelope spectrum as the objective function is established to achieve a sparse squared envelope spectrum for bearing diagnostics. These two proposed BD methods are robust to random transients and do not require characteristic frequency or impulse period as an input parameter for feature enhancement. The performance of the two proposed BD methods is verified on experimental datasets from two different railway axle bearing test rigs and compared with the state-of-the-art deconvolution methods. The results show that the two proposed methods can effectively enhance repetitive transient features in noisy vibration signals and accurately diagnose different faults of railway axle bearings.