Radiometric normalization is a fundamental and important preprocessing method for remote sensing applications using multitemporal satellite images due to uncertainties of at-sensor radiances caused by different sun angles and atmospheric conditions. In case the atmospheric model and ground measurements are unavailable during data acquisitions, relative normalization is an alternative method which minimizes the radiometric differences among images without the requirement of additional information. The keys to a successful relative normalization are the selection of pseudo invariant features (PIFs) from bitemporal images and the regression of selected PIFs for transformation coefficient determination. Previous studies on transformation coefficient determination adopted band-by-band regression. These studies have obtained satisfactory normalization results; however, they have not fully considered the spectral inconsistency problem caused by individual band regression. To alleviate this problem, this study proposed a constrained orthogonal regression, which enforces pixel spectral signatures to be as consistent as possible during radiometric normalization while band regression quality is preserved. In addition, instead of selecting one of the input images as reference for radiometric transformation, a common radiometric level located between bitemporal images is selected as the reference to further reduce possible spectral inconsistency. Qualitative and quantitative analyses of several bitemporal images acquired by the Landsat 8 sensor were conducted to evaluate the proposed method with the measurements of spectral distance and similarity. The experimental results demonstrate the superiority of the proposed method to related regression and radiometric normalization methods, in terms of spectral signature consistency.