Seasonal fluctuations of glacier speed largely reflect the change in subglacial hydrology and bed conditions throughout the year. Satellite observations provide rich optical and Synthetic Aperture Radar (SAR) images with a few-day revisit duration, which enables us to derive scene-pair velocity maps. However, the number and the quality of the scene-pair velocity maps often vary in different seasons, resulting in a temporally inhomogeneous data set that is challenging for glacier modeling. Here we present a novel optimization method to extract glacier velocity time series with a uniform and short sampling interval from the scene-pair maps. If the scene-pair velocity maps stem from Sentinel-1 and Sentinel-2 images, the sampling interval can be six days or even less, which can capture the seasonal ice-speed variation well. Our regularized optimization is based on how the glacier speed is mathematically averaged over time. This method extracts the time series on a pixel-by-pixel basis, and hence, a data cube of ice speed (i.e., ice speed of a glacier based on two horizontal coordinates and one temporal coordinates, all with uniform sampling intervals) can be derived. It considers two sources for assessing the uncertainty of the final results: (1) the inherited uncertainty from each scene-pair map and (2) the uncertainty of ice-speed variation at the period no scene-pair maps cover. With two manually tuned hyperparameters, this method can give us realistic error assessments. We validate this method by targeting several glaciers in High Mountain Asia, where external field or weather data exist. Our results show that this regularized optimization is capable of producing the glacier verlocity data cube that captures the seasonal variations and sub-seasonal events (e.g., evolution of a glacier surge) during the observation period (2017-2024) using roughly a few hundred optical and SAR images.