In hyperspectral images, endmembers characterizing one class of ground object may vary due to illumination, weathering, slight variations of the materials. This phenomenon is called intraclass endmember variability which is one of the important factors affecting the performance of unmixing. However, intra-class endmember variability is often ignored in unmixing, which causes a decrease in the accuracy of unmixing. How to deal with intra-class endmember variability is the focus. To address this problem, we propose a novel hyperspectral unmixing method based on Least Squares Twin Support Vector Machines (ULSTWSVM). ULSTWSVM uses multiple training samples (endmembers) to model a pure class, which takes intra-class endmember variability into account in unmixing. At the same time, ULSTWSVM obtains abundances by calculating the distances from the mixed pixels to the classification hyperplanes, which is simple and efficient. ULSTWSVM mainly comprises three steps: (1) to obtain the two non-parallel classification hyperplanes by solving two quadratic programming problems (QPPs) in least squares sense, (2) to calculate distances from the mixed pixels to classification hyperplanes, and (3) to normalize the distances and convert them to abundances. Experimental results on both synthetic and real hyperspectral data show that the proposed method outperforms the methods used for comparison.