In the satellite MVS task, our goal is to use an end-to-end coarse-to-fine framework to predict the height map H by leveraging the matching relationship between N − 1 adjacent views and the corresponding camera parameters. First of all, we extract the image features F i i = 0 N − 1 from the reference images I ₀ and source images I 0 i = 0 N − 1 . Then the cost volume V is constructed by the differentiable RPC warping based on hypothetical height planes D. Next, a regularization process is executed on cost volume by a 3D Unet. After regularization, the regularized cost volume V _re regresses a probability volume P by the softmax operation. Lastly, the final height map H is calculated by the hypothetical height planes D and probability volume P. At the inference stage, the trained model on the satellite MVS task must infer the depth maps of all views of all scenes. Finally, a depth map fusion method is used to obtain point clouds.

Our proposed A-SATMVSNet is an trainable framework, which consists of two import parts: feature extraction and cost volume construction. As shown in Figure 1, the N input images I i i = 0 N − 1 ∈ R H × W × 3 are sent to multi-scale feature extraction module. After feature extraction, the multi-scale feature maps F i = 0 N − 1 are fed into the cost volume C construction part in three stages. The cost volume C is constructed by the differentiable rpc warping (Section 2.5). Then, the obtain cost volume C are regularized to generate probability volumes P by the softmax operation. Finally, the height maps can be obtained through regression.

In this section, we mainly describe the proposed attention-aware multi-scale feature extraction module. There are many popular feature extraction modules such as UNet-based (Ronneberger et al., 2015; Isensee et al., 2018; Li et al., 2018; Oktay et al., 2018; Huang et al., 2020), feature pyramid network-based (Lin et al., 2017; Kim et al., 2018; Seferbekov et al., 2018; Zhao et al., 2021), resnet-based (He et al., 2016; Targ et al., 2016; Szegedy et al., 2017; Bharati et al., 2021), etc. All the above feature extraction modules perform well in multi-view stereo matching tasks. In our study, we propose a new feature extraction module based on (Cheng et al., 2020) where it is combined with an attention module. The basic module consists of an encoder and a decoder with skip connection. The module outputs a three-scale feature pyramid whose size is {1/16, 1/4, 1} of the input satellite image size, and the number of feature channels is 32, 16, and 8 respectively. In the encoder part, an attention module is designed. The attention module and feature extraction layer in encoder and decoder part are shown in Figure 2.

Figure 2A shows the architecture of the detailed feature extraction layer in encoder network with attention module. First, a convolution layer with 3 × 3 kernel size is used to extract features. After that, the feature map is sent to three different dilated convolution layers with dilation rate of 2, 3 and 4 respectively. Then, all the three output feature maps are sent to a 3 × 3 convolution layer with an attention module. Finally, the three output feature maps are concatenated to generate a new feature map, as the final feature map. The formulation of our triple dilate convolution is defined as follow: F out = F i n 1 ⊗ w 1 ⊙ I + w 1 ⊙ I , F i n 2 ⊗ w 2 ⊙ I + w 2 ⊙ I , F i n 3 ⊗ w 3 ⊙ I + w 3 ⊙ I (1) where ⊗ represents the multiply operation, ⊙ denotes the element-wise product, w _i represents the ith weights of dilate convolution.

Figure 2B shows the architecture of the attention module. The input feature map is defined as F _in . Two convolution layers with a kernel of 3 × 3 are employed to generate further features F i n _ 1 . Then, a sigmoid function is used to obtain attention weights defined as F _w . The final output feature is defined as F _out , which is calculated as: F out = F i n + F i n _ 1 ⊗ F w . (2)

Figure 2C is the architecture of the decoder network, which consists of a deconvolution layer with a kernel size of 3 × 3, stride size of 2 and a convolution layer with a stride size of 1.

The rational polynomial camera model (RPC) is extensively used in satellite imagery processing, which connects the image points and corresponding world coordinate points with cubic rational polynomial coefficients (Gao et al., 2021). We define the world coordinates as l a t n , l o n n , h e i n which represents the latitude, longitude and height. The corresponding normalized image coordinates are defined as s a m p n , l i n e n . P ^fwd and P ^inv are both cubic polynomials. The transformation between world coordinates and image coordinates are shown as bellow: l a t n = P 1 inv s a m p n , l i n e n , h e i n P 2 inv s a m p n , l i n e n , h e i n , l o n n = P 3 inv s a m p n , l i n e n , h e i n P 4 inv s a m p n , l i n e n , h e i n , s a m p n = P 1 fwd l a t n , l o n n , h e i n P 2 fwd l a t n , l o n n , h e i n , l i n e n = P 3 fwd l a t n , l o n n , h e i n P 4 fwd l a t n , l o n n , h e i n , P X , Y , Z = ∑ i = 0 m 1 ∑ j = 0 m 2 ∑ k = 0 m 3 c ijk ⋅ X i ⋅ Y i ⋅ Z i . (3)

In multi-view stereo matching task via satellite imagery, the RPC model is a widely used geometric model, which can provide a high accuracy to the rigorous sensor model (RSM) (Tao and Hu, 2001).

Currently, most state-of-the-art MVS methods warp the source views to a reference view to obtain per-view matching feature volumes by a homography matrix and a set of fronto-parallel depth hypotheses planes D. The definition of the differentiable homography based pin-hole camera model is as bellow: H i d = d K i T i T ref − 1 K ref − 1 , (4) where T and K denote camera extrinsic and intrinsic respectively. Compared with pin-hole camera model, the cubic rational polynomial camera (RPC) model is widely used in satellite domain, which has the advantage than all camera models, e.g., projective, affine and the linear pushbroom. A matrix alone cannot formulate the warping of the RPC model due to its complexity. In this regard, SatMVS proposes a rigorous and efficient RPC warping module that is fundamentally a high-order tensor transformation, which is fundamental to the structure of SatMVS. Using a set of hypothetical height planes in the world coordinate system, the RPC warping module projects images from different views to the reference view, instead of the fronto-parallel planes of a reference view, because the RPC model does not include explicit physical parameters for defining the front of a camera.

Firstly, SatMVS transforms the ternary cubic polynomial by using cubic polynomials to a quaternion cubic homogeneous polynomial f x 1 , x 2 , x 3 , x 4 = ∑ a i a j a k ⋅ x i x j x k , where a i a j a k i , j , k ∈ 1,2,3,4 . And X is expressed as a tensor, which consists of four variables x ₁, x ₂, x ₃, x ₄, i.e., 1 : X = x 1 , x 2 , x 3 , x 4 T . Besides, T is also expressed as the polynomial coefffcients, whose shape is 4 × 4 × 4. After the tensor contraction operation, the definition of the numerator and denominator of the RPC model can be defined as bellow: f X = T ijk X i X j X k . (5)

Extendedly, the formulation of the RPC model with a set of points is defined as bellow: f b m X = T ijk b X i b m X j b m X k b m , (6) where X b m represents the mth point in the bth batch and T b represents the coefficient tensor in the bth batch. Through element-wise division, the RPC warping of all the points in a batch can be calculated in one shot.

Previous methods usually aggregate the feature volumes to a cost volume by leveraging the cost metric (Hiep et al., 2009). The common practice is to use the variance-based cost metric (CM) to average N − 1 feature volumes. CM considers that the confidence values of the corresponding pixels between the corresponding feature volumes of each view are equally important. The formulation of variance-based cost metric is defined as bellow: C = C M V 1 , … , V 2 = ∑ i = 0 N − 1 V i − V ̄ 2 N , (7) where V ̄ represents the average volume among all feature volumes. However, equal importance is obviously not reasonable, because the satellite images taken by the Ziyuan-3 (ZY-3) satellite have varying shooting camera angles may affect the confidences in the feature volumes due to the matching errors caused by different conditions such as occlusion and non-Lambertian surfaces. If we utilize Eq. 4 to calculate the cost volume, it will affect the final height map estimation.

Therefore, as illustrated in Figure 3, we design an adaptive feature volume aggregation module to calculate an aggregation weighting volume for each feature volume to achieve unequally confidence aggregation. The definition of our module is defined as Eq. 5:

In this way, pixels that may cause the matching errors are suppressed, i.e., the confidences corresponding to pixels are allocated the lower weight, while those with critical feature information are given higher weight. We also formulate our adaptive feature volume aggregation module as follows: C d = 1 N − 1 ∑ i = 1 N − 1 1 + ω c i d ⊙ c i d , (8) where ⊙ denotes Hadamard multiplication and ω(…) is the pixel-wise attention maps adaptively yielded according to per-view cost volumes.

Cost volume regularization (regression to obtain height map) can be seen as a segmentation problem and is handled using the UNet commonly used for semantic segmentation tasks. Therefore, similar to the UNet-shape network used by the previous methods for cost volume regularization, we adopt a similar multi-stage 3D UNet to aggregate neighboring information from a large receptive field, which is composed of three stages (downsample, bottleneck, upsample). First, in the downsampling stage, we leverage the ordinary convolution to obtain the intermediate volume V ̃ ∈ R 64 × D × H 256 × W 256 by three times downsampling. Then we use a bottleneck to learn the high level depth features. We obtain the final regularized cost volume by multiple deconvolutions and skip connections. And the skip connections are used to transfer the corresponding scale intermediate volume V ̃ . The details for the network are shown in Table 1.

In this paper, we adopt the TLC SatMVS dataset proposed by SatMVS. The TLC SatMVS dataset consists of the triple-view images, and the height maps, which are generated by projecting the GT DMS with the corresponding RPC parameters of TLC cameras, which are mounted by Ziyuan-3 (ZY-3) satellite, as shown in Figure 4. And the GT DSMs are obtained from high-accuracy LiDAR observations and ground control point (GCP)-supported photogrammetric software. The dataset consists of 5,011 image patches with resolution 768 × 384.

Figure 5 shows the visualization results of the proposed method on TLC SatMVS Dataset. We restore the depth map from three images. It can be seen that in mountainous areas with large topographic relief, this method has certain effect on depth estimation of multi-view remote sensing images, which verifies the effectiveness of this method.

Table 2 shows the quantitative results on TLC SatMVS Dataset. We compare with traditional and deep learning based MVS methods for satellite images to demonstrate the effectiveness of our model. We have the following observations: 1) We can observe that our method achieves the best among current state-of-the-arts methods in the metrics (MAE:1.597, RMSE:2.036, <2.5m:82.68, Comp:84.32). 2) For traditional MVS method which adopts the pin-hole camera model, e.g., adapted COLMAP, our model outperforms it in all metrics (MAE, RMSE, < 2.5 m , < 7.5 m , Comp, Runtime). 3) Compared with RED-Net, CasMVSNet and UCS-Net, which adopt the pin-hole camera model, although our model has lower scores in < 7.5 m and Runtime, our model achieves SOTA results in other metrics. 4) Furthermore, we also compare the proposed method with the state-of-the-art models in satellite MVS domain, e.g., SatMVS(REDNet), SatMVS(CasMVSNet), SatMVS(UCS-Net). We can observe that our model achieves SOTA results in MAE, RMSE, < 2.5 m , Comp. Besides, we are very close to the current SOTA in terms of < 7.5 m , which exhibits that our method has comparable performance. 5) We can observe that our A-SATMVSNet has a competitive inference time in Table 2. Specifically, our method increases only a slight time in inference Runtime, while extremely outperforms other satellite domain methods in most of metrics, e.g., MAE, RMSE, <2.5m and Comp.

Figure 6 shows the visualization results of different compared methods on TLC SatMVS dataset. It can be seen that the details of the SatMVS(CasMVSNet) method at the corners are seriously missing. However, the results of the method in this paper are more realistic at detail, which are closer to the truth value.

The graph representations of MAE, RMSE, < 2.5 m, and < 7.5 m are visualized in Figure 7, and it is observed that all four metrics exhibit convergence as the epoch increases. This indicates that the model’s performance gradually improves with additional training data. Overall, the visualization of these metrics serves as a useful tool for monitoring and evaluating the performance of the model during training.

Additionally, we also exhibit the visualisation results of our proposed method for four different types of areas: a) darker areas, b) discontinuous areas, c) weakly textured areas, and d) areas with strong height variations, as shown in Figure 8. Our observations are as follows: In a), despite the overall darkness of the scene, our model effectively estimates the height map of the red-boxed area and accurately describes the undulations of the terrain. In b), the red-boxed area is clearly discontinuous with the surrounding area, but our method still produces accurate height estimations without any false noise heights. In c), even though the texture of terrain in the red-boxed area is not very distinct, our method effectively estimates the height map for each pothole. In d), the image contains significant sharp height shifts and some colour noise, but this does not affect the effectiveness of our model in estimating the height map.

We train and validate the large-size satellite images (5,120 × 5,120) with RED-Net, SatMVS(RED-Net) and our A-SATMVSNet on a NVIDIA RTX 3090. The results are shown in Table 3. We can observe that although using the large-size satellite images as input data, which affects the performance of our model, our method still maintain a relatively excellent performance compared with other methods. This confirms the effectiveness and advantage of our proposed model.