MBES refers to a type of highly integrated multibeam bathymetric sensor. It could help with full-coverage depth measurements at high resolution and determine the nature of seafloor surfaces in the deep sea. The basic principle of MBES is shown in Supplementary Figure A1. The MBES transducer is essentially a combination of an acoustic projector array and a perpendicular hydrophone array. The former emits acoustic pulses at a specific frequency, with a narrow opening along-track angle and a wide across-track angle, in a given swath (Mahmud and Yusof, 2006; Costa et al., 2009). The latter is built to listen to echo reflections with received beams. Thus, the seafloor strips, ensonified by the projectors, will intersect with those observed by the hydrophones, producing the beam footprints. When receiving across-track beams of certain time intervals one after another, the position and depth of the seafloor measurement could be calculated, given the angle of incidence and the two-way travel time of each beam (Zhao et al., 2020; Wu et al., 2021). In a complete transmission and reception period, the projector array runs only once to generate acoustic pulses, while the hydrophone array acquires multiple received beams with appropriate delays. As underwater vehicles proceed forward, a strip of water depth measurements reflecting bathymetric mapping at a specific width could be derived from the MBES, providing full coverage of the seafloor surface morphology, which would benefit the identification and location of seafloor categories with high precision, high density, and high efficiency.

Essentially, we utilized the Digital Elevation Model (DEM) to solve the numerical problem of visualizing the geospatial entities of the seafloor surfaces with a finite set of depth measurements from the MBES. The core interpolation calculation allows the production of a gridded multibeam bathymetric map. We applied a weighted average point-to-point interpolation to generate the digital elevation. Assuming that the elevation point is to be inserted in the center of each sliding sampled window, the elevation value E G is determined by approximating the weighted averaging of the surrounding elevation values within the window, which can be formulated as

where the number of neighboring elevation points in the sliding window is denoted as n , E i refers to the i th elevation value, and G i represents the corresponding weight. For the output elevation values, the sum of the products between the surrounding elevation values and their corresponding weights within the window is divided by the sum of all the weights. Each weight G i is defined as the reciprocal of the spatial distance D i between the surrounding elevation points and the center to be inserted,

The greater the spatial distance D i , the smaller the corresponding weight G i , and vice versa. An example of a seafloor strip before and after the interpolation is shown in Figure 1 , where the color bar denotes the water depth values. We can see from the experimental results that the high-resolution multibeam bathymetric mapping could reasonably depict the integrity of the seafloor surface stretching in the DEM, especially the stitching of the gaps in the edges.

The slope refers to the measurement that determines the steepness or degree of inclination in seafloor bathymetric mapping relative to the horizontal plane, which constitutes the fundamental index of benthic habitat and colonization at a variety of scales (Friedman et al., 2013). Multibeam bathymetric mapping can be approximated by a bivariate quadratic equation, and we compute the slope with the first derivative of the elevation values. The slope with the origin at the central point in the local coordinate system within the sliding window is hereby calculated as

The slope direction A could be defined as,

where S x and S y represent the slope with respect to x and y directions, which can take a variety of forms. We determine the slope value of the central point from the finite differential of the surrounding neighbors within the sliding window, as is shown in Supplementary Figure B.1 . The slope S x and S y of the horizontal and vertical directions could be denoted as,

where z 1 − z 4 are the elevation values in the sliding window, respectively, and Δ l is the grid length.

The curvature behaves as a quantitative measurement of the degree of distortion on the surface of the seafloor geomorphic changes, providing a possible assessment of uplift or depression (Shary, 1995). The profile curvature values stand for the stretching morphology of the seafloor surface, with positive curvature attesting to an upwardly concave and a negative curvature, indicating upwardly convex, and a value of zero indicating flat seafloor surfaces. It helps to delimit distinct habitat regions by identifying boundaries in seafloor morphology, delineating between favorable and unfavorable habitats for communities. The curvature is a second spatial derivative of the seabed terrain, which can be expressed as

where l and q are the first derivatives of the elevation values in the horizontal and vertical directions, respectively; r , s , t correspond to the derivative of the horizontal slope with respect to x direction, the derivative of the horizontal slope with respect to y direction, and the derivative of the vertical slope with respect to y direction, respectively.

where z ′ 1 x − z ′ 4 x , z ′ 1 y − z ′ 4 y are the first derivatives of the elevation values in the horizontal and vertical directions within the sliding window, as is calculated in Supplementary Figure C1 .

Surface roughness reflects the degree of the structural complexity of the seafloor surface stretching, which to some extent indicates its macrotopographic characteristics and undulation status, and can be defined as the ratio of the total seafloor surface of the sampled region to a projected plane to decouple measurements from the overall slope (Friedman et al., 2013). Each topographic seafloor surface stretching can be divided into non-overlapping virtual quadrats, and the surface roughness value is derived from each virtual square as

where S s and S p are the seafloor surface area and the horizontal projected area, respectively, in a given virtual quadrat. Let the slope at a given i th elevation point in the sliding window be S i , the corresponding surface roughness R s could then benefit from the calculation of this available topographic factor as follows:

The surface roughness R s of each virtual quadrat with n elevation points can then be expressed as

It was believed that such thematic maps of topological parameters reflecting the seafloor elements and types are effective in classifying seafloor categories in terms of their formation processes and evolution (Burrough and McDonnell, 1998). The topological parameters of high similarity would most likely be shared with the identical seafloor categories. Since most attempts to characterize seafloor elements are limited to a relatively restricted range of morphological attributes, while seafloor types represent characteristic patterns that repeat regardless of scales (MacMillan et al., 2000), we endeavored to utilize high-resolution multibeam bathymetric mapping to extract micro geomorphologic factors such as slope and curvature, as well as macro geomorphologic factors like surface roughness, and to assess the effectiveness of individual or joint morphological cues in distinguishing seafloor surface types. It should be noted that the formation of seafloor surfaces can be viewed from a variety of spatial scales, and the effect of scales involves geomorphology in a complex, hierarchical context. Thus, seafloor classification is related to the issue of scales in different geomorphological settings and the role that morphological cues play in seafloor surface stretching (De Boer, 1992).

The thematic maps of topological parameters for a few example MBES images are shown in Figure 2 , with the original images, the slope, the surface roughness, and the curvature, respectively, displayed from top to bottom. The slope of the seamount generally approached a large value with high-level relief amplitude; the slope of the trench bottom basin was relatively small with nearly flat surfaces; and the slope of the island slope deepwater terrace shifted frequently, representing the divergence of the degree of seafloor surface steepness. The surface roughness provides a macroscopic view of the complexity of seafloor surfaces and reflects the degree to which the seabed terrain is susceptible to erosion. Higher surface roughness values corresponded to more complex or eroded seafloor terrain, e.g., around the island slope deep water terrace. Conversely, flat seafloor surfaces experienced less erosion and exhibited lower roughness values. The curvature directly affected the net erosion, reflecting the degree of seafloor surface fragmentation. When the curvature value of the sea mount was relatively small, the degree of fragmentation was the lowest, and the curvature value of the island slope deep water terrace was relatively large, representing a high degree of fragmentation. The curvature directly affected the net erosion, reflecting the degree of seafloor surface fragmentation. When the curvature value of the sea mount was relatively small, the degree of fragmentation was the lowest, and the curvature value of the island slope deep water terrace was relatively large, representing a high degree of fragmentation.

We initially utilized basic clustering techniques (K-means) to agglomeratively assign elevation points with highly similar topological parameters into the same group and to deviate from the significantly inconsistent outlier elevation points. We could therefore locate and identify individual notions of landforms and geological structures at certain scales with specific physical attributes and translate them to the complete coverage of bathymetric mapping to estimate the potentially appropriate scales as a whole for reference. The individual and joint morphological cues in combinations have served as the input to assess the clustering performances, in terms of PA, MPA, and MIoU, as is shown in Supplementary Table E1 , where the first row is the clustering evaluation of only the bathymetric topographic mapping from MBES, and the second, third, and fourth rows are the evaluation results when introducing the additional morphological cues, respectively, including the slope, surface roughness and curvature. Among them, the clustering performance was superior when both slope and surface roughness were fed as inputs together with the original bathymetric mapping. The comparison of clustering performance with the individual and joint morphological cues for example MBES imaging is shown in Figure 3 , with the original example images, the clustering results from bathymetric mapping+slope, +surface roughness, +curvature, and the ground truth listed from left to right respectively. It was shown that some regions of the trench seamount group were quite easily misclassified as island slopes, leading to many mistakenly divided holes. Owing to the complexity and variability of seafloor surfaces, there exist large divergences even within identical seafloor types and possible similarities across distinct seafloor types, all of which would influence the discrimination process. We have tried to integrate the joint morphological cues into the deep learning-based models to improve the accuracy of distinguishing seafloor categories.

The basic Densely Connected Convolutional Networks (DenseNet) embraces the hypothesis that shorter connections exhibit high performance in a substantially deeper network manner (Huang et al., 2017; Jégou et al., 2017). The feature maps of all previous layers are used as inputs for each layer, and its own feature maps are introduced as inputs to all subsequent layers. Therefore, a basic DenseNet comprising L layers will result in L ( L + 1 ) / 2 direct connections in a feed-forward fashion. Let H i ( · ) be the non-linear transformation implemented in the i th layer, with the output of the i th layer denoted as x i . DenseNet proposes a dense connectivity pattern that introduces direct connections from each layer to all subsequent layers. Consequently, the i th layer receives the feature maps of all previous layers as the input

where [ x 0 , x 1 , … , x i − 1 ] refers to the concatenation of the feature maps produced in the previous layers. For ease of implementation, the multiple inputs of H i (·) could be concatenated into a single tensor. Since the concatenation operation may not be feasible if the size of the feature maps changes during down-sampling, DenseNet would be further divided into multiple dense blocks, with the transition layers between them for convolution and pooling.

We used DenseNet121 as the backbone network of our proposed scheme for seafloor surface classification. The non-linear transformation H i (·) was initially defined as a composite function of consecutive operations, i.e., Batch Normalization (BN), followed by a Rectified Linear Unit (ReLU) and a Convolution (Conv). The design of a 1×1 convolution was introduced as a bottleneck layer before each 3×3 convolution to improve computational efficiency. The DenseNet121 network configuration was made up of four dense blocks. Before entering the first dense block, the initial convolution layer comprised 2 k convolutions of size 7×7 with step size 2, and the number of feature maps in all other layers followed from the setting k . The transition layers took a 1×1 convolution, followed by a 2 × 2 Average pooling between two contiguous dense blocks. At the end of the last dense block, global Average pooling was performed and then a softmax classifier was applied. The number of feature maps in the four dense blocks was 6, 12, 24, and 16, respectively, and the corresponding size of features was 1 / 4 , 1 / 8 , 1 / 16 , 1 / 32 of the original input.

DenseNet121 transforms the input into a feature tensor by gradually reducing the spatial resolution and increasing the number of feature maps along a downsampling path. As for the design and the upsampling path, the Tiramisu model has had great success in the naive extension of DenseNet to fully convolutional networks, while mitigating the linear growth of the feature map explosion in very deep neural networks with very few parameters, replacing the convolution operation with a sequence of dense blocks and the transposed convolution referred to as transition-up (TU) blocks, with an approximately 10-fold reduction with respect to the state-of-the-art models (Jégou et al., 2017). In this paper, in order to explore the possibilities of developing smart capabilities in understanding the seafloor stretching morphology for underwater vehicles, we have updated the DenseNet architecture with an upsampling path of a more simplified transition-up process, i.e., the minimalistic transition-up blocks, which could transform the low-resolution features into high-resolution predictions by recovering details from early layers with blending semantics from deeper layers (Kreso et al., 2017). The design of minimalistic TU blocks is introduced to play the role of the upsampling path in DenseNet121. TU blocks blend the smaller and larger representations whose spatial resolutions differ by a factor of 2 from the upsampling and downsampling paths, respectively, via a skip connection. The blending procedure is repeated recursively by simple summation along the upsampling path, with skip connections arriving from the outputs of each dense block instead of the symmetric encoder-decoder network. The final TU block produces logits at the resolution of the DenseNet stem. The dense predictions at the input resolution are finally obtained by 4× bilinear upsampling. The minimalistic design helps lightweight semantic execution with a low memory footprint and low-dimensional feature tensors during upsampling and discourages overfitting to low-level textures, which potentially presents significant online computation capacities in distinguishing seafloor categories for underwater vehicles.

We adaptively refined the input feature maps along channels by seamlessly integrating the Convolutional Block Attention Module (CBAM) (Woo et al., 2018) into DenseNet121. The CBAM module sequentially infers channel-wise attention maps, which are multiplied by input feature maps. Unlike the Squeeze-and-Excitation (SE) module (Hu et al., 2018), we have tried to exploit the inter-channel relationships by employing both Average pooling and Max pooling in parallel. Given an intermediate feature map X of size H × W × C , with H , W , C being the height, width, and channel number of the feature map, respectively, the spatial dimension of the feature map is squeezed as follows:

where X a v g c and X max c are the outputs of the Average pooling and the Max pooling, respectively, with a size of 1 × 1 × C . The Average pooling aggregates the spatial dimension to suggest the extent of the seafloor surface stretching, and the Max pooling gathers clues of distinctive seafloor surface features to simultaneously infer finer channel-wise attention. Both descriptors allow the global receptive fields to be embedded.

An excitation operation, where the specific activations govern the excitation of the channels by the dependency, feeds the two descriptors into a shared multi-layer perceptron (MLP) with a hidden layer to produce the channel attention map. To reduce the parameter overhead, the hidden activation size is set to C / r , where r is the reduction ratio. The output in MLP is recovered to generate the feature vectors of size 1 × 1 × C . After the shared MLP is applied, the feature vectors are merged by the element-wise summation. In short, channel attention is computed as

where W 1 and W 2 respectively refer to the weights of the two layers, δ stands for the ReLU activation function, and σ denotes the sigmoid function. Finally, the channel attention output M c is multiplied with the initial feature map X to retrieve the newly refined features with calibration,

where ⊗ denotes the element-wise multiplication. The weight coefficient from the channel attention values is broadcast along the spatial dimension during the multiplication to adaptively screen the optimal feature map along the channels.

We have further embraced the idea of a kind of spatial pyramid pooling module (SPP) (He et al., 2015) into our DenseNet121 architecture since it may not sufficiently incorporate the momentous global contextual prior for the receptive fields of the seafloor surface stretching, especially on high-level layers. The basic module of the pyramid scene parsing network (PSPNet) is developed to help exploit and enhance the capability of global context-aware features through aggregation along with sub-regions from multiple receptive fields. We have proposed the introduction of a global context with a sub-region context that enriches to distinguish seafloor surface categories in a pyramidal manner, using both the Average pooling and the Max pooling, as is shown in Supplementary Figure H1 .

Let the number of channels from the channel attention module be C D ; the dimensionality reduction is first performed on the input feature maps by a 1×1 convolution. The Average pooling and the Max pooling simultaneously conclude the feature maps in sub-regions of pyramid scales, with the latter appropriately compensating for the former in detail, and then connect together at pyramid levels along the channel dimension. To maintain the weight of the global seafloor features, a 1×1 convolution layer is applied after each pyramid level. The low-dimensional feature maps are directly upsampled to obtain feature maps of the same size before pooling by bilinear interpolation. Multiple levels of pyramid pooling features are concatenated with the original feature maps before the pooling stage as the final globally enhanced seafloor features, and then output with 1 × 1 × C D / 4 convolution for the next upsampling.

We further proposed fusion strategies to merge with the morphological cues in the context of DenseNet so as to enhance the semantic understanding among seafloor surface types, as is shown in Supplementary Figures I1 and I2 . The first one is that we have attempted to superimpose the morphological features as the input of DenseNet together with the bathymetric seafloor mapping, calibrating the deep-level feature mapping with the help of the channel attention module, enhancing the global feature extraction from the spatial pyramid pooling, and restoring the high-resolution predictions in the up-sampling path for the pixel-level seafloor surface classification. The second strategy is to make an up-sampling of those morphological features through a 1 × 1 convolution as a branch to join with the deep-level feature mapping of the same dimensionality extracted from the DenseNet branch to jointly contribute as the input of the residual block for the subsequent seafloor type prediction. Due to the existence of the identity mapping in ResNet, the residual block could at least copy the previous layer to prevent degradation and simultaneously refine morphological details. In addition, we have evaluated the impact of multiple morphological cues on promoting the descriptiveness and distinguishability of seafloor surface classification.

In our simulation experiment, the developed scheme has been verified by the high-resolution multibeam bathymetric data from the NOAA Office of Ocean Exploration and Research (OER) for the expeditions EX1605L1, EX1605L2, and EX1605L3, with Kongsberg EM302 multibeam echosounders on board the research vessel Okeanos Explorer. The total time of the expedition is 1631.269 h, lasting for 59 days, from the 20th of April to the 10th of July 2016, with a track length of 26703.6897 km and an average speed of 16.33 km/h in the Mariana Trench Marine National Monument and the Commonwealth of the Northern Mariana Islands, as is shown in Supplementary Table K1 .

Meanwhile, the submersible ROV Deep Discoverer (D2), equipped with high-definition cameras and a lighting system, was connected to the camera platform Seirios and the research vessel via an umbilical cable, which provided the possibility of visual cues about the benthic habitat and colonization that are difficult to obtain in the deep sea (Cantwell, 2016). The detailed summary of the ROV Deep Discoverer dive log of EX1605L3 is listed in Supplementary Table L1 , and it includes the latitude and longitude, bottom time, and maximum depth.

First, we essentially utilized the manually labeled seafloor surface annotation as the standard reference so as to identify eight seafloor stretching categories via DenseNet. The normative standard of our manual annotation is listed in Table 1, where the descriptive morphological formation features are commonly known to systematically evaluate the seafloor surface categories (Nishizawa et al., 2009; Harris et al., 2014). In the beginning, we divided the original MBES images into overlapping sub-blocks based on their relatively independent physical attributes of morphological structures at the given scales. We normalized the above MBES images at multiple scales, with their corresponding morphological cues and manual labeling into the basic uniform size 256×256. Such transformed sub-blocks were varied with multiple processing steps, such as random flip, rotation, translation, etc., to promote the diversity of the samples. The selection of the basic uniform size satisfied a comprehensive view of most seafloor topography in our experiment, allowing for interpretation, classification, and validation under the given average swath width of MBES scans. Once a variety of scales with regard to geomorphological formations of seafloor surfaces have been used, normalization would be taken to adapt to the proposed model. In total, 11,720 sub-blocks were chosen, with 8200 samples for training and 3520 for testing, of which 697 samples were originally labeled to the island slope ridge category, 2765 samples to the island slope category, 1145 samples to the island slope deep water terrace category, 2682 samples to the trench seamount group category, 1690 samples to the trench edge slope category, 1240 samples to the trench bottom basin category, 840 samples to the island platform category, and 661 samples to the slope fault basin category. We could further accumulate and refine the seafloor surface annotation as the ground truth through the acquisition of more MBES images.

The configuration of the supercomputing solutions during the model building, training, and testing process was as follows: NVIDIA TITAN Xp graphics card and GeForce GTX 1080Ti graphics cards, an Intel Core i5-2410M CPU with a main frequency of 2.3GHZ, 32GB of memory cards, an Ubuntu 16.04 operating system, a Tensorflow 1.3.0 deep learning framework, a Python3.5 interpreter, data science libraries including Numpy and Pandas, and netCDF data viewers. For optimization, the best Adam optimizer was adopted, among which the exponential decay rate of the first-order moment estimation β 1 and the second-order moment estimation β 2 were 0.9 and 0.99, respectively, by using the cross entropy as the loss function, the learning rate was initially set to 0.001, with the batch size of 16. It should be noted that we examined the hyper-parameters in our simulation experiment, especially the learning rate and the batch size, to ensure the impact on the convergence of our developed model. When the batch size varied from 8 to 32 and the learning rate varied from 0.0005 to 0.01, it was demonstrated from our experimental results that the selected parameters exhibited quite comparable convergence for our proposed scheme.

We employed PA, MPA, and MIoU metrics to quantify semantic segmentation performance with the help of manual annotation. Assuming that there are k categories of seafloor surfaces, let n i j be the total number of image pixels that originally belonged to the i th category but have been incorrectly classified into the j th category, and n j i be the total number of image pixels that originally belonged to the j th category but have been incorrectly classified into the i th category, with n i i the total number of image pixels that belonged to the i th category and have been correctly classified into the i th category.

PA refers to the ratio between the amount of properly classified image pixels and the total number, which can be expressed as the following formula:

MPA refers to the ratio of the number of correctly classified image pixels on a per-category basis, which is then averaged over the total number of categories,

MIoU calculates the average IoU ratio across all categories, which describes the degree of overlap ratio between the intersection and union of categories,

We further started to evaluate the semantic segmentation performance of our proposed scheme. First, we verified the configuration of a variety of backbone networks, such as ResNet50, ResNet101, and DenseNet121, to determine whether it would be more effective to extract the possibly deeper level features for the seafloor surface stretching by means of the identical upsampling modules. As shown in Supplementary Figure J1 , the selection of DenseNet121 initially achieved comparable performance for semantic segmentation of seafloor surface stretching in terms of PA, MPA, and MIoU metrics.

We carried out a series of ablation studies to quantitatively investigate the extent to which the progress of semantic segmentation performance could benefit individually from the improvement of the channel attention module and spatial pyramid pooling in our proposed model. The performance verification for each step is listed in Tables 2 and 3 in our ablation studies, respectively, in terms of PA, MPA, and MIoU metrics. The channel attention module combined both global average pooling and global maximum pooling to optimize the generation of the deep-level feature descriptors. We made the comparative evaluation of the baseline Densenet121, with either the global average pooling or the global maximum pooling, as well as with both types of the pooling. In our experimental results, it has been demonstrated that the effectiveness of both types of pooling behaved better in parallel, where the maximum pooling supplied the possible losses derived from the Average pooling. Spatial pyramid pooling was added to the baseline Densenet121 with the channel attention module, using various pooling selections at multiple pyramid scales. It was shown that the Average pooling alone outperformed the maximum pooling alone, while the two complementary poolings in parallel improved the semantic segmentation accuracy more.

We also examined which types of morphological cues are more relevant to the semantic segmentation of seafloor surface categories, together with the features retrieved directly from MBES imagery via Densenet. Table 4 lists the evaluation of the semantic segmentation accuracy by merging multiple morphological cues into the DenseNet backbone network with the embedded channel attention module (C) and spatial pyramid pooling module (S) in the context of two types of feature fusion strategies, including slope, roughness, curvature, slope + roughness, slope + curvature, roughness + curvature, and slope + roughness + curvature. The first mode concatenated the individual or joint morphological cues with multibeam bathymetric seafloor mapping in advance to generate the multi-channel input for DenseNet, and the resulting fused feature maps would be adaptively optimized with the channel attention module, advanced into global feature representation with the spatial pyramid pooling, then restored to high-resolution predictions from up-sampling with the aid of transition-up blocks, outputting the subsequent semantic seafloor classification. In the second mode, the bathymetric seafloor mapping was individually input into DenseNet, with the channel attention module and spatial pyramid pooling employed. Concatenated feature maps of the same dimensionality were extracted from up-sampled morphological cues by convolution from another branch in parallel, and then commonly fed the feature fusion into the residual block to output the seafloor type prediction. From our experimental results, the second mode achieved the overall performance improvement compared to the baseline and the first mode, which to a certain extent plays a role in compensating for the loss in down-sampling, thereby improving the descriptiveness and distinguishability of seafloor surface categories. The morphological cues of slope + roughness exhibited better performance, while the curvature did not show a significant improvement in accuracy. It was inferred that the slope tends to indicate the degree of steepness in seafloor surface stretching, and the surface roughness might display the extent of erosion in the seafloor surface topography, all of which contribute to the semantic segmentation. Also, the curvature reflects the degree of fragmentation, which might not be seen as a very distinguishable index and might lead to misclassification to a large extent.

Furthermore, the semantic segmentation accuracy of each individual seafloor surface category was systematically evaluated against the classic Fully Convolutional Network (FCN) (Long et al., 2015) in terms of the IoU measure, as is shown in Table 5 . Since IoU describes the degree of overlap between the actual outputs and the expected outputs of each category prediction, the higher the IoU index value, the better the segmentation performance. It could be seen from our experimental results that our proposed scheme has made general progress in distinguishing individual seafloor surface categories from the rest via DenseNet with the channel attention module and spatial pyramid pooling strategies, which tended to enhance the feature mapping capability of the entire network and thus improve the segmentation accuracy accordingly.

The overall performance evaluation of our developed scheme has been illustrated in Table 6 , where the semantic segmentation accuracy calculation, from the baseline Densenet121 with the minimalistic transition-up (TU) blocks (DT), embedded with the channel attention module (DTC) and spatial pyramid pooling module (DTCS), to the coupled feature fusion with the morphological cues (DTCSF), was quantitatively measured step by step in terms of PA, MPA, and MIoU metrics.

We further performed the comparative evaluation with some state-of-the-art models for semantic segmentation of seafloor surface stretching, including FCN-8s, SegNet (Badrinarayanan et al., 2017), RefineNet (Lin et al., 2017), PSPNet (Zhao et al., 2017), DeepLab v3+ (Chen et al., 2018), and our developed scheme, as shown in Table 7 . The classic FCN-8s network integrates the multi-layer feature maps during down-sampling in FCN. The SegNet network calls the pooling index at the corresponding encoder in the decoder to upsample the feature map through the unpooling operation. RefineNet explicitly exploits all the information available along the down-sampling process to enable high-resolution prediction through long-range residual connections. PSPNet captures global context through different-region-based context aggregation by the pyramid pooling module to improve network performance. DeepLab v3+ makes use of an encoder-decoder to perform multi-scale information fusion while retaining the dilated convolutions and Atrous Spatial Pyramid Pooling (ASPP) layer of the original DeepLab series. It should be noted that the above segmentation results for seafloor stretching categories were initially generated from the average product of 5-fold cross-validation with our developed model by dividing into mutually exclusive subsets with nearly equal numbers of randomly selected samples. From the experimental results, it was demonstrated that our proposed scheme had achieved a significant improvement in semantic segmentation performance, with PA, MPA, and MIoU metrics reaching up to 89.87%, 82.01%, and 73.52%, respectively. The model also exhibited a high level of stability in terms of PA, MPA, and MIoU metrics with a series of cross-validation rounds.

The semantic segmentation of multibeam bathymetric seafloor mapping has been further visualized, as is shown in Figure 5 , where the example MBES image, the segmentation results of both FCN-8s and our proposed scheme, and the corresponding labels are listed from left to right, respectively, with the island slope ridge in red, the island slope in green, the island slope deepwater terrace in cyan, the trench seamount group in yellow, the trench edge slope in blue, the trench bottom basin in purple, the island platform in orange, and the slope fault basin in black. From our experimental results, it was demonstrated that our developed scheme visually outperformed the classic FCN-8s, and simultaneously enhanced the details in-between edges, with the ability to preserve the salient features and eliminate redundancy on a global scale, showing its superiority in the descriptiveness and distinguishability of the seafloor surface categories. Some semantic segmentation results of the example multibeam bathymetric seafloor mapping along waypoints of the expedition track around the Mariana Trench are shown in Figure 6 , where the location of the waypoint, the original example MBES images, the segmentation results, and the ground truth are listed from left to right.

Furthermore, we made an attempt to focus on a more delicate observation and a preliminary study of how the seafloor surface stretching functions as a submarine benthic habitat and what type of biogeographic pattern distribution of the benthic organisms are present in the extremely deep sea, with the help of both the acoustic sensor on board the research vessel Okeanos Explorer and the optical sensor mounted in the ROV Deep Discoverer. A total of 10,000 underwater images of the dominant resident biological species and their corresponding habitats, at dive depths within a range of 250-5000m, included Rimicaris, Austinograea, Symphurus thermophiles, Bathymodiolus, Phenacolepadidae, Shinkailepas, Thoridae Lebbeus, Lamellibrachia, etc., and were considered as our alternative underwater vision dataset for this preliminary study. The primary benthic species retrieved from the video of each dive by the ROV Deep Discoverer during the EX1605L3 expedition route are recorded in Supplementary Table M1 .

We established a global geographical link between the ROV dive path and the MBES bathymetric mapping route. Figure 7 shows the connection between the latitude and longitude of the example ROV dive paths and the location of the MBES imaging survey, including the original example MBES bathymetric mapping, the corresponding seafloor surface categories, the ROV dive paths on Eifuku Seamount and Daikoku Seamount, and the possible typical benthic habitats retrieved from visual cues along the paths, thus linking the seafloor surface topography with the primary benthic biogeographic patterns. We tried to statistically match the corresponding primary benthic habitats and species with the seafloor surface stretching by roughly retrieving the microgeographic cues from each ROV dive, and subdividing the benthic habitats with the visual cues from optical sensing. Figure 8 lists some examples of dominant benthic habitats and species that visually reflect the possible biogeographic patterns that respectively appeared and were distributed at distinct locations of the seamount above Figures 8A–F and the seamount below Figures 8G–L , which also makes it possible to provide an initial insight into the diversity and distribution of the benthic community.

It could be seen from the visual clues from the ROV dive on Eifuku Seamount that the benthic species, especially fish and octocoral fauna, were unexpectedly diverse, and the typical geomorphology discovered included the crater wall and the hydrothermal chimney structure near the summit, while the visual clues from the ROV dive on Daikoku Seamount demonstrated the high activity of the hydrothermal vents, the possible evidence of the recent eruption, the volcaniclastics, the sulfur pond and the thick volcanic smoke, the plume, and the flatfish communities, e.g., Symphurus thermophilus and Gandalfus yunohana. The extent of the seafloor surface stretching and the estimation of the primary benthic biogeographic patterns reflect the coupling variation of multivariate environmental variables in the deep sea. The associative study derived from the sparse observation statistics through both acoustic and optical sensors not only produces the possibility of capturing the potential relationships between the full coverage of seafloor mapping and the benthic habitats, even the benthic species assemblage maps, but also provides the opportunity to examine the predicted biogeographic patterns with better-described variations and uncertainties towards the distinct geographical characteristics of seafloor surfaces.