## 1 Introduction  
Remote sensing has proven itself to be essential in land change detection [1], crop growth monitoring [2], and forest succession monitoring [3]. Different kinds of satellite image data, including Landsat, MODIS, Sentinel, GF, and _ZY_, among others, are achieving great success. Building on this success, by integrating different types of satellite data, has the potential to solve ecological, natural resource, and human development challenges through large coverage remote sensing images, capitalizing on our ability to process large datasets [4, 5, 6]. However, among platforms, images differ due to variation in temporal, spatial, radiometric, and observation angles. Solving image data inconsistency is the most important issue that will allow the application of multiple remote sensing images. Radiometric value is essential in remote sensing, so radiometric normalization, which is a preprocessing procedure to eliminate radiometric differences among images, has been studied by several researchers with the development of many radiometric normalization algorithms [7, 8, 9, 10]. Radiometric normalization algorithms can be divided into two main categories based on the transformation of grayscale values to physical signals: absolute radiometric normalization requiring accurate sensor calibration parameters and atmospheric properties [11, 12], and relative radiometric normalization (RRN), which is an image\(-\)based approach using one image as a reference to normalize another image through its radiometric characteristics [13, 14, 15]. Due to parameter differences among sensors and the difficulty in collectingatmospheric parameters, the relative radiometric normalization method is generally used for image normalization that does not require extra parameters [16; 17]. Hence, in this process, the subject image tends to need similar radiometric conditions as the reference image. Nevertheless, images from different satellite sensors have inherent bandwidth and spectral response differences (among others) that limits the practical application of absolute radiometric normalization but facilitates relative radiometric normalization algorithms in many applications. Among these relative radiometric normalization techniques, linear regression using pseudo\(-\)invariant features (PIFs) has been promoted as an empirical and practical method to normalize radiometric distortion caused by non\(-\)land surface\(-\)related factors [18; 19]. This approach assumes that reference and subject images have a consistent response and adjusts the subject image's radiometric properties with paired pixels from the two images [20; 21]. The key to preprocessing is to extract accurate PIFs, allowing image data collected at different times or by different sensors to be normalized. Comparing seven relative radiometric normalization methods with Landsat MSS images at different times, [PERSON], [PERSON] [19] and [PERSON], [PERSON] [22] have shown that accurately extracted PIFs can improve image radiometric consistency for the same sensor type at different times. Accurate PIFs are essential to compare sensor images. [PERSON], [PERSON] [7], [PERSON], [PERSON]\(-\)[PERSON] [16], [PERSON], [PERSON] [17], and [PERSON], [PERSON] [23] normalized images from different sensor types (e.g., Landsat, Sentinel\(-\)2, Gaofen, and super\(-\)resolution images) using a PIF\(-\)based method and reported coherent image correction among varying time series and sensor types. The biggest difference between these studies is the processing of manually or automatically PIFs selection. Automatically extracting PIFs makes the radiometric correction processing more efficient and has been proposed by several researchers [24; 25; 26; 27; 28]. In these studies, PIFs are extracted automatically and sorted into a radiometric control set, radiometric control points, no\(-\)change set, and unchanged pixels [29; 30; 31; 32]. However, manual PIFs extraction is more accurate in some applications, especially when ground types are complicated.