There is an unmet need to accurately identify the locations of innervation zones (IZs) of spastic muscles, so as to guide botulinum toxin (BTX) injections for the best clinical outcome. A novel 3D IZ imaging (3DIZI) approach was developed by combining the bioelectrical source imaging and surface electromyogram (EMG) decomposition methods to image the 3D distribution of IZs in the target muscles. Surface IZ locations of motor units (MUs), identified from the bipolar map of their MU action potentials (MUAPs) were employed as a prior knowledge in the 3DIZI approach to improve its imaging accuracy. The performance of the 3DIZI approach was first optimized and evaluated via a series of designed computer simulations, and then validated with the intramuscular EMG data, together with simultaneously recorded 128-channel surface EMG data from the biceps of two subjects. Both simulation and experimental validation results demonstrate the high performance of the 3DIZI approach in accurately reconstructing the distributions of IZs and the dynamic propagation of internal muscle activities in the biceps from high-density surface EMG recordings.