This paper studies a statistical model for heteroscedastic (i.e., power fluctuating) signals embedded in white Gaussian noise. Using the Riemannian geometry theory, we propose an unified approach to tackle several problems related to this model. The first axis of contribution concerns parameters (signal subspace and power factors) estimation, for which we derive intrinsic Cramér-Rao bounds and propose a flexible Riemannian optimization algorithmic framework in order to compute the maximum likelihood estimator (as well as other cost functions involving the parameters). Interestingly, the obtained bounds are in closed forms and interpretable in terms of problem's dimensions and SNR. The second axis of contribution concerns the problem of clustering data assuming a mixture of heteroscedastic signals model, for which we generalize the Euclidean K-means++ to the considered Riemannian parameter space. We propose an application of the resulting clustering algorithm on the Indian Pines segmentation problem benchmark.