Calibrating an airborne linear-array multi-camera system on the master focal plane with existing bundled images

ABSTRACT

Integrating multiple hyperspectral linear-array cameras into a Linear-array Multi-camera System (LMS) enables an expanded Field of View (FOV) while preserving high spatial resolution. However, the constrained capabilities of tie points are weakened by the limited overlap and short baseline between sub-cameras, hindering the attainment of reliable solutions for the theoretically ideal relative calibration model. Additionally, the absence of stereoscopic imaging in the LMS introduces the challenge of acquiring dense Ground Control Points (GCPs) throughout the FOV during calibration. In this paper, for calibrating an advanced airborne hyperspectral LMS named Arial Hyperspectral System-1000 (AHS-1000), we propose an equivalent geometric model that transforms the three-camera configuration of the LMS into a single-camera representation on the master camera’s focal plane, which was successfully applied to the data processing system of AHS-1000. We extracted dense GCPs by matching each Level-1 sub-image of the LMS with the bundled images of an existing three-line-scanner aerial image dataset, whose orientation parameters have been accurately solved by bundle adjustment with sparse GCPs; the matched points of each sub-image were consistently projected onto the master focal plane as image observations for our equivalent model. Finally, we established a collinearity equation system and performed a combined strip adjustment to calibrate the distortion terms for each pixel of the equivalent camera. The results demonstrate the effectiveness and robustness of our equivalent modeling method, which outperforms the classic modeling method by 2.83 pixels in interior accuracy and 0.86 pixels in relative accuracy, as well as outperforming another similar equivalent modeling method by 1.22 pixels in interior accuracy and 0.18 pixels in relative accuracy. Additionally, our comparison utilizing bundled images, public DOM (Digital Orthophoto Maps) /DEM (Digital Elevation Models), and high-precision DOM/DEM from aerial three-line-scanner demonstrates that bundled images can achieve higher calibration accuracy.

KEYWORDS:

• Linear-array multi-camera system
• sensor calibration
• equivalent geometrical model
• existing bundled images

1. Introduction

To enhance the efficiency of data collection and processing, airborne hyperspectral sensors are anticipated to possess high spectral resolution, high spatial resolution, and a large Field of View (FOV) concurrently (Jia et al. 2020). As an alternative to intricate and challenging solutions in optical design aimed at overcoming spatial limitations of individual hyperspectral sensors, the assembly of multiple hyperspectral linear-array cameras into a Linear-array Multi-camera System (LMS) is often considered a more favorable approach for expanding the FOV while maintaining high spatial resolution (Hu et al. 2005; Zhang et al. 2019). Due to independent data collection, the sub-cameras often suffer from geometric misalignments, necessitating the accurate recovery of their Relative Orientation Parameters (ROP) for seamless image stitching. Furthermore, to ensure both seamless stitching and precise georeferencing, it is essential to perform calibration for the LMS, taking into account the variations in the Interior Orientation Parameters (IOP) of the sub-cameras caused by their inherent geometric tolerances.

Achieving overall optimization of both interior accuracy and relative accuracy during the LMS calibration is a challenging task. Theoretically, accomplishing reliable LMS calibration involves establishing relative orientation constraints between sub-cameras based on corresponding points, as well as imposing interior accuracy constraints on each sub-camera using Ground Control Points (GCPs) (Gruen and Zhang 2002; Hartley 1997; Nister 2004). Nonetheless, the limited overlap between adjacent sub-cameras presents difficulties in accurately matching corresponding points, and uncertainty in the object-space coordinates of tie points increases due to the small baseline, prompting the development of equivalent geometric models in many studies (Cheng et al. 2018; Jiang et al. 2017). These models treat the misalignments caused by ROP as part of errors in IOP and compensate for them by utilizing the constraints provided by GCPs. Despite this progress, quantitative comparisons between these equivalent models and ideal models remain scarce. Additionally, it is essential to simultaneously formulate and solve the distortion parameters for all sub-cameras, rather than pursuing separate step-by-step approaches. In the case of the 5-camera system’s geometric stitching on the HY-1C satellite, Cao et al. (2020) first calibrated the IOP and the Exterior Orientation Parameters (EOP) of the master camera and then determined the IOP of the slave cameras. However, there is often a strong correlation between the IOP and EOP for linear-array imagery (e.g. the position of the array along the flight direction correlates with the pitch angle); consequently, this sequencing approach can result in the EOP obtained from bundle adjustment being locally optimal only for the master camera, rather than representing the overall optimal solution for the multi-camera system.

Another challenge in LMS calibration is the specific requirements for GCP density and distribution. Unlike Three-Line-Scanner (TLS) imaging systems, which can perform interior orientation calibration with a limited number of GCPs and utilize matching points from multi-view images, LMS lacks stereoscopic imaging conditions, necessitating denser GCP coverage across the entire FOV for accurate calibration (Leprince, Muse, and Avouac 2008). An alternative approach is to establish calibration strips with large overlaps; however, this method could significantly increase the complexity of data acquisition and calibration processing (Lenz et al. 2014). Therefore, in situations where high-cost ground calibration fields are unavailable, leveraging existing data or products with comparable or higher resolution imagery to generate dense digital GCPs becomes a natural consideration. While digital GCPs have been widely employed in the on-orbit calibration of spaceborne sensors (Chen et al. 2015; Tao et al. 2014; Wang et al. 2017), their utilization in aerial LMS calibration remains underexplored. The effectiveness and comparative analysis of different strategies for obtaining digital GCPs, such as extracting them from original images after bundle adjustment or from Digital Orthophoto Maps (DOM) and Digital Elevation Models (DEM) products, still lack research and clarity.

In this paper, we focus on the calibration problem of the Arial Hyperspectral System-1000 (AHS-1000), which is a recently modified version of a Chinese-made airborne hyperspectral LMS (Wang et al. 2016). Building upon the concept of equivalent geometric modeling, we transform the three-camera configuration of the LMS into an equivalent representation of a single camera by designating the intermediate camera as the master and projecting the slave cameras’ linear arrays onto its focal plane. Furthermore, since there was a recent aerial photogrammetry mission conducted by a TLS, we set the calibration data collection area to be consistent with the range of the TLS mission. This allows us to utilize the TLS images after bundle adjustment as geo-referencing data. Specifically, we first project the original sub-images onto the plane with the mean terrain height of the imaging area using the Position and Orientation System (POS) data, resulting in Level-1 images that overcome the geometric deformation caused by airborne platform jitter. Then, we extract dense corresponding points from TLS images for each Level-1 sub-image by image matching, the corresponding points located at forward, nadir, and backward images of TLS are further utilized to generate digital GCPs. Under the constraints of these GCPs, we reproject the matched points for the three Level-1 sub-images consistently onto the master focal plane and use their projected coordinates as image observations for the equivalent geometric model. Finally, we establish a collinearity equation system based on these observations and perform a combined strip adjustment, solving the EOP of the LMS and calibrating the distortion terms for each pixel of the equivalent camera. The proposed method has been applied to the AHS-1000 data processing system and validated for its effectiveness.

The experimental results validate the effectiveness and robustness of our proposed approach, which outperforms the ideal relative orientation calibration method by 2.83 pixels in interior accuracy and 0.86 pixels in relative accuracy, as well as outperforms the equivalent modeling method used by Cao et al. (2020) by 1.22 pixels in interior accuracy and 0.18 pixels in relative accuracy. Furthermore, we performed a confirmatory analysis between the proposed georeferencing strategy, which utilizes existing bundled images, and the strategy employing DOM/DEM from public data sources or aerial TLS data. The comparison revealed that the proposed strategy exhibited a more concentrated distribution of GCP residuals, resulting in a more accurate calibration outcome. Additionally, we evaluated the georeferencing accuracy of the LMS after calibration, showing a reduction in checkpoint error along the sample direction from the original 1.52 pixels to 0.66 pixels. Our main contributions can be summarized as follows:

1. By projecting the linear arrays of the slave cameras onto the focal plane of the master camera, we model the LMS as an equivalent single camera for geometric calibration, which facilitates holistic optimization of both interior and relative accuracy.

2. We apply the proposed equivalent geometric model to a newly developed airborne hyperspectral LMS, AHS-1000, and experimentally compare it with the ideal relative calibration model and a similar equivalent modeling method, thus validating its advantages and stability.

3. By utilizing the existing bundled images and their orientation parameters for LMS calibration, we qualitatively demonstrate that this georeferencing strategy outperforms the accuracy obtained from relying on DOM and DEM.

The remainder of this paper is organized as follows: Section 2 reviews the related work. Section 3 introduces the imager AHS-1000 and the experimental datasets. Section 4 details the proposed calibration method. Section 5 presents and discusses the results. Section 6 draws the conclusion.

2. Basic theory of algorithms

2.1. LMS calibration

LMS calibration primarily addresses internal distortions in sub-cameras and relative misalignments between sub-images resulting from boresight angles and lever-arm effects. Generally, recent research has placed greater emphasis on enhancing relative accuracy to ensure the production of seamless and accurate geo-rectified image products. Researchers have proposed two main categories of approaches: the misalignment compensation methods from image space and the misalignment calibration methods from object space. The image-space methods typically employ the shift model (Li, Liu, and Wang 2010), affine transformation model (Jacobsen 1999), or piecewise polynomial model (Meng et al. 2013) for image alignment and seamless stitching. However, this approach heavily relies on establishing reliable correspondences between adjacent sub-images, which is challenging when dealing with limited overlap or poorly textured regions.

The object-space misalignment calibration methods can be further classified into the classic relative orientation modeling method and the equivalent modeling method. The classic method is typically performed after completing the internal distortion calibration of sub-cameras. It focuses on rigorously calibrating the accurate angles and spatial offsets between adjacent cameras by simultaneously utilizing image-space correspondences and constraints from object space (Lichti et al. 2020). This method is achieved by minimizing the reprojection error of correspondences, which are projected from the original image coordinate system into the image coordinate system of reference camera based on initial relative parameters. It has been proved effective for frame multi-camera system calibration (Guo et al. 2021; Rameau et al. 2022). However, the narrow overlap design of LMS systems makes it difficult to obtain correspondences as reliably as frame cameras do, and the small baseline makes it even harder to calibrate the exact relative parameters under high flight altitude. While this method has been reported to achieve good calibration results on simulated aerial hyperspectral LMS data (Li et al. 2021), the challenges in obtaining reliable and evenly distributed matching points within narrow overlap areas could still introduce calibration biases.

In comparison, the equivalent modeling method aims to represent multiple imagers as a unified approximate imaging model, bypassing the need for rigorous ROP calibration between adjacent imagers. Instead, it focuses on indirectly addressing the internal distortions and relative orientation of each image based on this approximate model. Drawing upon the calibration research of spaceborne sensors with multiple linear arrays mounted on the same focal plane, such as ALOS PRISM (Kocaman and Gruen 2008) and ZY-3 (Jiang et al. 2014), some researchers constructed the approximate model based on considering the pointing angle of each pixel on a normalized plane. Jiang et al. (2017) conducted calibration on dual cameras using the pointing angle model and compensated for internal distortions and misalignment using a polynomial model. Nevertheless, this study overlooked the resolution changes brought about by tilting as the angle between the two camera optical axes was only 2°.

In a recent study conducted by Cao et al. (2020), which reported on the calibration of the HY-1C ultraviolet imager, they also employed the concept of equivalent modeling by representing the 5-camera system as a distorted virtual camera. However, there are significant disparities between their work and ours in terms of the calibration model and adjustment approach. Firstly, we diverge from their assumption of a shared optical center for the sub-cameras and incorporate the lever arm in our LMS modeling. Secondly, instead of adopting a two-step calibration process for the master camera and slave cameras, we integrate all sub-image GCPs into a single bundle adjustment procedure, which offers a more advantageous means of obtaining a globally optimal solution for the EOP.

2.2. Control data for sensor calibration

Control data, which is a set of GCPs, plays a crucial role in sensor calibration (Shan et al. 2023; Wang et al. 2012). Traditional control data are typically measured manually from the calibration field (Aguilar, Saldaña, and Aguilar 2013). However, collecting a sufficient number of GCPs to cover the entire range of an LMS, can be a time-consuming and labor-intensive task (Cheng et al. 2018). To overcome these limitations, other researchers propose to use internal constraints created by tie points among multiple overlapped images (Pi et al. 2019; Pi, Li, and Yang 2020; Yang et al. 2020). However, relying solely on tie points compromises the accuracy of the exterior calibration results and may still require sparse GCPs to achieve optimal external accuracy (Cao et al. 2021).

Another commonly used control data for sensor calibration is digital GCP obtained through matching with digital photogrammetric products. Many spaceborne sensor calibration researches derived digital GCPs from public geographic sources, such as the Shuttle Radar Topography Mission (SRTM) DEM or images from map service providers like Google Earth (Cao et al. 2022; Tao et al. 2014). However, the public geographic data contain local geometric deformation (Ruzgienė, Xiang, and Gečytė 2011), leading to inaccuracies in high-resolution airborne LMS calibration. Other studies requested large-scale and more accurate DOM and DEM from official geographic institutions, such as the National Geomatics Center of China (Jiang et al. 2017; Zhang et al. 2011; Zheng et al. 2015). However, height error would be introduced since many digital GCPs are typically feature points located on ground objects that may have changed or been removed in DEM. While the 0.5 m satellite DOM and the 1” Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM have been reported to succeed in supporting calibration on simulated aerial hyperspectral LMS data (Li et al. 2021), the practical efficacy of this georeferencing strategy on the AHS-1000 requires validation.

Considering that the existing DOM and DEM products may be limited in terms of spatial resolution or data quality, instead of directly using these final products to extract GCPs, the original images and the intermediate block bundle adjustment results can provide more originated and precise information for georeferencing (Zhang and Tao 2017). proposed to use the bundle adjustment result of existing images to obtain GCPs for aerial triangulation of newly collected images. However, the application of this idea is still limited to improving automatic triangulation efficiency and has not been studied and applied to LMS calibration yet. Besides, they only evaluated the advantage of existing bundled images over DOM/DEM based on experience without performing a comparison experiment.

3. Materials

3.1. AHS-1000 imaging system

We apply the proposed calibration method to AHS-1000, a modified version of the airborne imaging spectrometer introduced by Wang et al. (2016). As shown in Figure 1, following the previous design, AHS-1000 comprises three subsystems designed to collect data in three spectral ranges: Visible-Near Infrared (VNIR), Short-Wave Infrared (SWIR), and Long-Wave Infrared (LWIR). Each subsystem consists of three cameras of the same spectral range, assembled using a parallel light tube optical system to ensure that their primary optical axes intersect at the same point.

Figure 1. Camera configuration and working principle of the AHS-1000 imager.

The main parameters of the three types of sub-camera in AHS-1000 are shown in Table 1. The AHS-1000 is equipped with a total of 1400 spectral bands, allowing for comprehensive spectral coverage. The VNIR, SWIR, and LWIR subsystems achieve Ground Sampling Distances (GSDs) of 0.5 m, 1 m, and 1.2 m, respectively, at a flight altitude of 2000 m. Each sub-camera within the subsystems has a unique viewing angle while maintaining a consistent FOV of approximately 14.7°, forming a composite FOV of around 40°. The VNIR sub-camera utilizes Charge-Coupled Device (CCD) sensor arrays, while the SWIR and LWIR sub-cameras employ Mercury-Cadmium-Telluride Focal-Plane Arrays (MCT FPA).

Table 1. Main parameters of the three types of sub-camera of AHS-1000 (Wang et al. 2016).

Item	VNIR	SWIR	LWIR
Spectral range ($\mathit{\mu }\text{m}$)	0.4–0.95	0.95–2.5	8–12.5
FOV (°)	14.7	14.7	14.7
Detector/Array length (pixel)	CCD/1024	MCT/512	MCT/256
Pixel size ($\mathit{\mu }\text{m}$)	32	25	30
Focal length (mm)	128	50	50
Number of bands	512	512	376

The AHS-1000 incorporates the POS/AV 610, a system produced by Applanix, which records positions and attitudes at a rate of 100 Hz during data collection. The POS achieves horizontal and vertical positioning accuracies of approximately 1.5 m and 3 m, respectively, while demonstrating an attitude angle determination accuracy of around 0.005 deg. Given the identical assembly design of the sub-cameras in the three subsystems of AHS-1000, we have employed the proposed calibration method to calibrate each subsystem. Our specific focus was on the VNIR subsystem to showcase the effectiveness of the proposed calibration method. In the VNIR subsystem, the middle camera’s images overlap with approximately 140 pixels from the right camera and 100 pixels from the left camera. Considering that the 512 bands of the VNIR image share the same optical mechanism and the image texture should be clear for matching or measuring, we selected the 72nd band with a wavelength of 0.56 µm for the experiments.

3.2. Experimental dataset

In formulating the calibration plan for the AHS-1000, we scheduled a data collection mission in July 2020 to cover the Yushu area in Jilin Province, China, with the goal of minimizing the cost of acquiring dense GCPs. This decision was made based on the availability of an existing image dataset from Yushu that we had performed bundle adjustment in May 2020. This existing dataset was collected by another airborne mapping camera called AMS-3000 (Xi and Duan 2020), which is similar to the three-line-scanner ADS series (Hinsken et al. 2002).

Figure 2 shows the coverage and overview of the experimental datasets. Three distinct AHS-1000 datasets, labeled as datasets A, B, and C, are presented. Each dataset consists of six image strips. The blue rectangle in Figure 2(a) illustrates the coverage of these datasets, while Figure 2(b) provides an overview. Dataset A is designated as the calibration dataset and is used for performance comparisons with different methods. Datasets B and C serve for extension georeferencing validation of the proposed method. Datasets A and B primarily cover urban and rural areas, whereas dataset C predominantly encompasses fields and rivers. The image GSDs for datasets A, B, and C are 0.7 m, 0.7 m, and 1.0 m, respectively. For georeferencing validation experiments, we measured 28 GCPs using a Global Navigation Satellite System (GNSS) receiver in Real Time Kinematic (RTK) mode, which were called RTK GCPs. The locations of these GCPs are shown in Figure 2(a), 16 of them are designated as circular markers and are utilized as checkpoints to evaluate the accuracy of the bundle adjustment process.

Figure 2. (a) Coverage of the experimental datasets and GCPs locations; (b) overview of the AHS-1000 datasets.

In Figure 2(a), the green rectangle represents the coverage of the existing bundled images obtained from the AMS-3000. These images consist of four strips, each comprising three images captured in the forward, nadir, and backward directions. To enhance the accuracy of the exterior orientation parameters, 8 RTK GCPs were collected, followed by a bundle adjustment process by using 5 of the RTK GCPs as control points while the rest as checkpoints. The resulting checkpoints demonstrate horizontal and vertical accuracies of 0.06 cm and 0.14 cm, respectively, satisfying the requirements for 1:1000 scale production (Xi and Duan 2020).

4. Methodology

The proposed method revolves around the core concept of constructing an equivalent geometric model based on a master focal plane to perform a joint calibration of the LMS. In this method, the middle camera is designated as the master, while the other cameras serve as slaves. The linear arrays of the slave cameras are projected onto the focal plane of the master camera. By employing a post-compensation of the reprojection errors of image point observations after bundle adjustment, the internal distortion of the sub-cameras and the misalignment, coupled with resolution changes induced by the boresight angles and lever arm, are jointly calibrated as equivalent distortion on the master focal plane. The overall technical process, as illustrated in Figure 3, can be divided into three key modules: 1) Level-1 sub-image generation, 2) GCP matching, and 3) Calibration based on master focal plane.

Figure 3. The overall technical process of the proposed method.

4.1. The equivalent single-camera representation

To ensure the stability of geo-referencing accuracy, an imaging system consisting of three cameras needs to follow two requirements: synchronized exposure and rigid assembly. These two requirements determine that the sub-images have the same lines and the boresight angles and spatial offsets are constant between sub-cameras. Due to the boresight angle and lever arm between adjacent cameras, ground objects in overlapping areas are observed at different times. We label the left, middle, and right cameras as 1, 2, 3, in that order, and the strict geometric model can be described by Equation (1)(1) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\text{break}\left({\mathit{R}}_{S,i}^M{\left[\begin{array}{c}x-x_0+\mathrm{\Delta }\mathit{x}\\ \mathit{y}-y_0+\mathrm{\Delta y}\\ -\mathit{f}\end{array}\right]}_{\mathit{i}}+{\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]}_{\mathit{i}}\right)$(1) :

(1) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\text{break}\left({\mathit{R}}_{S,i}^M{\left[\begin{array}{c}x-x_0+\mathrm{\Delta }\mathit{x}\\ \mathit{y}-y_0+\mathrm{\Delta y}\\ -\mathit{f}\end{array}\right]}_{\mathit{i}}+{\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]}_{\mathit{i}}\right)$(1)

where ${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_{\mathrm{E}\mathrm{C}\mathrm{E}\mathrm{F}}^{\mathrm{T}}$ is the Earth-Centered Earth-Fixed (ECEF) coordinate of a ground point; ${\left[\begin{array}{ccc}{X}_S\left(t_j\right)& Y_S\left(t_j\right)& Z_S\left(t_j\right)\end{array}\right]}^{\mathrm{T}}$ represents the perspective center of LMS in the ECEF frame at time $t_j$; $\mathit{R}\left(t_j\right)$ is the rotation matrix from the camera coordinate system to the ECEF coordinate system at time $t_j$;${\mathit{R}}_u$ is the rotation matrix between the LMS and the Inertial Measurement System (IMS); ${\mathit{R}}_{S,i}^M$ is the rotation matrix from the $\mathit{i}$-th sub-camera to the master camera, $\mathit{i}=1,2,3$; ${\left[\begin{array}{ccc}{T}_x& T_y& T_z\end{array}\right]}_{\mathit{i}}^{\mathrm{T}}$ is the translation vector from the $\mathit{i}$-th camera to the master camera; ${\left(\mathit{x},\mathit{y}\right)}_i$ is the coordinate of the ground point in focal plane frame of the $\mathit{i}$-th camera; ${\left(x_0,y_0\right)}_{\mathit{i}}$ is the coordinate of principal point of camera $\mathit{i}$; ${\left(\mathrm{\Delta x},\mathrm{\Delta y}\right)}_i$ is the internal distortion parameters of the $\mathit{i}$-th camera; $f_i$ is the focal length of camera $\mathit{i}$; $\mathit{\lambda }$ is the scale factor. The 2nd camera is set as the master camera, and the 1st and 3rd cameras serve as slave. Thus, ${\mathit{R}}_{S,2}^M$ should be set as the identity matrix and ${\left[\begin{array}{ccc}{T}_x& T_y& T_z\end{array}\right]}_2^{\mathrm{T}}$ should be equal to 0.

The translational offsets between the LMS and the GNSS antenna, and the boresight misalignment between the axes of the IMS and the LMS are estimated during bundle adjustment (Gruen and Zhang 2002; Zhang et al. 2014). Our method aims at calibrating the internal distortions of LMS. Traditionally, to calibrate the LMS, control points on each sub-image and corresponding points serving as constrain between master and slave images are required. However, for airborne LMS like AHS-1000, the jitter of the platform and the limited overlapping area make it difficult to match adequate and accurate corresponding points. Besides, the baseline between the adjacent imagers is too small to support reliable forward intersection for getting coordinates of tie points. To avoid matching points between the master and slave cameras, an equivalent geometric model based on the master focal plane is established.

As shown in Figure 4(a), $P_M-\mathit{XYZ}$ and $P_S-\mathit{XYZ}$ are the master camera coordinate system and slave camera coordinate system, respectively. The transformation from the slave camera to the master camera is defined by ${\mathit{R}}_S^M$ and $\mathit{T}$. $o_M$ and $o_S$ are the principal point of the master camera and slave camera, respectively. $f_M$ and $f_S$ are the focal length of the master camera and slave camera, respectively. At the exposure time $t_j$, pixel Q of the slave camera can be projected to point $\mathit{Q}{ }^{\mathrm{\prime }}$ in the master focal plane. Thus, the coordinate $\left(x_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}},y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}\right)$ of each projected pixel of slave camera in the focal plane of master camera can be expressed as:

Figure 4. The sketch map of the imaging geometry between the master and slave cameras: (a) overview of imaging geometry, (b) enlarged rotation relations between camera S and M.

(2) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\frac{a_1{x}_S+a_2{y}_S-a_3{f}_S+T_x}{a_7{x}_S+a_8{y}_S-a_9{f}_S+T_z}\\ \\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\frac{a_4{x}_S+a_5{y}_S-a_6{f}_S+T_y}{a_7{x}_S+a_8{y}_S-a_9{f}_S+T_z}\end{array}$(2)

where $a_m\left(\mathit{m}=1,\dots ,9\right)$ are the elements of rotation matrix ${\mathit{R}}_S^M$; $\left(x_S,y_S\right)$ is the coordinate of a pixel in slave camera focal plane, and $\left(x_S,y_S\right)=\left(\mathit{x}-x_0+\mathrm{\Delta x},\mathit{y}-y_0+\mathrm{\Delta y}\right)$.

Since the slave and master cameras are assembled using a parallel light tube optical system to ensure that their primary optical axes intersect at the same point and their X axes parallel, we only consider the roll and a small yaw angle between adjacent cameras. Thus, the rotation matrix ${\mathit{R}}_{S,i}^M$ can be expressed as Equation (3)(3) ${\mathit{R}}_{S,i}^M={\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\mathit{\beta }& -sin\mathit{\beta }\\ 0& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i{\left[\begin{array}{ccc}cos\mathit{\tau }& -sin\mathit{\tau }& 0\\ sin\mathit{\tau }& cos\mathit{\tau }& 0\\ 0& 0& 1\end{array}\right]}_i\approx {\left[\begin{array}{ccc}1& -\mathit{\tau }& 0\\ \mathit{\tau }cos\mathit{\beta }& cos\mathit{\beta }& -sin\mathit{\beta }\\ \mathit{\tau }sin\mathit{\beta }& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i$(3) :

(3) ${\mathit{R}}_{S,i}^M={\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\mathit{\beta }& -sin\mathit{\beta }\\ 0& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i{\left[\begin{array}{ccc}cos\mathit{\tau }& -sin\mathit{\tau }& 0\\ sin\mathit{\tau }& cos\mathit{\tau }& 0\\ 0& 0& 1\end{array}\right]}_i\approx {\left[\begin{array}{ccc}1& -\mathit{\tau }& 0\\ \mathit{\tau }cos\mathit{\beta }& cos\mathit{\beta }& -sin\mathit{\beta }\\ \mathit{\tau }sin\mathit{\beta }& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i$(3)

where ${\beta }_i$ and ${\tau }_i$ are the roll angle and small yaw angle between the $\mathit{i}$-th camera and master camera, respectively.

For the linear-array camera, $x_S\approx 0$; Equation (2)(2) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\frac{a_1{x}_S+a_2{y}_S-a_3{f}_S+T_x}{a_7{x}_S+a_8{y}_S-a_9{f}_S+T_z}\\ \\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\frac{a_4{x}_S+a_5{y}_S-a_6{f}_S+T_y}{a_7{x}_S+a_8{y}_S-a_9{f}_S+T_z}\end{array}$(2) can be written as Equation (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) after combining Equation (3)(3) ${\mathit{R}}_{S,i}^M={\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\mathit{\beta }& -sin\mathit{\beta }\\ 0& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i{\left[\begin{array}{ccc}cos\mathit{\tau }& -sin\mathit{\tau }& 0\\ sin\mathit{\tau }& cos\mathit{\tau }& 0\\ 0& 0& 1\end{array}\right]}_i\approx {\left[\begin{array}{ccc}1& -\mathit{\tau }& 0\\ \mathit{\tau }cos\mathit{\beta }& cos\mathit{\beta }& -sin\mathit{\beta }\\ \mathit{\tau }sin\mathit{\beta }& sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right]}_i$(3) and simplification:

(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4)

where $tan{\psi }_{\mathit{ys}}=y_s/f_S$; ${\psi }_y$ is the look angle of each pixel of slave camera.

As an empirical mathematical model, the look angle model is widely used in the internal parameter calibration of line array cameras (Breton et al. 2002; Cao et al. 2021; Jiang et al. 2014). It unifies lens distortion, principal point offset and focal length bias with high-order polynomials to avoid over-parameterization. Based on the laboratory calibration values, internal parameter $tan{\psi }_{\mathit{ys}}$ of each sub-camera can be expressed as Equation (5)(5) $tan{\psi }_{\mathit{ys}}=tan{\psi }_{\mathit{ys},\mathit{r}}+a_0+a_1\mathit{s}+a_2{s}^2+\dots +a_j{s}^j,\mathit{j}\le 5$(5) .

(5) $tan{\psi }_{\mathit{ys}}=tan{\psi }_{\mathit{ys},\mathit{r}}+a_0+a_1\mathit{s}+a_2{s}^2+\dots +a_j{s}^j,\mathit{j}\le 5$(5)

where $tan{\psi }_{\mathit{ys},\mathit{r}}$ is obtained from laboratory calibration values; $\mathit{s}$ is the index number of the pixel; $a_0,\dots ,a_j$ are the model parameters for absorbing the influence of various geometric distortion factors, including lens distortion, principal point offset and focal length bias.

Equations (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) and (5(5) $tan{\psi }_{\mathit{ys}}=tan{\psi }_{\mathit{ys},\mathit{r}}+a_0+a_1\mathit{s}+a_2{s}^2+\dots +a_j{s}^j,\mathit{j}\le 5$(5) ) establish the relationship between the internal parameters of each sub-camera and the coordinates in the master focal plane, but they are relatively complex and are not conducive to effective calibration. Considering that Equation (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) is an inverse proportional function with respect to $tan{\psi }_{\mathit{ys}}$, we use Taylor expansion principle to simplify Equation (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) into a polynomial with respect to $tan{\psi }_{\mathit{ys}}$. In the case of AHS-1000, we found that sub-pixel approximation level can be achieved by performing first-order and third-order Taylor expansions on $x_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}$ and $y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}$ at the principal point, respectively. Thus, Equation (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) can be written as:

(6) $\begin{array}{rl}& \mathit{x}{S}^{\mathrm{\prime }}=f_M[c_{\mathit{s}0}+c_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)]\\ & \mathit{y}{S}^{\mathrm{\prime }}=f_M[d_{\mathit{s}0}+d_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)\\ & +d_2{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^2\\ & +d_3{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^3]\end{array}$(6)

where $c_{\mathit{s}0}$, $c_1$, $d_{\mathit{s}0}$, $d_1$, $d_2$ and $d_3$ are the coefficients after Taylor expansion, and they can describe the geometric deformation caused by misalignment between adjacent cameras; ${\psi }_{\mathit{ys}0}$ is the look angle of principle point.

In Equation (5)(5) $tan{\psi }_{\mathit{ys}}=tan{\psi }_{\mathit{ys},\mathit{r}}+a_0+a_1\mathit{s}+a_2{s}^2+\dots +a_j{s}^j,\mathit{j}\le 5$(5) , the order of $\mathit{s}$ depends on the accuracy of laboratory calibration values, as well as the severity of radial and the decentering distortions. When $\mathit{s}$ takes the 5th order, it will lead to terms up to the 15th order in Equation (6)(6) $\begin{array}{rl}& \mathit{x}{S}^{\mathrm{\prime }}=f_M[c_{\mathit{s}0}+c_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)]\\ & \mathit{y}{S}^{\mathrm{\prime }}=f_M[d_{\mathit{s}0}+d_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)\\ & +d_2{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^2\\ & +d_3{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^3]\end{array}$(6) , which may lead to over-parameterization. In the case of AHS-1000, taking first-order for $\mathit{s}$ is enough since the laboratory calibration values can achieve sub-pixel level. Thus, Equation (6)(6) $\begin{array}{rl}& \mathit{x}{S}^{\mathrm{\prime }}=f_M[c_{\mathit{s}0}+c_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)]\\ & \mathit{y}{S}^{\mathrm{\prime }}=f_M[d_{\mathit{s}0}+d_1\left(tan{\psi }_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)\\ & +d_2{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^2\\ & +d_3{\left(tan{\mathit{\psi }}_{\mathit{ys}}-tan{\psi }_{\mathit{ys}0}\right)}^3]\end{array}$(6) can be further simplified as Equation (7)(7) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[u_0+u_1\mathit{s}]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[v_0+v_1\mathit{s}+v_2{s}^2+v_3{s}^3]\end{array}$(7) :

(7) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[u_0+u_1\mathit{s}]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[v_0+v_1\mathit{s}+v_2{s}^2+v_3{s}^3]\end{array}$(7)

where $\left({\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}},{\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}\right)$ is the projected coordinate of each pixel of slave camera in the master focal plane; $u_0$, $u_1$, $v_0$, $v_1$, $v_2$ and $v_3$ are compensation parameters for misalignment between adjacent cameras and various geometric distortion factors of each sub-camera.

Thus, we can obtain the internal calibration terms for our equivalent geometric model, which can be expressed as Equation (8)(8) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=f_M[u_{0,\mathit{i}}+u_{1,\mathit{i}}\mathit{s}]\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=f_M[v_{0,\mathit{i}}+v_{1,\mathit{i}}\mathit{s}+v_{2,\mathit{i}}{s}^2+v_{3,\mathit{i}}{s}^3]\end{array}\mathit{i}=1,2,3$(8) :

(8) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=f_M[u_{0,\mathit{i}}+u_{1,\mathit{i}}\mathit{s}]\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=f_M[v_{0,\mathit{i}}+v_{1,\mathit{i}}\mathit{s}+v_{2,\mathit{i}}{s}^2+v_{3,\mathit{i}}{s}^3]\end{array}\mathit{i}=1,2,3$(8)

where $\left(\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}},\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\right)$ is the equivalent distortion of each sub-camera in the master focal plane; $u_{0,\mathit{i}}$, $u_{1,\mathit{i}}$, $v_{0,\mathit{i}}$, $v_{1,\mathit{i}}$, $v_{2,\mathit{i}}$ and $v_{3,\mathit{i}}$ are the calibration parameters; $\mathit{i}$ is the number of sub-camera.

In Equation (8)(8) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=f_M[u_{0,\mathit{i}}+u_{1,\mathit{i}}\mathit{s}]\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=f_M[v_{0,\mathit{i}}+v_{1,\mathit{i}}\mathit{s}+v_{2,\mathit{i}}{s}^2+v_{3,\mathit{i}}{s}^3]\end{array}\mathit{i}=1,2,3$(8) , $f_M$ is the focal length of master camera, which is only used for setting a reference plane. Since the master camera is also calibrated based on Equation (8)(8) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=f_M[u_{0,\mathit{i}}+u_{1,\mathit{i}}\mathit{s}]\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=f_M[v_{0,\mathit{i}}+v_{1,\mathit{i}}\mathit{s}+v_{2,\mathit{i}}{s}^2+v_{3,\mathit{i}}{s}^3]\end{array}\mathit{i}=1,2,3$(8) , and its focal length bias is absorbed in $u_{0,2}$, $u_{1,2}$, $v_{0,2}$, $v_{1,2}$, $v_{2,2}$ and $v_{3,2}$, $f_M$ will not be calibrated.

The influence of $\left(\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}},\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\right)$ can be seen in Figure 4(a). The projected sensor array of the sub-camera mainly suffers three system errors: 1) an offset of the projected principal point $o_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}$, expressed by $\left(f_M{u}_{0,\mathit{i}},f_M{v}_{0,\mathit{i}}\right)$; 2) a rotation angle $\mathit{\theta }$ with respect to the y-axis in master focal plane, caused by $f_M{u}_{1,\mathit{i}}$; 3) a resolution change along the y-axis, indicating variations in the size of the projected pixels due to $f_M{v}_{1,\mathit{i}},f_M{v}_{2,\mathit{i}}$ and $f_M{v}_{3,\mathit{i}}$. To make Equation (8)(8) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=f_M[u_{0,\mathit{i}}+u_{1,\mathit{i}}\mathit{s}]\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=f_M[v_{0,\mathit{i}}+v_{1,\mathit{i}}\mathit{s}+v_{2,\mathit{i}}{s}^2+v_{3,\mathit{i}}{s}^3]\end{array}\mathit{i}=1,2,3$(8) more concise, we partially define the calibration parameters, which is written as Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) :

(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9)

where $\left(\mathrm{\Delta }{x}_{0,\mathit{i}},\mathrm{\Delta }{y}_{0,\mathit{i}}\right)$ is the projected principal point offset of camera $\mathit{i}$ in the master focal plane; ${\theta }_i$ is the rotation angle of the $\mathit{i}$-th camera’s projected linear-array with respect to the y-axis in the master focal plane; ${\bar{y}}_{\mathit{i}}=y_i-y_{0,\mathit{i}}$, $y_i$ and $y_{0,\mathit{i}}$ are the y-axis coordinate of the projected pixel and that of projected principal point, respectively; $k_{0,\mathit{i}},k_{1,\mathit{i}}$ and $k_{2,\mathit{i}}$ are the deformation compensation coefficients.

It should be noted that Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) is accurate enough for compensating deformations in AHS-1000 caused by misalignment between adjacent cameras and various geometric distortion factors, because the yaw angle $\mathit{\tau }$(about 0.0004 rad) and baseline (below 9 cm) are small between adjacent camera. If the laboratory calibration values are not accurate, or the assembly process is not accurate enough, higher order in Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) should be considered.

If the boresight angles and spatial offsets between adjacent cameras are not measured in laboratory, the initial projected coordinates of each pixel of sub-camera in the master focal plane can be obtained by only considering the roll angle and set the spatial offsets as 0. The roll angle $\mathit{\beta }$ can be estimated according to the number of overlapping pixels. As shown in Figure 4(b), we applied a central symmetry transformation to the camera focal plane along the optical axis. The initial projected coordinate of $\mathit{Q}$ can be expressed as Equation (10)(10) $\begin{array}{rl}& x_{\mathit{s},\mathit{i}}=0\\ & y_{\mathit{s},\mathit{i}}=f_M\tan \left(\mathit{\alpha }+\mathit{\beta }\right)\mathit{with}\\ & \mathit{\alpha }=arctan\left(\frac{y_Q-y_{0,\mathit{s}}}{f_s}\right),\\ & \mathit{\beta }=0.5\mathit{FOV}+arctan\left(\frac{{\mathit{y}}_{\mathit{ov}}-y_{0,\mathit{M}}}{f_M}\right)\end{array}$(10) :

(10) $\begin{array}{rl}& x_{\mathit{s},\mathit{i}}=0\\ & y_{\mathit{s},\mathit{i}}=f_M\tan \left(\mathit{\alpha }+\mathit{\beta }\right)\mathit{with}\\ & \mathit{\alpha }=arctan\left(\frac{y_Q-y_{0,\mathit{s}}}{f_s}\right),\\ & \mathit{\beta }=0.5\mathit{FOV}+arctan\left(\frac{{\mathit{y}}_{\mathit{ov}}-y_{0,\mathit{M}}}{f_M}\right)\end{array}$(10)

where $y_Q$ is y-coordinate of $\mathit{Q}$; $y_{0,\mathit{s}}$ is y-coordinate of the principal point of slave camera; FOV is the field of view angle of the sub-camera; $y_{\mathit{ov}}$ is the y-coordinate of the last overlapping pixel; $y_{0,\mathit{M}}$ is y-coordinate of the principal point of master camera.

Based on Equation (1)(1) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\text{break}\left({\mathit{R}}_{S,i}^M{\left[\begin{array}{c}x-x_0+\mathrm{\Delta }\mathit{x}\\ \mathit{y}-y_0+\mathrm{\Delta y}\\ -\mathit{f}\end{array}\right]}_{\mathit{i}}+{\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]}_{\mathit{i}}\right)$(1) , Equation (7)(7) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[u_0+u_1\mathit{s}]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}={\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+f_M[v_0+v_1\mathit{s}+v_2{s}^2+v_3{s}^3]\end{array}$(7) and Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) , the equivalent geometric model is established as Equation (11)(11) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\left[\begin{array}{c}{\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}\\ {\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\\ -f_M\end{array}\right]$(11) :

(11) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\left[\begin{array}{c}{\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}\\ {\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\\ -f_M\end{array}\right]$(11)

where $\left({\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}},{\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}\right)$ is the projected coordinate of image point in master focal plane; $\left(\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}},\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\right)$ is the equivalent distortion terms defined by Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) .

4.2. Geometric correspondence determination

4.2.1. Level-1 sub-image generation

As shown in Figure 3, the original Level 0 sub-images exhibit significant geometric deformation, causing the ground objects in images to lose their shapes. Using these Level 0 sub-images directly for image matching would result in a low matching success ratio. To address this issue, data preprocessing is performed to generate rectified Level 1 sub-images by resampling the original sub-images onto a mean elevation plane.

During the geo-rectification process, all sub-images, including both the master and slave images, are processed using the same filtered EOP recorded by the POS. The geo-referencing of the images is performed using Equation (1)(1) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\text{break}\left({\mathit{R}}_{S,i}^M{\left[\begin{array}{c}x-x_0+\mathrm{\Delta }\mathit{x}\\ \mathit{y}-y_0+\mathrm{\Delta y}\\ -\mathit{f}\end{array}\right]}_{\mathit{i}}+{\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]}_{\mathit{i}}\right)$(1) . As depicted in Figure 3, the derived Level 1 sub-images do not exhibit the significant deformation caused by the jitter. It is important to note that the Level 1 sub-images are only used for matching corresponding points, and subsequent processes like bundle adjustment rely on the coordinates of the matching points on the master focal plane.

4.2.2. GCP matching

To estimate distortion terms at each pixel, a TLS utilizes stereo geometric constraints derived from tie points between multi-view images, often relying on a few sparsely distributed GCPs. However, in the case of the AHS-1000 VNIR subsystem, which functions as a one-line scanner, the availability of dense GCPs becomes crucial for achieving precise distortion estimation. The objective of GCP matching is to obtain dense digital control points that cover the entire FOV of the subsystem using the existing images.

A sequential matching process is used to obtain evenly distributed control points across the images. Given the average elevation, the range in the existing images corresponding to each AHS-1000 strip is first determined through forward and back-projection. As shown by the blue arrows in Figure 5, the second strip of AHS-1000 images corresponds to the yellow areas in existing bundled images, and the nadir images of the TLS dataset it covers will be matched to the test image sequentially. The matched points on the nadir images are further matched with the forward and backward images and selected as GCPs. The object-space coordinates of the GCPs are calculated through space intersection using the orientation parameters of the bundled images.

Figure 5. The matching strategy of obtaining dense digital GCPs from existing bundled images.

A dense point matching algorithm with two steps of expansion (Zhang et al. 2016) is adopted to acquire dense corresponding points which can smoothly operate on linear-array images without prior information. The number of GCPs can be adjusted by setting different matching intervals. 25,357 GCPs are obtained with a 33-pixel interval on dataset A in our case. Additionally, matching tie points between strip images of AHS-1000 is crucial for establishing constraints on strip relative orientation and aiding the bundle adjustment process, as shown by the red arrows in Figure 5.

4.3. Calibration based on master focal plane

After obtaining dense digital GCPs, a strip bundle adjustment is performed to refine the EOP of all scan lines and calculate the reprojection errors of the GCPs. To calibrate the equivalent distortion on the master focal plane, the image coordinates of the GCPs on the Level 1 images need to be transformed to the master focal plane as observations. This transformation process begins by mapping the points from the Level 1 images back to their corresponding pixels on the original sub-images. Subsequently, these points on the slave images are projected onto the slave focal plane and, ultimately, onto the master focal plane using Equation (11)(11) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\left[\begin{array}{c}{\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}\\ {\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\\ -f_M\end{array}\right]$(11) . To reduce the complexity of bundle adjustment, the Lagrange Interpolation Model (Gruen and Zhang 2002; Hinsken et al. 2002), which selects certain scan lines as orientation fixes at regular time intervals, is employed and only the EOP of those orientation fixes are solved in the least squares adjustment.

Since the equivalent distortion parameters defined by Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) are not self-calibrated during the bundle adjustment, the reprojection errors at GCP observations can provide valuable information about image distortion. By analyzing the distribution pattern of the reprojection errors, a post-compensation can be applied to fit the image distortion parameters for each sub-camera based on Equation (9)(9) $\begin{array}{c}\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{x}_{0,\mathit{i}}+{\bar{y}}_{\mathit{i}}\sin \left({\theta }_i\right)\\ \mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}=\mathrm{\Delta }{y}_{0,\mathit{i}}+k_{0,\mathit{i}}{\bar{y}}_{\mathit{i}}+k_{1,\mathit{i}}{\bar{y}}_{\mathit{i}}^2+k_{2,\mathit{i}}{\bar{y}}_{\mathit{i}}^3\end{array}\mathit{i}=1,2,3$(9) .

Given that the accuracies of the matched control points may vary, we incorporate a strategy in the adjustment process to calculate the reprojection errors of each control point after each iteration and update their respective weights accordingly. This approach ensures that control points with higher accuracy have a greater influence on the adjustment results, effectively improving the overall accuracy and reliability. Based on Equation (11)(11) ${\left[\begin{array}{c}\mathit{X}\\ \mathit{Y}\\ \mathit{Z}\end{array}\right]}_{\mathrm{ECEF}}=\left[\begin{array}{c}{X}_S\left(t_j\right)\\ Y_S\left(t_j\right)\\ Z_S\left(t_j\right)\end{array}\right]+\mathit{\lambda }\mathit{R}\left(t_j\right){\mathit{R}}_u\left[\begin{array}{c}{\hat{x}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{x}_{\mathit{s},\mathit{i}}\\ {\hat{y}}_{\mathit{S}{ }^{\mathrm{\prime }}}+\mathrm{\Delta }{y}_{\mathit{s},\mathit{i}}\\ -f_M\end{array}\right]$(11) , the distortion terms for each pixel on the master focal plane are determined by individually fitting curves to the residual errors of the GCPs from each sub-camera. To improve the accuracy of the curve fitting process, an additional step of noise point filtering can be applied, which is inspired by a density analysis method (Zhu et al. 2015). Our filtering approach is to divide the residual points into grids and analyze the point density in each. Noise grids have densities below a specific threshold, determined by the median of point density in each grid column, while grids exceeding this threshold are considered valid. Subsequently, we calculate the distance between the points in the noise grids and their nearest points in the valid grids. Points with distances less than half of the grid height are labeled as valid, while any remaining noise points within the noise grids are removed.

5. Results and discussion

5.1. Experimental settings

To evaluate the performance of the proposed calibration method, we conducted a comparative analysis with two methods: the classic relative orientation modeling method as utilized by Li et al. (2021), and the similar equivalent modeling method described by Cao et al. (2020). The classic method solves boresight angles and spatial offsets between sub-cameras during bundle adjustment, with the geometric constraint of the corresponding points between adjacent sub-images. Cao et al.’s method is performed based on the assumption that all the sub-cameras share the same optical center and a two-step calibration process for solving look-angles of each pixel of the master camera and slave cameras. For the corresponding point matching, a two-step expansion algorithm was employed (Zhang et al. 2016). In both the proposed method and the comparison methods, the same set of digital GCPs matched from the existing bundled images was used for calibration.

To assess the interior calibration accuracy of the calibrated images, we used geo-referencing accuracies of the final directly stitched image in the sample direction (${\sigma }_{\mathrm{s}\mathrm{a}\mathrm{mple}}^{\mathrm{i}\mathrm{n}}$), line direction (${\sigma }_{\mathrm{line}}^{\mathrm{i}\mathrm{n}}$), and planimetry (${\sigma }^{\mathrm{in}}$) as evaluation indicators. Meanwhile, the following steps were performed for evaluating the relative orientation accuracy. Firstly, the corresponding points were matched between adjacent sub-images after rectification; secondly, using the points on the master image as references, the corresponding points on slave images were projected to the ground and subsequently reprojected onto the master image. The Root Mean Square Error (RMSE) of the reprojection errors of these corresponding points in sample direction (${\sigma }_{\mathrm{s}\mathrm{a}\mathrm{mple}}^{\mathrm{r}\mathrm{e}}$), line direction (${\sigma }_{\mathrm{line}}^{\mathrm{r}\mathrm{e}}$) and planimetry (${\sigma }^{\mathrm{re}}$) were calculated as the relative accuracy.

Furthermore, we compared our geo-referencing strategy using bundled images with two others. One group inspired by the strategy of Li et al. (2021) utilized a 0.5 m public DOM from Google Earth with approximately 2 m geo-referencing accuracy (Goudarzi and Landry 2017; Pulighe, Baiocchi, and Lupia 2016) and the public 1” ASTER DEM with around 13 m vertical accuracy (Okolie et al. 2024; Tachikawa et al. 2011). The other group employed products of AMS-3000, including the 0.2 m DOM with 0.2 m geo-referencing accuracy and 15 m DEM with 2 m height accuracy, verified based on RTK GCPs.

Both the proposed method and the comparison methods were implemented using the C++ language. The experiments were conducted on a computer with an Intel i7 -10,750 H CPU and 32GB RAM.

5.2. Overall comparison

As shown in Figure 6(a), 25783 GCPs are obtained from AMS-1000 bundled images. These GCPs are uniformly distributed across the primary range of dataset A, facilitating pixel-by-pixel calibration for the AHS-1000. The distribution of tie points prepared for the classic modeling method is shown in Figure 6(b). The lengths of the six strips in dataset A all exceed 60,000 rows, and we have cropped a portion of each strip to clearly show the distribution of tie points. These tie points are distributed within the narrow overlap region of adjacent cameras, and the matching accuracy achieves sub-pixel level. The numbers of tie points in each strip image between adjacent cameras are listed in Table 2.

Figure 6. The distributions of GCPs and tie points in dataset A: (a) distribution of GCPs, (b) distribution of tie points in portions of the six strip images.

Table 2. The number of tie points between adjacent cameras in each strip image of dataset A.

	Number of tie points
	Strip 1	Strip 2	Strip 3	Strip 4	Strip 5	Strip 6	Total
Camera L-M	513	290	466	392	376	348	2385
Camera R-M	781	622	746	683	608	679	4119

5.2.1. Quantitative evaluation

The overall comparison of interior georeferencing accuracy and relative georeferencing accuracy is shown in Table 3. Cameras L, M, and R correspond to the left, middle, and right cameras, respectively.

Table 3. Comparison of interior and relative orientation accuracies.

Method	Interior error (pixel)			Relative error (pixel)
				Camera L-M			Camera R-M			Average
	${\sigma }_{\mathrm{s}\mathrm{a}\mathrm{mple}}^{\mathrm{i}\mathrm{n}}$	${\sigma }_{\mathrm{line}}^{\mathrm{i}\mathrm{n}}$	${\sigma }^{\mathrm{in}}$	${\sigma }_{\mathrm{s}\mathrm{a}\mathrm{mple}}^{\mathrm{r}\mathrm{e}}$	${\sigma }_{\mathrm{line}}^{\mathrm{r}\mathrm{e}}$	${\sigma }^{\mathrm{re}}$	${\sigma }_{\mathrm{s}\mathrm{a}\mathrm{mple}}^{\mathrm{r}\mathrm{e}}$	${\sigma }_{\mathrm{line}}^{\mathrm{r}\mathrm{e}}$	${\sigma }^{\mathrm{re}}$	${\bar{\sigma }}_{\mathrm{sample}}^{\mathrm{re}}$	${\bar{\sigma }}_{\mathrm{line}}^{\mathrm{re}}$	${\bar{\sigma }}^{\mathrm{re}}$
Classic modeling	2.63	2.72	3.78	2.15	1.54	2.64	0.65	0.43	0.78	1.40	0.99	1.71
Cao et al.’s	1.85	1.12	2.17	0.95	0.47	1.06	0.82	0.56	0.99	0.89	0.52	1.03
Equivalent modeling (ours)	0.68	0.64	0.95	0.93	0.33	0.99	0.54	0.44	0.69	0.75	0.39	0.85

Our equivalent modeling method exhibited superior relative accuracy when compared to the other two methods. Specifically, the relative accuracy between camera L and M in our method was 0.99 pixels, whereas the classic method achieved 2.64 pixels and Cao et al.’s method achieved 1.06 pixels. Furthermore, our method consistently achieved sub-pixel level accuracies for both camera L-M and camera R-M, while Cao et al.’s method did not reach sub-pixel level accuracy for camera L-M. The performance of the classic method varied greatly, with a planimetry relative accuracy of 0.78 pixels for camera R-M but 2.64 pixels for camera L-M. This discrepancy can be attributed to the differing number of overlapping pixels between the sub-cameras, with camera R-M having approximately 40 more overlapping pixels than camera L-M. Consequently, the classic method had approximately 1734 fewer corresponding points between camera L and M compared to the 4119 points between camera R and M in total. Furthermore, the baseline between the adjacent imagers is very small (below 9 cm). The base-height ratio can not support high-precision relative calibration which is based on the stereo constraints of tie points. It is worth noting that the relative accuracy of 0.43 pixels along the line direction of the classic method for camera R-M was the best, which was 0.01 pixels better than ours. Considering that our method was better than the classic method on relative accuracy along line direction for camera L-M by 1.21 pixels, the 0.01 pixels improvement furtherly indicates that the performance of the classic method heavily depends on the reliability of correspondences in overlapping area and is limited by the small baseline. On average, our method outperformed the classic method by 0.86 pixels and Cao et al.’s method by 0.18 pixels in relative accuracy. The higher average relative accuracy achieved by both our method and Cao et al.’s method compared to the classic method underscores the effectiveness of the equivalent modeling approach in attaining consistent calibration results for LMS among sub-cameras. This advantage is particularly notable because it does not depend on correspondences between sub-cameras and can even benefit from a small baseline, distinguishing it from the classic method. According to Equation (4)(4) $\begin{array}{c}{x}_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[\frac{-\mathit{\tau }}{sin\mathit{\beta }}+\frac{T_{x}sin\mathit{\beta }/f_S-\mathit{\tau }\left(cos\mathit{\beta }-T_z/f_S\right)}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\\ y_{\mathit{S}{\mathit{ }}^{\mathrm{\prime }}}=f_M\left[cot\mathit{\beta }+\frac{1-T_{z}cos\mathit{\beta }/f_S+T_{y}sin\mathit{\beta }/f_S}{tan{\psi }_{\mathit{ys}}{sin}^2\mathit{\beta }-cos\mathit{\beta }sin\mathit{\beta }+T_{z}sin\mathit{\beta }/f_S}\right]\end{array}$(4) , the smaller the baseline is, the higher the approximate accuracy the polynomials can achieve. On the other hand, it can be seen that all the metrics of our method were better than Cao et al.’s method. This can be attributed to our one-step calibration process, which can achieve global optimal calibration result especially for aerial LMS suffering from significant platform jitter. The strong correlations between interior parameters and exterior parameters may result in different calibration result if using the two-step calibration process as utilized by Cao et al.’s method when different calibration images are employed.

Regarding interior georeferencing accuracy, our method also outperformed the classic method and Cao et al.’s method, achieving a planimetry accuracy of 0.95 pixels compared to 3.78 pixels and 2.17 pixels. Indeed, over 25,000 GCPs were matched from the existing images for bundle adjustment of the three methods. For the classic method, due to the fact that the tie points between adjacent sub-images were much less than digital GCPs, the tie points had 100 times more weight than those of the GCPs to make the constraint between sub-images works. However, the uneven distribution of tie points along the sample or line direction, coupled with the challenge of determining the precise object-space coordinates due to the small baseline, resulted in a suboptimal relative calibration outcome. Differently, our proposed method and Cao et al.’s method mitigated the dependence on tie points in the overlapping area. The 1.22-pixel improvement in interior accuracy achieved by our method compared to Cao et al.’s method demonstrates that our approach effectively addressed the misalignment caused by boresight angles and spatial offsets. Therefore, our master focal plane equivalent model and one-step calibration process have proven effective in obtaining more precise and reliable calibration results.

5.2.2. Visual alignment performance in overlapping areas

To further illustrate the effectiveness of our method in improving relative accuracy, we visually compared and analyzed the alignment performance in overlapping areas of the rectified sub-images produced by the three methods. We performed the rectification based on the three methods with the same strip image and superimposed the left and right sub-images onto the middle sub-image, as shown in Figure 7(a). For ease of comparison, the middle images of the three methods are displayed in the red channel while the left and right images are displayed in green. Four blocks of the partially enlarged images, denoted as areas 1, 2, 3, and 4, are shown in Figure 7(b). Among the four areas, areas 2 and 4 lie in the overlapping area of camera R and M, while areas 1 and 3 lie in that of camera L and M.

Figure 7. Visual alignment performance comparison in overlapping areas of different methods: (a) overview of the test strip image and location of comparison blocks, (b) enlarged comparison blocks.

In Figure 7(b), it is evident that the classic method exhibits an offset of around 3 pixels along the line direction in areas 1 and 3, which is more significant than the approximately 1-pixel offset observed in areas 2 and 4. This suggests that the ROP solved between camera L and M had lower accuracy compared to those between camera R and M. As previously discussed in Section 5.2.1, the narrower overlap between camera L and M and the small baseline resulted in poorer relative calibration results. The performance of Cao et al.’s method is better than the classic method in areas 1 and 3 but slightly worse in areas 2 and 4. These observations indicate that Cao et al.’s method reduces the influence of narrow overlapping areas, but its achievable accuracy is limited due to the two-step calibration process and the assumption of a shared optical center.

For enlarged block images of ours, the visual performances in all four areas were better than that of the comparison methods; the offset was barely visible to the naked eye. The misalignment between the sub-cameras was effectively reduced, which indicated that the proposed method succeeded in compensating for errors caused by the boresight angles and spatial offsets between adjacent cameras. Furthermore, our method outperformed Cao et al.‘s method in mitigating misalignment in areas 2, 3 and 4, underscoring the advantages of our one-step calibration process and stringent equivalent assumption. Additionally, our method demonstrated consistent performance across cameras R, L, and M, showing resilience to narrow overlaps.

5.3. Advantage of using bundled images for geometric control

Figures 8 and 9 depict a comparison of residual error distributions at GCP observations obtained from public DOM/DEM, AMS-3000 DOM/DEM and bundled images before and after calibration. To facilitate comparison, the residual errors before and after the calibration process are represented by red and blue colors, respectively.

Figure 8. Residual error distributions at GCP observations obtained from public DOM/DEM (top row), AMS-3000 DOM/DEM (middle row), and bundled images (bottom row) before calibration. Residual errors of using public DOM/DEM along sample direction (a) and line direction (b). Residual errors of using AMS-3000 DOM/DEM along sample direction (c) and line direction (d). Residual errors of using bundled images along sample direction (e) and line direction (f).

Figure 9. Residual error distributions at GCP observations obtained from public DOM/ASTER DEM (top row), AMS-3000 DOM/DEM (middle row), and bundled images (bottom row) after calibration. Residual errors of using public DOM/DEM along sample direction (a) and line direction (b). Residual errors of using AMS-3000 DOM/DEM along sample direction (c) and line direction (d). Residual errors of using bundled images along sample direction (e) and line direction (f).

In Figures 8 and 9, the x-axis is the sample coordinate of the projecting sensor array on the master focal plane, while the y-axis represents the residual errors in sample or line directions after the bundle adjustment. To ensure clear visibility of trends and prevent the blending of residual errors in the overlapping area, pixels from both the master camera and the slave cameras are equally sampled in the overlapping region. Regarding the distribution of residual errors at GCP observations of the DOM/DEM and bundled images before calibration, depicted in Figure 8, a polynomial distorted pattern in the sample direction and a linear distortion pattern in the line direction are observed. These patterns can be divided into three segments corresponding to the left, middle, and right cameras, resulting from the inaccurate initial projection coordinates of pixels on the master focal plane.

Comparing the residual error distributions before and after calibration, our bundle images group and the AMS-3000 DOM/DEM group show significant advantages over the public DOM/DEM group. Firstly, our bundle images group and the AMS-3000 DOM/DEM group demonstrate a higher concentration of residual errors near the piecewise curves, whereas the public DOM/DEM group falls short in this aspect. Secondly, they avoid the distorted curve pattern observed in the public DOM/DEM-based method after calibration. The discrete distribution pattern of the public DOM/DEM leads to inaccurate curve fitting parameters, resulting in poor calibration accuracy using the post-compensation method. Although the misalignment error between sub-cameras is reduced in public DOM/DEM group, as shown in Figures 9(a,b), the bending trend of the distribution at the ends indicates that the scanners at the ends of the linear-array are not accurately calibrated. This discrepancy arises from calculating the coordinates of GCPs in the presence of non-negligible elevation errors, which means the 0.5 m satellite DOM and the 1” ASTER DEM are not accurate enough for AHS-1000 calibration.

As shown in Figures 9(c,b), distortions are generally eliminated after calibration using the AMS-3000 DOM/DEM, indicating that higher accuracy products can effectively improve the calibration result. While our bundled images group shows an advantage in concentrating residual errors near piecewise curves over the AMS-3000 DOM/DEM, this advantage is less pronounced compared to the public DOM/DEM, as depicted in Figures 9(e,f). Indeed, the RMSE of GCPs after calibration in our bundled images group is 0.47 m, outperforming the 0.55 m RMSE of the AMS-3000 products group and the 1.02 m RMSE of the public DOM/DEM group. The impact of the 2 m vertical error in AMS-3000 DEM is not as significant as expected, likely due to the relatively flat terrain of dataset A, where elevations range from 160.31 m to 243.23 m, with a mean elevation of 197.68 m. However, the absence of readily available suitable DOM and DEM in the database for camera calibration may necessitate new production tasks, thereby hindering efficiency and increasing costs. Given that each project saves the bundle adjustment results, utilizing bundled images for geo-referencing stands out for ensuring the highest calibration accuracy while also improving efficiency.

5.4. Georeferencing accuracy validation

To further assess the stability and accuracy of the calibration result, experiments were conducted with and without the use of RTK GCPs, as shown in Figure 2(a), on datasets B and C. When all the RTK GCPs were utilized as checkpoints, the reprojection errors were computed. In cases where several RTK GCPs were utilized as control points while the others were used as checkpoints, bundle adjustment was performed and the reprojection errors of the checkpoints were subsequently calculated based on the refined EOP. Both experiments were based on the calibration result obtained from dataset A. The RMSEs and the standard deviations of the checkpoints after calibration are listed in Table 4. For ease of comparison, the results of the bundle adjustment without calibration are also listed in Table 4.

Table 4. The georeferencing accuracy for bundle adjustment of AHS-1000 VNIR subsystem before and after calibration.

Dataset	Number of control points	Number of checkpoints	Group	Calibration	RMSE (pixel)		Standard deviation (pixel)
					sample	line	sample	line
B	0	12	1	before	3.50	26.68	2.49	6.06
				after	2.87	23.43	1.49	6.97
	5	7	2	before	1.52	0.58	0.93	0.43
				after	0.66	0.78	0.47	0.52
C	0	16	3	before	4.07	9.56	2.91	4.16
				after	2.70	8.78	1.60	4.10
	7	9	4	before	1.20	0.85	0.76	0.64
				after	0.68	0.72	0.43	0.49

As shown in Table 4, the use of sparse RTK GCPs for bundled adjustment in groups 2 and 4 resulted in higher accuracy compared to groups 1 and 3 without RTK GCPs, as expected. It is worth noting that the RMSE along the line direction in group 1 was significantly larger (over 20 pixels) than those in group 3 (less than 10 pixels), which can be attributed to potentially unfavorable flight conditions during data acquisition. Besides, the RMSE along the sample direction in groups 2 and 4 exhibited more significant improvement compared to the improvement along the line direction, respectively. This improvement difference can be attributed to the adjustment of the EOP, which effectively compensated for errors in the line direction. Figure 8 visually demonstrates the successful compensation of relative errors along the line direction through our proposed calibration method. Moreover, there was an improvement in standard deviation after calibration, indicating enhanced convergence during bundle adjustment.

Analyzing the bundle adjustment results with sparse RTK GCPs, the experiment in group 2 exhibited a notable improvement in RMSE along the sample direction from 1.52 pixels to 0.66 pixels after calibration, while group 4 showed improvement from 1.20 pixels to 0.68 pixels along the sample direction. These results provide evidence that our proposed method effectively calibrated the equivalent distortion terms, leading to improved accuracy.

6. Conclusions

In this paper, an equivalent geometric calibration method based on the master focal plane for the airborne LMS was proposed. We first projected the sensor array of slave cameras to the master focal plane. Then, original sub-images were rectified into Level 1 sub-images for correcting geometric deformation, which were used to match with existing bundled images to obtain corresponding points that cover the whole FOV of LMS. The ground coordinates of these points were calculated with the accurate orientation parameters of existing bundled images and used as control points in strip bundle adjustment. Finally, the equivalent distortions on the master focal plane were calibrated using a post-compensation method to fit the reprojection errors of GCPs. In the experiment, the VNIR sub-system of AHS-1000, which consisted of three independent cameras, was calibrated using existing bundled images of another airborne stereo mapping camera AMS-3000. The results demonstrate that our equivalent modeling method surpasses the classic method by 2.83 pixels in interior accuracy and 0.86 pixels in relative accuracy, and outperforms another similar equivalent modeling method by 1.22 pixels in interior accuracy and 0.18 pixels in relative accuracy. Besides, the correspondence offset along the line direction in narrow overlapping areas of the proposed method was reduced to the sub-pixel level, which was better than the comparison methods with visible offset. The comparison experiment between the georeferencing strategy using existing bundled images and DOM/DEM from public data sources or aerial TLS indicated that dense digital GCPs from existing bundled images can provide more accurate calibration result. The georeferencing accuracy validation experiment further revealed that the proposed calibration method can improve georeferencing accuracy along the sample direction from 1.52 pixels to 0.66 pixels.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The experimental dataset is constrained by a confidentiality agreement, preventing the public disclosure of the dataset at this stage.

Additional information

Funding

This work is supported by National High-score Major Special Projects [grant number 50-H31D01-0508-13/15], China Yunnan province Major Science and Technology Special Plan Project [grant number 202202AF080004], China Yunnan province Ma Hongqi Academician Free Exploration Project [grant number 202005AA160011], and China Huaneng Lancang River Hydropower Inc. Science and Technology Project [grant number HY2021/D11].

Notes on contributors

Hao Zhang

Hao Zhang received the B.E. degree from Wuhan University, Wuhan, China, in 2018. He is currently pursuing the PhD degree with the Wuhan University, China. His research interests include simultaneous localization and mapping, photogrammetry and machine learning.

Yansong Duan

Yansong Duan received the PhD degree from Wuhan University, Wuhan, China, in 2016. He is currently an associate professor with the School of Remote Sensing and Information Engineering, Wuhan University. He works in the fields of digital photogrammetry, multi-view images 3D reconstruction and machine learning.

Qi Zhou

Qi Zhou received the B.E. degree from Wuhan University, Wuhan, China, in 2021. She is currently a master at Wuhan University, China. Her research interests include image matching, remote sensing and photogrammetry.

Qi Chen

Qi Chen received the PhD degrees from Wuhan University, Wuhan, China, in 2015. He is currently an associate professor. His research interests include space and aerial photogrammetry and pattern recognition from remote sensing imagery with deep-learning techniques.

Bowen Cai

Bowen Cai received the B.E. degree in remote sensing science and technology from Wuhan University in 2018. He is currently pursuing his PhD degree with the School of Remote Sensing and Information Engineering, Wuhan University. His research interests include urban remote sensing, social ecology and remote sensing applications.

Pengjie Tao

Pengjie Tao is currently an associate professor. His research interests include satellite photogrammetry, registration of optical images and LiDAR points, and multi-view images 3D reconstruction.

Zuxun Zhang

Zuxun Zhang is a professor and major in digital photogrammetry and remote sensing. He was selected as a member of International Academy of Sciences for Europe and Asia in 1995 and a member of Chinese Academy of Engineering in 2003. He has led a number of national and ministerial major research projects and has achieved substantial outcomes in aviation image mapping automation.