In the wireless communication system, the pilot is essential for channel estimation. This paper study the pilot design of the 1-bit massive multi-input multi-output system with the flat Rayleigh fading channel. This paper aims to design the pilot matrix minimizing the estimation error of the Bussgang linear minimum mean squared error (BLMMSE) channel estimator, of which the estimation error is given by a closed-form formula. An approximate BLMMSE channel estimation error is derived in the low signal-to-noise ratio (SNR) region. The lower bound of the approximate channel estimation error for all possible pilot matrices is derived. The analysis shows that the pilot matrix with pairwise orthogonal and equal-norm row vectors and equal-norm column vectors is optimal in terms of approximate channel estimation error. The hypersphere-based steepest descent (HSD) algorithm is proposed to design the pilot matrix by optimizing the performance of the BLMMSE estimator. The numerical results show that the proposed HSD algorithm has a noticeable improvement in channel estimation compared with the DFT pilot and the random pilot. The HSD algorithm is also valid for the near maximum likelihood estimator and SNR mismatch.