We report the orbital X-ray variability of the high-mass X-ray binary (HMXB) GX 301-2. GX 301-2 underwent a spin-up process in 2018-2020 with the period evolving from ~685-670 s. The energy-resolved pulse profiles of the pulsar at 1-60 keV varied from single-peaked and sinusoidal shapes to multipeaked ones across different orbital phases. Pulse fractions evolving over the orbit had negative correlations with the X-ray flux. The broad-band X-ray energy spectrum of the pulsar can be described with a partially covered negative-positive cut-off power-law continuum model. Near the periastron passage of the pulsar we found strong variation in the additional column density ($N_{\mathrm{H}_{2}}$), which correlated with variation of the flux. Curves of growth for both Fe K`and Fe Klines were plotted to investigate the distribution of matter around the neutron star. We also found evidence for two cyclotron absorption lines in the phase-averaged spectra in GX 301-2, with one line of 30-42 keV and the other line varying over 48-56 keV. The centroid energies of both lines show a similar relationship with X-ray luminosity: positive correlation in the lower luminosity range, and a negative relation above a critical luminosity of $10^{37}\, \rm erg\, s^{-1}$. We estimate the surface magnetic field of the neutron star in GX 301-2 to be ~(0.5-2) × 10¹³ G. The two cyclotron line energies have a nearly fixed ratio of ~1.63 while having a low strength ratio (~0.05), suggesting that these two features may actually be one line.