# Comparison of automated crater catalogs for Mars from
benedixDerivingSurfaceAges2020 and leeAutomatedCraterDetection2021

###### Abstract

Crater mapping using neural networks and other automated methods has increased
recently with automated Crater Detection Algorithms (CDAs) applied to
planetary bodies throughout the solar system. A recent publication by
benedixDerivingSurfaceAges2020 showed high performance at small scales
compared to similar automated CDAs but with a net positive diameter bias in
many crater candidates. I compare the publicly available catalogs from
benedixDerivingSurfaceAges2020 and leeAutomatedCraterDetection2021 and show
that the reported performance is sensitive to the metrics used to test the
catalogs. I show how the more permissive comparison methods indicate a higher
CDA performance by allowing worse candidate craters to match ground-truth
craters. I show that the benedixDerivingSurfaceAges2020 catalog has a
substantial performance loss with increasing latitude and identify an image
projection issue that might cause this loss. Finally, I suggest future
applications of neural networks in generating large scientific datasets be
validated using secondary networks with independent data sources or training
methods.

Earth and Space Science University of Toronto 60 St. George St., Toronto,
Ontario, Canada Christopher<EMAIL_ADDRESS>
The apparent performance of automatic crater detection algorithms is sensitive
to the choice of metric used to gauge the performance.

Feature Detection Algorithms trained to find spatial features should be used
with appropriately projected images where possible.

As neural networks extend crater detection beyond practical human limits,
validating their results with independent networks is critical.

## Plain Language Summary

I have compared two computer programs that use neural networks to identify
craters on solar system bodies. The crater catalogs created by these programs
have been previously compared against an existing human-made catalog to
measure their performance, but each comparison used a different method making
it challenging to understand the difference between the catalogs. In this
paper, I use the comparison methods from these two independent papers, discuss
where the catalogs agree and disagree, and suggest why they disagree. I
emphasize the need for more accurate and consistent methods to validate
automatic crater-mapping programs, especially when it is impractical for
humans to check the craters found by automated methods.

##  1 Introduction

Crater mapping using neural networks and other automated methods has increased
in recent years with automated Crater Detection Algorithms (CDAs) applied to
planetary bodies throughout the solar system including Mercury and the Moon
[<]e.g.,¿[]Silburt2019 Mars and Pluto
[<]e.g.,¿[]leeAutomatedCraterDetection2021,Lee2018a,benedixDerivingSurfaceAges2020
and large asteroids [<]e.g.,¿[]latorreTransferLearningRealtime2023. These
methods approach or exceed the performance of human experts at identifying
craters from imagery and topography data and have recently produced catalogs
that are too large to be validated by humans in a reasonable time [Lagain,
Benedix . (2021)]. In this work, I compare the results from two automatically
generated catalogs for Mars and two additional catalogs described below. I
calculate the performance metrics against a shared ground-truth and show how
the different test criteria originally developed for these catalogs [Lee
(2019), Benedix . (2020)] leads to different apparent quality of each catalog.

In the next section, I describe the catalogs and discuss some of the stated
limitations of the neural networks and catalogs provided. I then reproduce the
comparison algorithms from leeAutomatedCraterDetection2019 and
benedixDerivingSurfaceAges2020 and use both methods to calculate the
performance of each crater catalog in this comparison. I compare the
performances of these catalogs as a function of crater location and size and
discuss the differences in performance. I also identify a possible calculation
error in benedixDerivingSurfaceAges2020 and suggest a correction. Finally, I
give my conclusions on this comparison and suggest ways to mitigate errors and
improve future catalogs.

##  2 Methods

###  2.1 Catalogs

In this work, I compare two publicly available crater catalogs generated by
neural networks, Benedix from benedixDerivingSurfaceAges2020 and DeepMars2
from leeAutomatedCraterDetection2021. I also compare a catalog generated as
part of leeAutomatedCraterDetection2021 (irH1K) but not separately published,
and a new catalog based on benedixDerivingSurfaceAges2020 (B_Rescaled) that I
created as part of this work to correct a proposed error in the original
catalog.

The Benedix catalog was generated by benedixDerivingSurfaceAges2020 from
infrared images of the planet using the “You Only Look Once” [Redmon . (2016),
Redmon Farhadi (2018)] neural network, reporting 328,833 detected craters at
diameters between 1km and 200km, covering the surface between $\pm$65 degrees
latitude. Following [Benedix . (2020)], the trained network was used to
identify 90 million small craters [Lagain, Benedix . (2021)], and the same
authors have sought to validate the Robbins2012a catalog using computer-aided
tools [Lagain, Bouley . (2021)]. The Benedix catalog was downloaded for this
work from the supplementary material of benedixDerivingSurfaceAges2020.

The DeepMars2 catalog was generated by leeAutomatedCraterDetection2021 using
digital terrain models (DTMs) of the surface and co-registered infra-red
imagery using a ResUNET network [Diakogiannis . (2019)]. Two separate networks
were trained using the same Robbins2012a catalog as the ground-truth with
randomly selected swatches from across the planet that were re-projected onto
an orthographic projection. In addition to the combined catalog that used DTMs
and IR detections, a high-resolution IR catalog was also developed in
leeAutomatedCraterDetection2021 but truncated to 3.5km diameter and above to
match the DTM limitations. In the following sections, we relax this constraint
and include this infra-red imagery catalog truncated to a minimum of 1km in
diameter (labeled as irH1k). The DeepMars2 catalog was downloaded for this
work from the associated BorealisData repository for the
leeAutomatedCraterDetection2021 paper at leeCraterCatalogsSoftware2020. The
irH1k catalog was created as part of leeAutomatedCraterDetection2021 but is
made available on the BorealisData repository associated with this work at
leeReplicationDataComparison2023.

Finally, as part of this work, a possible error was identified in the Benedix
catalog, and a correction was applied to create the B_Rescaled catalog. The
correction involved rescaling the calculation of the pixel-to-geophysical
size, resulting in smaller crater sizes at high latitudes in the new catalog.
The B_Rescaled catalog is available from the BorealisData repository
associated with this work at leeReplicationDataComparison2023.

Newer models have been developed since benedixDerivingSurfaceAges2020 and
leeAutomatedCraterDetection2021, but typically examine a small area of Mars as
validation [Yang Cai (2022), Giannakis . (2023)] and provide no publicly
accessible crater catalogs. There have also been improvements in surface
imagery since Robbins2012a with improved geo-location of images [Robbins
(2017)] and higher resolution visible imagery (http://murray-
lab.caltech.edu/CTX/) from the Mars Reconnaissance Orbiter Context Camera
[Malin . (2007)].

In the rest of this paper, I refer to the algorithms and catalogs from
leeAutomatedCraterDetection2021, benedixDerivingSurfaceAges2020, and
Robbins2012a. To avoid misattribution, I refer to the catalogs by their names
above (Benedix, DeepMars2, irH1K, B_Rescaled, and Robbins) in the text and
figures, and use the acronyms L19 [Lee (2019)] and B20 [Benedix . (2020)] to
refer to the comparison methodology developed by those authors.

###  2.2 Ground Truth Comparison methodology

Both leeAutomatedCraterDetection2021 and benedixDerivingSurfaceAges2020
compare their results against the Robbins catalog as a source of ‘ground-
truth’. In each paper, the comparison required matching every crater from the
CDA catalog with a crater in the Robbins catalog and required the definition
of a matching algorithm. The L19 algorithm uses the following criteria to
determine a positive crater match to the ground-truth:

$\displaystyle F_{\mathrm{D}}$
$\displaystyle=\frac{D_{\mathrm{C}}-D_{\mathrm{G}}}{min(D_{\mathrm{C}},D_{\mathrm{G}})}$
$\displaystyle\leq 0.25,$ $\displaystyle F_{\mathrm{Y}}$
$\displaystyle=\kappa_{Y}\frac{Y_{\mathrm{C}}-Y_{\mathrm{G}}}{min(D_{\mathrm{C}},D_{\mathrm{G}})}$
$\displaystyle\leq 0.25,$ $\displaystyle F_{\mathrm{X}}$
$\displaystyle=\kappa_{X}\frac{X_{\mathrm{C}}-X_{\mathrm{G}}}{min(D_{\mathrm{C}},D_{\mathrm{G}})}$
$\displaystyle\leq 0.25.$ (1)

Where $D$,$X$,$Y$ are the diameter (km), longitude (degrees), and latitude
(degrees) for the catalog crater($C$) and ground-truth crater ($G$).
$\kappa_{Y}=\frac{2\pi\,R_{\mathrm{mars}}}{360}$ and
$\kappa_{X}=\kappa_{Y}\cos{Y_{R}}$ convert from distance in degrees to
distance in kilometers for a spherical planet.

A positive match is identified in L19 only when all three criteria
($F_{\mathrm{D}}$,$F_{\mathrm{Y}}$,$F_{\mathrm{X}}$) are below 0.25. That is,
when a candidate crater and ground-truth crater are both within 25% of the
smaller crater’s diameter, and the centers of the two craters are within 25%
of the smaller crater’s diameter.

Based on benedixDerivingSurfaceAges2020, the B20 method uses the following
criteria to determine a positive crater match:

$\displaystyle F_{\mathrm{D}}$
$\displaystyle=\frac{D_{\mathrm{C}}-D_{\mathrm{G}}}{min(D_{\mathrm{C}},D_{\mathrm{G}})}$
$\displaystyle\leq 0.5,$ $\displaystyle F_{\mathrm{Y}}$
$\displaystyle=\frac{Y_{\mathrm{C}}-Y_{\mathrm{G}}}{|Y_{\mathrm{G}}|}$
$\displaystyle\leq 0.02,$ $\displaystyle F_{\mathrm{X}}$
$\displaystyle=\frac{X_{\mathrm{C}}-X_{\mathrm{G}}}{|X_{\mathrm{G}}|}$
$\displaystyle\leq 0.02.$ (2)

A positive match is identified in B20 if the diameters are within 50% and the
crater centers are within 2% of their absolute location in longitude or
latitude. Note that the B20 method doesn’t compare location relative to the
diameter as L19 does, using the location in degrees in the denominator
instead. B20 also allows matches when the craters do not overlap, while L19
requires an overlap by construction (the center of one crater must be within
the circle representing the matching crater).

###  2.3 Comparison Metrics

I calculated the performance of the catalogs with respect to the ground-truth
catalogs using the standard metric calculations
[<]e.g.,¿[]leeAutomatedCraterDetection2019,Goodfellow2016 that are also used
in leeAutomatedCraterDetection2021 and benedixDerivingSurfaceAges2020.

Each candidate catalog is constructed from a number $N_{C}$ total craters,
compared to $N_{G}$ ground-truth craters. Comparing these catalogs identifies
a number $N_{T}$ matching craters (‘true positive’) and $N_{F}$ non-matching
craters (‘false positive’). The performance metrics are all calculated using
these crater counts. Recall is the fraction of craters in the ground-truth
catalog that also exists in the candidate catalog: $R=N_{T}/N_{G}$. Precision
is the fraction of craters in the candidate catalog that exists in the ground-
truth catalog $R=N_{T}/N_{C}$. Since these metrics measure two distinct
properties of the catalog, a blended metric ($F_{1}$ score) is also typically
used, defined as the harmonic mean of Precision and Recall: $F_{1}=2PR/(P+R)$.

All three metrics are typically expressed as percentages. For example,
benedixDerivingSurfaceAges2020 reported their catalog has a Recall of 77% and
Precision of 75% (their table 1), leeAutomatedCraterDetection2021 reported
their catalog has a Recall of 80% and Precision of 87% (their table 1).
However, since these catalogs have different bounds on diameter and different
spatial coverage, a direct comparison of these numbers is not practical.

For the CDA craters with a matching ground-truth crater (i.e., craters part of
the $N_{T}$ count above), I also calculated the percentage errors between
matching pairs of craters based on equations 2 and 1 above. The percentage
errors are then aggregated and plotted as histograms following figure 3 of
benedixDerivingSurfaceAges2020. Since the L19 method calculated errors
relative to diameter, these errors cannot be directly compared with the B20
method, but the diameter errors used the same calculation and can be compared.

For each matching crater, I also calculated the Intersection-Over-Union (IOU)
metric that measures the overlap (intersection) of two features relative to
the total exposed surface area (union)
[<]e.g.,¿[]fairweatherAutomaticMappingSmall2022,benedixDerivingSurfaceAges2020.
The IOU calculation takes the pair of matching craters and calculates the
ratio

$IOU=\frac{\textrm{Total intersecting area of the crater pair}}{\textrm{Total
exposed area of the crater pair}}$

Perfectly overlapping craters have a maximum value of $IOU=1$, while non-
overlapping craters have a value of $IOU=0$. Values between 0 and 1 include
mismatched sizes and locations but require some overlap. For each of the
matched crater pairs, I calculated the IOU metric assuming circular craters
since all of the catalogs used in this comparison report crater parameters
assuming circular features, even where the underlying craters are degraded or
elliptical.

Finally, I calculated the Recall and Precision as a function of location and
size to show the performance of each network across the planet and as a
function of crater diameter. In principle, an unbiased network should perform
equally well across the planet but will perform worse as the crater size
approaches the native image scale and crater features become harder to detect.

Note that the main text of benedixDerivingSurfaceAges2020 also describes an
‘uncorrelated crater’ that is calculated in their table 1 as the number of
ground-truth craters not found by the CDA (usually called a ‘false negative’),
but it is described as indicating “where the algorithm found a crater that is
not present in the [ground-truth].”. Similarly, their ‘false positive’ craters
are described as being found by the ground-truth but not by the CDA (their
section 3.3). These definitions are reversed compared to normal usage, and in
their [Benedix . (2020)] supplementary material S3 and their table 1 the
calculations follow the standard definitions given above.

##  3 Results

Four catalogs are used in this section. The DeepMars2 catalog from
leeAutomatedCraterDetection2021, the irH1k infra-red channel catalog from
leeAutomatedCraterDetection2021, the Benedix catalog from
benedixDerivingSurfaceAges2020, and my attempted correction to the Benedix
catalog in the B_Rescaled catalog. Table 1 summarizes the performance of each
catalog measured against the Robbins ground-truth catalog.

The table shows the performance for both the B20 and L19 methodology for
craters from 1.5km to 10km in diameter within 65 degrees from the equator (the
limits used in benedixDerivingSurfaceAges2020). The values for the Benedix
catalog differ by 2-3% compared to benedixDerivingSurfaceAges2020, but these
numbers are sensitive to the numerical algorithm used to match craters. For
example, I could increase the Recall and Precision by 2% by allowing each
ground-truth crater to match multiple candidate craters, but it’s not clear
that benedixDerivingSurfaceAges2020 allowed this to occur.

The performance values for the DeepMars2 are consistent with
leeAutomatedCraterDetection2021, with a 1% difference in both Recall and
Precision caused by not removing ${\sim}100$ overlapping craters in the
ground-truth catalog for this work.

The B20 and L19 crater matching methods produce very different results. The
effect is most noticeable in the Benedix catalog, where performance drops by
35 percentage points using the stricter L19 method (a factor of two in CDA
performance), but irH1k also drops by 10% when changing from the B20 method to
the stricter L19 method. In contrast, the DeepMars2 is largely unaffected
because most non-overlapping craters were removed during post-processing [Lee
Hogan (2021)].

Method | Test | R | P | F1 | TP | FP | FN | test count
---|---|---|---|---|---|---|---|---
B20 | DeepMars2 | 20 | 95 | 33 | 32,995 | 1,672 | 130,416 | 34,667
Benedix | 72 | 70 | 71 | 117,485 | 50,875 | 45,926 | 168,360
B_Rescaled | 72 | 73 | 72 | 117,239 | 43,033 | 46,172 | 160,272
irH1k | 76 | 74 | 75 | 124,190 | 44,215 | 39,221 | 168,405
L19 | DeepMars2 | 20 | 92 | 32 | 32,013 | 2,654 | 131,398 | 34,667
Benedix | 36 | 35 | 35 | 58,777 | 109,583 | 104,634 | 168,360
B_Rescaled | 47 | 48 | 47 | 76,410 | 83,862 | 87,001 | 160,272
irH1k | 65 | 63 | 64 | 105,717 | 62,688 | 57,694 | 168,405
|  |  |  |  |  |  |  |
Table 1:  Table of metrics for each catalog with respect to the Robbins2012a
catalog for craters within 65 degrees of the equator and between 1.5km and
10km in diameter. Recall, Precision, and F1 are rounded percentages, while all
other numeric columns are integer values. A total of 163,411 ground-truth
craters were in the comparison. R=Recall, P=Precision, F1=harmonic mean score,
TP=true positive count,FP=false positive count, FN=false negative count (maybe
listed as uncorrelated craters in benedixDerivingSurfaceAges2020),

Figure 1 shows the size frequency distribution using the same bins as
benedixDerivingSurfaceAges2020 (20 diameter bins per log-decade). The plot
shows that most catalogs have similar distribution shapes, even with different
precision and recall metrics. The Salamuniccar2012 catalog is truncated above
25km, and it shows discretization artifacts at lower diameters, as does
DeepMars1 from leeAutomatedCraterDetection2019. The irH1k catalog only
processed images up to 1.3 degrees across (a maximum crater size of about
40km).

Figure 1: Size Frequency Distribution of all craters in each catalog. Benedix
from benedixDerivingSurfaceAges2020, DeepMars1 from
leeAutomatedCraterDetection2019 , DeepMars2 from
leeAutomatedCraterDetection2021, irH1k from the high resolution IR channel
[Lee Hogan (2021)], Robbins from Robbins2012a, and Salamuniccar from
Salamuniccar2012.

Figure 2 shows the distribution of errors for craters matched between a
candidate catalog and the Robbins catalog using the B20 method (equation 2).
The shapes of these figures agree reasonably well with
benedixDerivingSurfaceAges2020. The location error is relatively unbiased with
a prominent central peak, and the diameter errors shown here match the
reported positive bias[Benedix . (2020)]. The DeepMars2 catalog has a lower
total count but similar biases to the irH1K catalog. The diameter biases have
different signs between the benedixDerivingSurfaceAges2020 and
leeAutomatedCraterDetection2021 catalogs, possibly due to a difference in
network training.

Figure 2: Distribution of errors in longitude (left), latitude (center), and
crater Diameter (right) calculated according to the B20 method (equation 2)
using 500 bins in each histogram. Calculations were made relative to the
ground-truth catalog such that positive numbers imply a positive bias in the
CDA catalog.

Figure 3 shows the calculation of errors using the L19 method (equation 1).
Note that all errors are now calculated relative to the diameter of the
smaller crater in a matched pair since the L19 method included that diameter
as the denominator in each calculation. The shape of the diameter plot is
essentially the core of the same subplot in figure 2. In contrast, the
longitude and latitude figures are functionally different, showing the
position error relative to the crater diameter rather than in absolute
location as in figure 2.

Figure 3: Distribution of errors in longitude (left), latitude (center), and
crater Diameter (right) calculated according to the L19 method (equation 1).
The horizontal extent of the longitude and latitude plot appear larger here
but are expressed relative to the crater diameter rather than the absolute
value of longitude or latitude used in figure 2. Calculations were made
relative to the ground-truth catalog such that positive numbers imply a
positive bias in the CDA catalog

Figure 4 shows the distribution of Intersection-over-Union (IOU) values for
matching craters using both methods (L19 in solid lines, B20 in dashed lines)
and for each of catalogs relative to Robbins. We did not calculate IOU metrics
in training the L19 or L21 CDA because the algorithm matched crater rims with
narrow widths. benedixDerivingSurfaceAges2020 reported that validation and
training were based on IOU comparisons with the ground truth and required a
high threshold for a positive detection. The plot shows $IOU>0.3$ for all
matches in the L19 method (solid lines) compared to the theoretical minimum of
$IOU\approx 0.25$ given the constraints in equation 1. Since the B20 method
counts some non-overlapping craters as matches, there is a population of
matching craters with no intersection and a value $IOU=0$, marked with a
horizontal line for each catalog at the left edge of the plot.

Figure 4: Intersection-over-Union values for the four crater catalogs. Solid
lines show the IOUs of crater pairs matched using the L19 method. Dashed lines
show the same catalogs matched using the B20 method, and horizontal lines at
IOU=0 highlight the population of non-overlapping craters ($IOU=0$) found
using the B20 method. The histograms used 100 bins.

Figure 5 shows the Recall and Precision metrics as a function of longitude,
latitude, and diameter based on matched-craters using the L19 method. The
DeepMars2 catalog has the largest difference between Recall and Precision
because it removes all craters below 3.5km in diameter and thus misses 100,000
small craters, lowering the Recall across all longitudes and latitudes. All
four catalogs have consistent performance in longitude, fluctuating about 10%
overall. In latitude, the Benedix catalog has a significant performance drop
with increasing latitude in both hemispheres, while the irH1k has a small drop
in performance with latitude with worse performance in the northern
hemisphere. All four catalogs lose performance as the crater diameter
approaches the pixel scale of the source imagery, as expected. For DeepMars2,
this drop-off occurs at 3.5km. For the other (IR) catalogs, the drop-off
occurs near 2km.

Figure 5: Recall (solid) and Precision (dashed) plotted as a function of
longitude, latitude, and diameter for catalogs matched using the L19 method.
Longitude and latitude are aggregated into 5-degree bins, and diameter is
aggregated into 20 bins per log-decade.

##  4 Discussion

To generate the comparisons, I needed to make several assumptions to reproduce
the B20 method and crater-matching catalogs. Neither Benedix nor DeepMars2
include the corresponding craters from Robbins. Instead, I built the crater-
matching algorithms following the description in
benedixDerivingSurfaceAges2020 (B20, equation 2 above) and
leeAutomatedCraterDetection2019 (L19, equation 1 above).

Using these two comparison methods, I identified craters that matched with a
unique ground-truth crater and used these matches to generate the statistics
presented here. The results in table 1 deviate by a few percent compared to
benedixDerivingSurfaceAges2020, but these numbers can be adjusted by a percent
or two in either direction by changing the interpretation of parameters in the
matching algorithm. Still, the metrics are comparable within this work because
they are all calculated with the same methodology and code. The SFD (figure 1)
and the error distribution (figure 2) match well with the corresponding
figures from benedixDerivingSurfaceAges2020, providing some confidence that
the implementation of the B20 method matches the original. Similarly, the
metrics and biases reported in leeAutomatedCraterDetection2021 agree with the
values calculated here using the L19 method.

The B20 algorithm produced higher apparent performance because it included
craters that don’t overlap with their targets, accounting for 50% of Benedix
matches, 3% of DeepMars2, and 16% of irH1k. The L19 algorithm required
overlapping craters for a valid match and thus produced higher IOU
distributions in figure 4 but lower overall performance. Craters that matched
according to L19 also matched according to B20 except for the few craters
along the equator or prime meridian where 2% of absolute longitude is smaller
than 25% of the crater diameter.

The performance of each catalog is reasonably consistent with longitude and
has the expected drop as the crater diameter approaches the image pixel scale.
There is a drop in performance with latitude in the catalogs, but a much
larger drop in performance with latitude in Benedix suggests a more
fundamental error caused by neglecting to project images onto orthographic
projections before processing.

The source of the error is not obvious from the available catalogs without
access to the processing algorithms. However, of the 58,578 Benedix catalog
craters matched by B20 but not by L19, only 13,383 would satisfy the diameter
constraints of the L19 method. If the crater diameters in Benedix are affected
by a projection error, or equivalently using the wrong pixel-to-kilometer
scaling factor, they can be corrected in the existing dataset because the
catalog includes pixel scale parameters and location information. Scaling the
diameters by cosine latitude improves the catalog performance using the L19
method by allowing an additional 23,309 craters to match the ground-truth
catalog.

This hypothesis is supported by the raw pixel size information provided in the
Benedix catalog. Figure 6(top) shows the density map of craters as a function
of the latitude and the ratio of the height (‘dy’) to width (‘dx’) of the
rectangle encompassing the crater in pixel values. This ratio should be near 1
for orthographically projected images for circular craters at all latitudes,
but has a cosine latitude dependence for Plate Carèe (unprojected, or
latitude-longitude) images. In unprojected images, circular craters appear
stretched in longitude at higher latitudes, and the ratio dy/dx decreases for
circular features, as shown in the figure. This issue suggests the image
swatches were presented to the CDA as unprojected images where circular
physical features appear elliptical in the images, and the CDA was trained to
find horizontally stretched ellipses rather than circles.

The pixel level information can also be used to determine how the conversion
to geophysical units (in this case, kilometers) was done by
benedixDerivingSurfaceAges2020. Three image scales were used by
benedixDerivingSurfaceAges2020, and these scales accurately map all of the
report sizes in the ‘x’ dimension from the pixel size to the reported physical
size of the crater in kilometers. However, in these images, ‘x’ corresponds to
longitude and should have a different scaling with latitude, whereas the ‘y’
axis corresponds to the pixel scale in latitude and does not vary with
location. Thus, the crater’s pixel size in the ‘y’ axis should have been used
instead of the size in ‘x’. As a result, the Benedix catalog has an increasing
positive bias in crater diameter as a function of latitude (figure 6(bottom)).

Figure 6: (top) Heat map of crater population as a function of latitude and
the ratio of height-to-width (dy/dx) of the pixel-level rectangle in the
Benedix catalog. Darker shades represent a higher density of craters, and the
black line shows the cosine of latitude. (bottom) Scaled diameter bias as a
function of latitude, calculated as the median $F_{D}$ from equation 1 (solid)
or equation 2 (dashed) relative to its value at the equator. The drop off in
apparent bias at high latitude is because the high biases exceed the threshold
in each method, and the remaining craters have lower bias.

My proposed fix for this error is to apply the image pixel scaling to the ‘dy’
pixel values instead of ‘dx’, which produces the B_Rescaled catalog shown in
each table and figure. This new catalog has improved global metrics, reduced
errors, and increased IOU values relative to Benedix. The Recall and Precision
metrics also remain higher with latitude and increase in all dimensions
compared to the original catalog. However, this improvement does not
definitively prove that the Benedix catalog has a projection issue, and other
issues might cause the observed bias.

##  5 Conclusions

I have examined crater catalogs generated with neural networks [Lee Hogan
(2021), Benedix . (2020)] whose results suggested similar performance but were
calculated using different methodologies. Using the comparison method from
leeAutomatedCraterDetection2019 I found that the
benedixDerivingSurfaceAges2020 catalog has significantly lower Precision
compared to leeAutomatedCraterDetection2021 generating 30 times more false
crater detections, but did contain more craters in total and with a higher
Recall. When compared to the infrared channel of
leeAutomatedCraterDetection2021 (the irH1k catalog) that found craters down to
1km, the benedixDerivingSurfaceAges2020 catalog has lower performance in all
metrics.

I also identified a possible scaling error in the
benedixDerivingSurfaceAges2020 catalog related either to image projection or
using the incorrect pixel-to-kilometer scaling constant. I proposed a
correction that rescaled to crater sizes by the correct scaling coefficient
and generated the B_Rescaled catalog. This rescaled catalog improved in all
metrics over the original catalog.

Since the validation of crater catalogs becomes increasingly difficult with
their increasing size, the ground-truth comparison methods must be clearly
defined and repeatable without ambiguity. In leeAutomatedCraterDetection2019,
I attempted to provide enough information to reconstruct the catalogs and
provided detailed metric calculations and source code [Lee (20181)]. This
catalog also included the matching craters from Robbins2012a.
leeAutomatedCraterDetection2021 also provided software and catalogs necessary
to repeat the analysis leeAutomatedCraterDetection2021. Still, we did not
identify the matching ground-truth craters or the ‘true positives’ from CDA
catalogs. Similarly, benedixDerivingSurfaceAges2020 provided an accessible
catalog of craters but with no software or identifications of ground-truth
craters, and a more narrative description of the comparison methods.

Since these CDAs are now used to study craters smaller than existing catalog
limits, there are no ground-truth catalogs. It has become necessary to use
independent neural networks or other tools to provide supporting validation of
new catalogs. In leeAutomatedCraterDetection2021, we used multiple networks
trained on different input datasets to increase the overall performance. In
Lee2018a, I used multiple networks on the same high-resolution DTM dataset to
achieve the same result. In lagainMarsCraterDatabase2021, new web-based tools
were developed to validate an existing catalog and found several thousand new
craters and possible mis-identifications. These additional validation tools
from secondary networks or manual validation efforts should be considered
necessary in future studies to provide confidence in new crater catalogs.

Finally, since algorithm errors can arise, the source code and detailed
methodology for the entire CDA should be an essential part of future research
contributions. The source code and processed catalogs used in the analysis
here are provided at leeReplicationDataComparison2023.

## Open Research Section

The two published crater catalogs compared here can be found at
benedixDerivingSurfaceAges2020 and leeAutomatedCraterDetection2021 and links
therein. Software and catalog comparison data for this work is available at
https://doi.org/10.5683/SP3/AUJMHR.

###### Acknowledgements.

The author thanks benedixDerivingSurfaceAges2020 and Robbins2012a for
releasing their crater catalogs in easily accessible formats. The Natural
Sciences and Engineering Research Council of Canada funded part of this
project. The author thanks S. J. Robbins and an anonymous reviewer for
valuable comments that help clarify the article.

## References

  * Benedix . (2020) benedixDerivingSurfaceAges2020Benedix, GK., Lagain, A., Chai, K., Meka, S., Anderson, S., Norman, C.Tan, T.  202003\.  Deriving Surface Ages on Mars Using Automated Crater Counting Deriving Surface Ages on Mars Using Automated Crater Counting.  Earth and Space Science73.  10.1029/2019EA001005
  * Diakogiannis . (2019) Diakogiannis2019Diakogiannis, FI., Waldner, F., Caccetta, P. Wu, C.  2019\.  ResUNet-a: A Deep Learning Framework for Semantic Segmentation of Remotely Sensed Data ResUNet-a: A deep learning framework for semantic segmentation of remotely sensed data.  Arxiv.
  * Fairweather . (2022) fairweatherAutomaticMappingSmall2022Fairweather, JH., Lagain, A., Servis, K., Benedix, GK., Kumar, SS. Bland, PA.  202207\.  Automatic Mapping of Small Lunar Impact Craters Using LRO-NAC Images Automatic Mapping of Small Lunar Impact Craters Using LRO-NAC Images.  Earth and Space Science97.  10.1029/2021EA002177
  * Giannakis . (2023) giannakisDeepLearningUniversal2023Giannakis, I., Bhardwaj, A., Sam, L. Leontidis, G.  202304\.  Deep Learning Universal Crater Detection Using Segment Anything Model (SAM) Deep learning universal crater detection using Segment Anything Model (SAM) ( arXiv:2304.07764).  arXiv.  10.48550/arXiv.2304.07764
  * Goodfellow . (2016) Goodfellow2016Goodfellow, I., Bengio, Y. Courville, A.  2016\.  Deep Learning Deep Learning.  MIT Press.
  * Lagain, Benedix . (2021) lagainTharsisMantleSource2021Lagain, A., Benedix, GK., Servis, K., Baratoux, D., Doucet, LS., Rajšic, A.Miljković, K.  202111\.  The Tharsis Mantle Source of Depleted Shergottites Revealed by 90 Million Impact Craters The Tharsis mantle source of depleted shergottites revealed by 90 million impact craters.  Nature Communications1216352.  10.1038/s41467-021-26648-3
  * Lagain, Bouley . (2021) lagainMarsCraterDatabase2021Lagain, A., Bouley, S., Baratoux, D., Marmo, C., Costard, F., Delaa, O.Gamblin, O.  202108\.  Mars Crater Database: A Participative Project for the Classification of the Morphological Characteristics of Large Martian Craters Mars Crater Database: A participative project for the classification of the morphological characteristics of large Martian craters.  WU. Reimold C. Koeberl (), Large Meteorite Impacts and Planetary Evolution VI Large Meteorite Impacts and Planetary Evolution VI ( 550, 0).  Geological Society of America.  10.1130/2021.2550(29)
  * Latorre . (2023) latorreTransferLearningRealtime2023Latorre, F., Spiller, D., Sasidharan, ST., Basheer, S. Curti, F.  202304\.  Transfer Learning for Real-Time Crater Detection on Asteroids Using a Fully Convolutional Neural Network Transfer learning for real-time crater detection on asteroids using a Fully Convolutional Neural Network.  Icarus394115434.  10.1016/j.icarus.2023.115434
  * Lee (20181) leeCraterCatalogsSoftware2018Lee, C.  20181\.  Crater Catalogs and Software for ”Automated Crater Detection on Mars Using Deep Learning”. Crater Catalogs and software for ”Automated crater detection on Mars using deep learning”.
  * Lee (20182) Lee2018aLee, C.  20182\.  Martian Crater Identification Using Deep Learning Martian Crater Identification Using Deep Learning.  American Geophysical Union Fall Meeting American Geophysical Union Fall Meeting ( P41D-3768).
  * Lee (2019) leeAutomatedCraterDetection2019Lee, C.  2019\.  Automated Crater Detection on Mars Using Deep Learning Automated crater detection on Mars using Deep Learning.  Planetary and Space Science.  10.1016/j.pss.2019.03.008
  * Lee (20231) leeReplicationDataComparisonLee, C.  20231\.  Replication Data for: Comparison of Automated Crater Catalogs for Mars from Benedix et al. (2020) and Lee and Hogan (2021). Replication Data for: Comparison of automated crater catalogs for Mars from Benedix et al. (2020) and Lee and Hogan (2021).  borealisdata.caDataverse.  10.5683/SP3/AUJMHR
  * Lee (20232) leeReplicationDataComparison2023Lee, C.  2023207\.  Replication Data for: Comparison of Automated Crater Catalogs for Mars from Benedix et al. (2020) and Lee and Hogan (2021). Replication Data for: Comparison of automated crater catalogs for Mars from Benedix et al. (2020) and Lee and Hogan (2021).  Borealis.  10.5683/SP3/AUJMHR
  * Lee Hogan (2020) leeCraterCatalogsSoftware2020Lee, C. Hogan, J.  2020\.  Crater Catalogs and Software for ”Automated Crater Detection with Human Level Performance”. Crater Catalogs and software for ”Automated crater detection with human level performance”.  10.5683/SP2/CFUNII
  * Lee Hogan (2021) leeAutomatedCraterDetection2021Lee, C. Hogan, J.  202102\.  Automated Crater Detection with Human Level Performance Automated crater detection with human level performance.  Computers and Geosciences147104645.  10.1016/j.cageo.2020.104645
  * Malin . (2007) malinContextCameraInvestigation2007Malin, MC., Bell III, JF., Cantor, BA., Caplinger, MA., Calvin, WM., Clancy, RT.Wolff, MJ.  2007\.  Context Camera Investigation on Board the Mars Reconnaissance Orbiter Context Camera Investigation on board the Mars Reconnaissance Orbiter.  Journal of Geophysical Research: Planets112E5.  10.1029/2006JE002808
  * Redmon . (2016) redmonYouOnlyLook2016Redmon, J., Divvala, S., Girshick, R. Farhadi, A.  201606\.  You Only Look Once: Unified, Real-Time Object Detection You Only Look Once: Unified, Real-Time Object Detection.  2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) ( 779–788).  10.1109/CVPR.2016.91
  * Redmon Farhadi (2018) redmonYOLOv3IncrementalImprovement2018Redmon, J. Farhadi, A.  201804\.  YOLOv3: An Incremental Improvement YOLOv3: An Incremental Improvement ( arXiv:1804.02767).  arXiv.  10.48550/arXiv.1804.02767
  * Robbins (2017) Robbins2017Robbins, SJ.  2017\.  Progress on Old and New: D $>$ 1km Crater Catalogs for Mars and Moon. Progress on old and new: D $>$ 1km crater catalogs for Mars and Moon.  Planetary Crater Consortium 8 Planetary Crater Consortium 8 ( 1–2).  10.1111/maps.12895.
  * Robbins Hynek (2012) Robbins2012aRobbins, SJ. Hynek, BM.  2012\.  A New Global Database of Mars Impact Craters Larger than 1 Km: 2. Global Crater Properties and Regional Variations of the Simple-to-Complex Transition Diameter A new global database of Mars impact craters larger than 1 km: 2. Global crater properties and regional variations of the simple-to-complex transition diameter.  Journal of Geophysical Research E: Planets11761–21.  10.1029/2011JE003967
  * Salamunićcar . (2012) Salamuniccar2012Salamunićcar, G., Lončarić, S. Mazarico, E.  2012\.  LU60645GT and MA132843GT Catalogues of Lunar and Martian Impact Craters Developed Using a Crater Shape-based Interpolation Crater Detection Algorithm for Topography Data LU60645GT and MA132843GT catalogues of Lunar and Martian impact craters developed using a Crater Shape-based interpolation crater detection algorithm for topography data.  Planetary and Space Science601236–247.  10.1016/j.pss.2011.09.003
  * Silburt . (2019) Silburt2019Silburt, A., Ali-Dib, M., Zhu, C., Jackson, A., Valencia, D., Kissin, Y.Menou, K.  2019\.  Lunar Crater Identification via Deep Learning Lunar crater identification via deep learning.  Icarus31727–38.  10.1016/j.icarus.2018.06.022
  * Yang Cai (2022) yangHighResolutionFeaturePyramid2022Yang, S. Cai, Z.  2022\.  High-Resolution Feature Pyramid Network for Automatic Crater Detection on Mars High-Resolution Feature Pyramid Network for Automatic Crater Detection on Mars.  IEEE Transactions on Geoscience and Remote Sensing601–12.  10.1109/TGRS.2021.3104925


# Enhancing Next-Generation Urban Connectivity: Is the Integrated HAPS-
Terrestrial Network a Solution?

Afsoon Alidadi Shamsabadi, Animesh Yadav, Halim Yanikomeroglu Afsoon Alidadi
Shamsabadi and Halim Yanikomeroglu are with Non-Terrestrial Networks Lab.,
Department of Systems and Computer Engineering, Carleton University, Ottawa,
ON, K1S 5B6, Canada (e-mail: {afsoonalidadishamsa, halim}@sce.carleton.ca).
Animesh Yadav is with the School of Electrical Engineering and Computer
Science, Ohio University, Athens, OH, 45701, USA (e-mail: yadava@ohio.edu).

###### Abstract

Located in the stratospheric layer of Earth’s atmosphere, high altitude
platform station (HAPS) is a promising network infrastructure, which can bring
significant advantages to sixth-generation (6G) and beyond wireless
communications systems by forming vertical heterogeneous networks (vHetNets).
However, if not dealt with properly, integrated networks suffer from several
performance challenges compared to standalone networks. In harmonized spectrum
integrated networks, where different tiers share the same frequency spectrum,
interference is an important challenge to be addressed. This work focuses on
an integrated HAPS-terrestrial network, serving users in an overlapped urban
geographic area, and formulates a fairness optimization problem, aiming to
maximize the minimum spectral efficiency (SE) of the network. Due to the
highly nonconvex nature of the formulated problem, we develop a rapid
converging iterative algorithm that designs the multiple-input multiple-output
(MIMO) beamforming weights and the user association scheme such that the
propagated inter- and intra-tier interference is managed. Simulation results
demonstrate the proposed algorithm’s superiority over standalone terrestrial
networks and scenario where only the beamforming weights are optimized.

###### Index Terms:

HAPS, vHetNets, interference, user association, MIMO beamforming, spectral
efficiency

##  I Introduction

The rapid increase in mobile traffic and the emergence of novel use cases have
led to a pressing need for network architectures that can cope with the
forthcoming network trend in sixth-generation (6G) and beyond wireless
communications systems. In this context, non-terrestrial networks (NTNs) have
emerged as a promising architecture for incorporating powerful technologies
into communication networks[1]. Particularly, High altitude platform stations
(HAPS) are potential network platforms for supporting mobile users and
operators in facilitating various use cases such as computing, sensing,
positioning, and more. HAPS advantages such as lower latency than satellites,
and wider coverage area than terrestrial base stations, makes it a promising
complementary platform to existing wireless communications
infrastructure[2][3]. HAPS can play the role of a super macro base station
[4], and be integrated with conventional terrestrial base stations to form
vertical heterogeneous networks (vHetNets). Clearly, vHetNets improve the
coverage of telecommunications networks due to HAPS and can significantly
improve the network capacity, especially in hotspot urban areas where
conventional terrestrial networks cannot meet the increasing traffic demand.

However, harmonized spectrum integrated networks, where different tiers share
the same frequency band, are prone to performance degradation due to inter-
tier interference propagation[5]. To address this issue, efficient
interference management and coordination techniques must be implemented in
such multi-tier networks. Several techniques have been proposed to mitigate
the interference challenge[6, 7, 8, 9, 10, 11, 12]. In [6], the authors
proposed a HAPS-based null broadening approach for interference suppression in
an integrated network. The authors in [7] maximized the HAPS coverage by
decreasing the interference in the coexistence with terrestrial networks by
optimizing the beamwidth and beam direction of HAPS antenna. Beamforming was
employed in [8] to mitigate interference from HAPS to terrestrial cells. In
[9], the authors proposed an interference coordination method for an
integrated HAPS-terrestrial network considering the traffic load distribution.
In [10], the authors developed a user association scheme for an integrated
HAPS-terrestrial network based on the deep Q-learning (DQL) approach
considering delayed channel state information (CSI). The authors in [11]
developed an alternative user association and beamforming algorithm in an
integrated HAPS-terrestrial-satellite network, where HAPS plays the role of
the relay in satellite-UE connection. In [12], we developed an iterative
algorithm to design the subcarrier and transmit power allocation to user
equipment (UEs) to handle the interference in vHetNets.

Compared with [12], this work differs in four ways. First, we consider HAPS
and terrestrial macro base stations (MBSs) are equipped with multiple-input
multiple-output (MIMO) antennas, and hence, deploy beamforming scheme to form
narrow high-gain beams toward served UEs. Secondly, we assume all the UEs are
using the same time-frequency resource. Therefore, a proper design of the
beamforming weights is required to generate user-specific beams toward the UEs
and to control the allocated power to each UE, which is an important factor in
the interference management of vHetNets[12]. Thirdly, we also design a user
association scheme as it plays an important role in the interference
mitigation of the vHetNets, mainly due to the high transmit power and line-of-
sight (LoS) connection between HAPS and UEs. Lastly, we employed a realistic
channel model. Particularly, between the HAPS and UEs, a three-dimensional
(3D) Rician fading channel with a dominant LoS connection and non-LoS (NLoS)
parts is employed [14] and shadowing factor is included in the channel between
the MBSs and UEs.

The main contributions of this paper can be summarized as follows:

  * •

We formulated a max-min fairness (MMF) problem that maximizes the network’s
minimum spectral efficiency (SE) and designs the joint user association scheme
and beamforming weights.

  * •

Since the formulated problem is a mixed-integer non-linear program (MINLP), we
employed reformulation linearization techniques (RLT) to recast the problem to
an approximate convex format. Accordingly, we proposed a rapid converging
algorithm based on successive convex approximation (SCA) to solve the
approximate problem.

  * •

We provided simulation results to compare the performance of vHetNets, under
proposed algorithm, with baseline schemes and standalone terrestrial networks.
The results validated the performance improvement, achieved by integrating
HAPS with terrestrial networks, compared to the standalone terrestrial
networks.

The remainder of the paper is organized as follows. Section II presents the
system model and formulates the optimization problem, Section III details the
problem relaxation and explains the proposed algorithm, Section IV presents
and discusses the obtained numerical results, and finally, Section V concludes
the paper.

##  II System Model and Problem Formulation

Figure 1: Network architecture.

This paper considers a vHetNet deployed in an urban area, consisting of one
HAPS and $B$ MBSs, providing service to $U$ single antenna UEs within a
geographically overlapped urban region. The base stations111In this work, the
term base station refers to MBSs and HAPS, collectively. and UEs are indexed
by $b\in\\{1,\ldots,B+1\\}$ and $u\in\\{1,\ldots,U\\}$, respectively, with the
index $b=B+1$ reserved for HAPS. We assume that the vHetNet operates in the
downlink channel and all the UEs share the same time-frequency resource. The
considered system operates within the sub-6 GHz band, and each base station
$b$ is equipped with atotal number of $N_{b}=N^{\text{V}}_{b}\times
N^{\text{H}}_{b},$ antenna elements, where $N^{\text{V}}_{b}$ and
$N^{\text{H}}_{b}$ refer to the number of antenna elements in vertical and
horizontal axis, respectively. We refer to each base station’s coverage area
as a cell.

In this study, as illustrated in Fig. 1, each UE is allocated a dedicated spot
beam created by appropriate beamforming weights at the serving base station.
The beamforming matrix at base station $b$ is represented by
$\mathbf{W}^{b}=[\mathbf{w}^{b}_{1},\dots,\mathbf{w}^{b}_{U}]\in\mathbb{C}^{N_{b}\times
U}$, where
$\mathbf{w}^{b}_{u}=[w^{b}_{1,u},\dots,w^{b}_{N_{b},u}]^{T}\in\mathbb{C}^{N_{b}}$,
is defined as the beamforming weight vector of UE $u$ at base station $b$.
Additionally, we assume that the channel coefficients between each base
station $b$ and $U$ UEs form a channel matrix
$\mathbf{H}^{b}=[\mathbf{h}^{b}_{1},\dots,\mathbf{h}^{b}_{U}]\in\mathbb{C}^{N_{b}\times
U}$, where the channel vector between base station $b$ and UE $u$ is denoted
as
$\mathbf{h}^{b}_{u}=[h^{b}_{1,u},\dots,h^{b}_{N_{b},u}]^{T}\in\mathbb{C}^{N_{b}}$.
We assume that perfect CSI is available at each base station.

Considering the characteristics of urban terrestrial networks where the
connections between MBS and UE have dominant NLoS and negligible LoS
components, we consider small-scale fading, free-space path loss (FSPL), and
shadowing for the channel between each MBS and UE. In this regard, the channel
coefficient between antenna element $r\in\\{1,\ldots,N_{b}\\}$ of MBS $b$, and
UE $u$, denoted as $h^{b}_{r,u}$, can be formulated as follows[13]:

$h^{b}_{r,u}=\frac{\hat{h}^{b}_{r,u}\xi^{b}_{u}}{\sqrt{{PL}_{b,u}}},~{}\forall
r,~{}\forall u,~{}\forall~{}b\in\\{1,\dots,B\\},$ (1)

where $\hat{h}_{r,u}^{b}\sim\mathcal{NC}(0,1)$, normal random variable with
zero mean and unit variance, represents the small-scale Rayleigh fading
channel gain, and $\xi^{b}_{u}=10^{{{\xi_{u}^{{}^{\prime}b}}}/10}$ denotes the
log-normal shadowing gain, where ${\xi_{u}^{{}^{\prime}b}}$ is the Gaussian
random variable with zero mean and standard deviation $\sigma_{\xi}$ (in dB).
${PL}_{b,u}$ refers to the path loss between base station $b$ and UE $u$,
calculated as ${{PL}_{b,u}}=({4\pi f_{c}d_{b,u}}/{c})^{2},~{}\forall
u,~{}\forall b,$ where $f_{c}$ is the carrier frequency (in Hz), $d_{b,u}$ (in
m) is the distance between the base station $b$ and UE $u$, and $c$ represents
the speed of light in free space.

Considering the position of HAPS, located at an altitude 20 Km above the
ground, $\mathbf{h}^{B+1}_{u}$, representing the channel vector between HAPS
and UE $u$, can be modeled as a 3D Rician fading channel with a dominant LoS
and additional NLoS parts as [14]

$\mathbf{h}_{u}^{B+1}=\frac{1}{\sqrt{{PL}_{(B+1),u}}}(\sqrt{\frac{1}{1+K_{u}}}\hat{\mathbf{h}}_{u}+\sqrt{\frac{K_{u}}{1+K_{u}}}\overline{\mathbf{h}}_{u}),~{}\forall
u,$ (2)

where ${PL}_{(B+1),u}$ represents the FSPL between HAPS and UE $u$, $K_{u}$ is
the Rician factor for UE $u$, and
$\hat{\mathbf{h}}_{u}\in\mathbb{C}^{N_{B+1}}$ represents the NLoS component of
the channel vector with its elements from a normal random distribution with
zero mean and unit variance, $\mathcal{NC}(0,1)$. Accordingly,
$\overline{\mathbf{h}}_{u}$ is the LoS component of the channel vector as

$\overline{\mathbf{h}}_{u}=\mathbf{a}(\theta_{u},\phi_{u})\otimes\mathbf{b}(\theta_{u},\phi_{u}),~{}\forall
u,$ (3)

where the symbol $\otimes$ denotes the Kronecker product of two vectors, and
$\mathbf{a}(\theta_{u},\phi_{u})$ and $\mathbf{b}(\theta_{u},\phi_{u})$ are
steering vectors defined as

$\displaystyle\mathbf{a}(\theta_{u},\phi_{u})=[1,e^{j2\pi
d_{h}},\dots,e^{j2\pi(N^{\text{H}}_{B+1}-1)d_{h}}]^{T},~{}\forall u,$ (4a)
$\displaystyle\mathbf{b}(\theta_{u},\phi_{u})=[1,e^{j2\pi
d_{v}},\dots,e^{j2\pi(N^{\text{V}}_{B+1}-1)d_{v}}]^{T},~{}\forall u,$ (4b)

where $d_{h}=d_{\text{x}}\cos{\theta_{u}}\sin{\phi_{u}}/\lambda$ and
$d_{v}=d_{\text{y}}\cos{\theta_{u}}\cos{\phi_{u}}/\lambda$.
$\theta_{u}\in[0,\pi/2]$ and $\phi_{u}\in[-\pi,\pi)$ are the elevation and
azimuth angles of UE $u$. $d_{\text{x}}$ and $d_{\text{y}}$ are the antenna
elements spacing in $x$ and $y$ directions (as shown in Fig. 1), and
$\lambda=c/f_{c}$ represents the wavelength.

Now, we formulate the MMF problem. To this end, first, we define the binary
user association matrix, $\mathbf{A}\in~{}\mathbb{B}^{(B+1)\times U}$, whose
elements are defined as

${}{a_{b,u}}=\begin{cases}{1,}&{\text{if UE $u$ is associated to base station
$b$}},\\\ {0,}&{\text{otherwise.}}\end{cases}$ (5)

For the considered system model, the received signal of UE $u$ can be
expressed as

$y_{\text{u}}=\sum_{b=1}^{B+1}{{\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{u}}}s_{u}+\sum_{b=1}^{B+1}{\sum_{\begin{subarray}{c}k=1\\\
k\not=u\end{subarray}}^{U}{{\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{k}}s_{k}}}+n_{u},$
(6)

where $s_{u}$ refers to the transmitted data signal of UE $u$. In equation
(6), the first and second terms refer to the desired and interference signals,
received by UE $u$, respectively. Moreover, $n_{u}$ denotes the additive white
Gaussian noise (AWGN) at UE $u$, with zero mean and variance of
$\sigma^{2}_{n}$. $(\cdot)^{H}$ refers to the conjugate transpose operation.
Henceforth, the signal-to-interference-plus-noise ratio (SINR) for UE $u$ can
be defined as

$\gamma_{u}=\cfrac{\sum_{b=1}^{B+1}{\left|\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{u}\right|^{2}}}{\sum_{b=1}^{B+1}\sum_{\begin{subarray}{c}k=1\\\
k\not=u\end{subarray}}^{U}{\left|\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{k}\right|}^{2}+\sigma^{2}_{n}},~{}\forall
u,$ (7)

where $\mathbf{w}^{b}_{u}\in\mathbb{C}^{N_{b}}$ represents the beamforming
vector for UE $u$ at base station $b$, whereas
$\mathbf{h}^{b}_{u}\in\mathbb{C}^{N_{b}}$ represents the channel vector
between base station $b$ and UE $u$. It is noteworthy that the equations (6)
and (7) are applicable under the assumption that the beamforming weights for
non-associated UEs are zero at each base station. Now, we formulate the
fairness optimization problem for joint user association and beamforming
matrix design. The objective is to maximize the minimum SE of the network,
that is equivalent to maximizing the minimum SINR of the network. The initial
formulation of the optimization problem can be represented mathematically as

$\displaystyle\underset{\mathbf{A},~{}\mathbf{W},{{\gamma}_{m}}}{\text{maximize}}\,\,$
$\displaystyle~{}~{}\gamma_{m}\,$ s.t.
$\displaystyle\gamma_{u}\geq{\gamma_{m}},~{}\forall u,$
$\displaystyle\|\mathbf{w}^{b}_{u}\|^{2}_{2}\leq
P^{\text{max}}_{b}a_{b,u},\,~{}\forall u,~{}\forall b,\qquad$
$\displaystyle\sum_{u=1}^{U}a_{(B+1),u}F\log_{2}(1+\gamma_{u})\leq{R_{\text{BH}}},\qquad$
$\displaystyle{\|\mathbf{W}^{b}\|^{2}_{F}}\leq P^{\text{max}}_{b},\,~{}\forall
b,\qquad$ $\displaystyle\sum_{b=1}^{B+1}{a_{b,u}}\geq 1,\,~{}\forall u,$
$\displaystyle\mathbf{A}\in\mathbb{B}^{(B+1)\times U},$
$\displaystyle\mathbf{W}^{b}\in\mathbb{C}^{N_{b}\times U},~{}\forall b,\qquad$

where the objective function, ${\gamma}_{m}$, represents minimum SINR of
network, and $\mathbf{W}$ denotes the collection of all beamforming matrices
$\mathbf{W}^{b}\in\mathbb{C}^{N_{b}\times U},~{}\forall b$. $F$ represents the
allocated bandwidth to each UE, $R_{\text{BH}}$ denotes the achievable rate of
the backhaul link between HAPS and the ground station, and
$P_{b}^{\text{max}},~{}\forall b,$ is the maximum available transmit power at
base station $b$. $\|\cdot\|_{F}$ refers to the matrix Frobenius norm
operation.

Problem (II) aims to jointly design the user association and beamforming
matrices, under a set of constraints (II)-(II). Constraint (II) ensures that
the MMF is maintained within the network, while constraint (II) enforces that
at each base station, the beamforming weights for non-associated UEs are zero.
Constraint (II) ensures that the sum of the data rate for the UEs, associated
with HAPS, does not exceed the achievable rate of the HAPS backhaul link.
Constraint (II) limits the total transmit power of each base station to its
maximum available transmit power. Finally, constraint (II) guarantees that
each UE gets associated with at least one base station222In this study, there
is no constraint on the number of associated UEs with each base station..
Problem (II) is an MINLP problem due to non-convex constraints (II) and (II),
and binary variables in (II). The presence of binary variable matrix
$\mathbf{A}$ makes the problem challenging to be solved for the global
solution. In the following section, we propose a rapid converging SCA-
based[15] iterative algorithm to obtain the suboptimal solution for problem
(II). The approach utilizes RLT to transform MINLP into an equivalent
tractable problem, which is then solved iteratively until convergence to the
suboptimal solution.

##  III Problem Relaxation and Proposed Algorithm

This section presents a methodology to transform problem (II) into a
mathematically tractable form by applying various approximations and
techniques to deal with the non-convex parts. These techniques are summarized
in the following of this section. First, we deal with the non-convex
constraint (II). Due to the fractional formula on the left-hand side, the
constraint is non-convex. To tackle this, we define a few set of slack
variables $\beta_{u},~{}\forall u$, and $\alpha_{u},~{}\forall u$, and replace
the left-hand side by its lower bound approximation $\alpha_{u}$. Two new
constraints, as (9a), should be defined to ensure that $\alpha_{u}$ is the
lower bound approximation for $\gamma_{u}$:

$\displaystyle\alpha_{u}\beta_{u}\leq\sum_{b=1}^{B+1}{\left|\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{u}\right|^{2}},~{}\forall
u,$ (9a)
$\displaystyle\beta_{u}\geq\sum_{b=1}^{B+1}\sum_{\begin{subarray}{c}k=1\\\
k\not=u\end{subarray}}^{U}{\left|\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{k}\right|}^{2}+\sigma^{2}_{n},\,~{}\forall
u.$

Note that, equation (III) represents a convex inequality constraint but
equation (9a) is non-convex due to the norm at the right hand side. To handle
this, we define slack variable matrices $\mathbf{P}\in\mathbb{R}^{(B+1)\times
U},$ and $\mathbf{Q}\in\mathbb{R}^{(B+1)\times U},$ with their elements,
respectively, defined as

$\displaystyle
p_{b,u}\leq\Re{\left(\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{u}\right)},~{}\forall
u,~{}\forall b,\qquad$ (10a) $\displaystyle
q_{b,u}\leq\Im{\left(\left(\mathbf{h}^{b}_{u}\right)^{H}\mathbf{w}^{b}_{u}\right)},~{}\forall
u,~{}\forall b,\qquad$

where, $\Re(\cdot)$ and $\Im(\cdot)$ represent the real and imaginary parts of
a complex number, respectively. Consequently, equation (9a) can be
reformulated as

$\alpha_{u}\beta_{u}\leq{\sum_{b=1}^{B+1}{(p_{b,u})^{2}+(q_{b,u})^{2}}},~{}\forall
u.$ (11)

In order to be able to convert equation (11) to a convex form, we rewrite the
equation as

$\alpha_{u}\leq\sum_{b=1}^{B+1}{\frac{(p_{b,u})^{2}+(q_{b,u})^{2}}{\beta_{u}}},~{}\forall
u,$ (12)

which can be converted to equivalent convex form by replacing the right hand
side with its approximated first-order Taylor series as

${}\alpha_{u}\leq\sum_{b=1}^{B+1}\frac{2{p^{(n)}_{b,u}}}{{\beta^{(n)}_{u}}}\left(p_{b,u}-{p^{(n)}_{b,u}}\right)+\frac{2{q^{(n)}_{b,u}}}{{\beta^{(n)}_{u}}}\left(q_{b,u}-{q^{(n)}_{b,u}}\right)\\\
+\frac{\left({p^{(n)}_{b,u}}\right)^{2}+\left({q^{(n)}_{b,u}}\right)^{2}}{\beta^{(n)}_{u}}\left(1-\frac{\beta_{u}-{\beta^{(n)}_{u}}}{{\beta^{(n)}_{u}}}\right),~{}\forall
u,$ (13)

where superscript $(n)$ refers to the value of corresponding variable in the
$n$th iteration of the SCA process. Next, we tackle constraint (II) by
acknowledging that at the convergence, $\gamma_{u},\forall u,$ will converge
to $\gamma_{m}$[16], and by employing Jensen’s
inequality333$\mathbb{E}\\{\log(f(x))\\}\leq\log(\mathbb{E}\\{f(x)\\})$, where
$\mathbb{E}$ denotes the expected value.. Consequently, we can replace
constraint (II) as

$\sum_{u=1}^{U}a_{(B+1),u}\gamma_{m}\leq{2^{R_{\text{BH}}/{F}}-1}.$ (14)

Furthermore, we define slack variables $\eta_{u}=a_{(B+1),u}\gamma_{m},\forall
u,$ and define four new constraints as below to ensure that the actual
equality holds:

$\displaystyle 0\leq\eta_{u}\leq\gamma_{\text{max}}a_{(B+1),u},~{}\forall u,$
(15a) $\displaystyle
0\leq\gamma_{m}-\eta_{u}\leq\gamma_{\text{max}}(1-a_{(B+1),u}),~{}\forall u,$

where $\gamma_{\text{max}}$ is enough large number denoting the upper bound
for $\gamma_{m}$. Consequently, constraint (14) will be transformed as

$\sum_{u=1}^{U}\eta_{u}\leq{2^{R_{\text{BH}}/{F}}-1}.$ (16)

After applying the aforementioned relaxations, the SCA optimization problem at
$n$th iteration will be as

$\displaystyle\underset{\mathbf{W},\mathbf{A},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{P},\mathbf{Q},\boldsymbol{\eta},\gamma_{m}}{\text{maximize}}\,\,$
$\displaystyle\gamma_{m}\,$ (17a) s.t.
$\displaystyle\alpha_{u}\geq\gamma_{m},~{}\forall u,\,$ (17d)
$\displaystyle\eqref{eq:P1_const2},\eqref{eq:P1_const3},\eqref{eq:P1_const4},\eqref{eq:P1_const5},\eqref{eq:P1_const6},$
$\displaystyle\eqref{eq:betaigeqIi},\eqref{Real and Imag},\eqref{first
Taylor},\eqref{eq:BH_1},\eqref{eq:BH_2},\qquad$

where $\boldsymbol{\alpha}=[\alpha_{1},\dots,\alpha_{U}]^{T}$,
$\boldsymbol{\beta}=[\beta_{1},\dots,\beta_{U}]^{T}$, and
$\boldsymbol{\eta}=[\eta_{1},\dots,\eta_{U}]^{T}$. Problem (17a) is a second
order convex programming (SOCP) type problem which can be solved using solvers
such as Mosek. The pseudocode of the proposed iterative joint user association
and beamforming weights design (JUBD) algorithm is outlined in Algorithm 1.
According to the algorithm, (17a) is solved iteratively until the difference
ratio of the objective function value is less than the tolerance $10^{-4}$, or
maximum number of SCA iterations, $N_{\text{iter}}$, is reached, whichever
first. Since the objective function of problem (17a) is bounded above due to
the maximum transmit power constraint (II), and each iteration of Algorithm 1
generates a sequence of improved objective function values, Algorithm 1
converges to a stationary point of problem (II)[15]. In Algorithm 1, the
initial values are chosen from the feasible set. The worst-case computational
cost of the algorithm is mainly determined by the branch-and-bound method to
solve MIP at each SCA iteration, and is given as
$\mathcal{O}(2^{U\left(B+1\right)})$.

Algorithm 1 Proposed iterative JUBD algorithm for vHetNets.

1: Input: $U,~{}B,~{}N_{b},~{}\forall b,~{}\mathbf{H}^{b},~{}\forall
b,~{}P_{b}^{\text{max}},~{}\forall
b,~{}\sigma^{2}_{n},~{}F,~{}R_{\text{BH}},~{}N_{\text{iter}}=10$.

2: Output: $\mathbf{W}^{*},~{}\mathbf{A}^{*}$.

3: Initialize
$\boldsymbol{\beta}^{(0)},~{}\mathbf{P}^{(0)},~{}\mathbf{Q}^{(0)},$ and set
$n:=0.$

4: Repeat until $|\gamma_{m}^{(n)}-\gamma_{m}^{(n-1)}|/\gamma_{m}^{(n-1)}\leq
10^{-4}$ or until $n\leq N_{\text{iter}}:$ \- solve (17a) to find
$~{}\mathbf{A}^{(n)*},~{}\mathbf{W}^{(n)*},~{}\boldsymbol{\alpha}^{(n)*},~{}\boldsymbol{\beta}^{(n)*},~{}\mathbf{Q}^{(n)*},$$\mathbf{P}^{(n)*},~{}\boldsymbol{\gamma}_{m}^{(n)*}$.
\- update
$\boldsymbol{\beta}^{(n+1)}:=\boldsymbol{\beta}^{(n)*},~{}\mathbf{P}^{(n+1)}:=\mathbf{P}^{(n)*},$$\mathbf{Q}^{(n+1)}:=\mathbf{Q}^{(n)*},~{}n:=n+1$.

##  IV Numerical Results

TABLE I: Simulation Parameters. Parameter | Value
---|---
Center frequency ($f_{c}$) |  $2.545$ GHz
Shadowing standard deviation ($\sigma_{\xi}$), Rician factor ($K_{u}$) |  $8$, $10$
Number of antenna elements ($N_{\text{b}},~{}b\in\\{1,\dots,B\\}$, $N_{B+1}$) |  $4\times 4$, $8\times 8$
HAPS antenna elements spacing ($d_{x},~{}d_{y}$) | $\lambda/2$
Maximum transmit power ($P_{b}^{\text{max}},~{}b\in\\{1,\dots,B\\}$, $P_{B+1}^{\text{max}}$) |  $43$ dBm, $52$ dBm
UE allocated bandwidth ($F$), HAPS backhaul rate ($R_{\text{BH}}$) |  $1$ MHz, $20$ Gbps
AWGN variance ($\sigma^{2}_{n}$) |  $-100$ dBm
Figure 2: Convergence behaviour of the proposed Algorithm 1.

(a) CDF of minimum spectral efficiency for 5 Scenarios.

(b) Impact of HAPS antenna configuration in Scenario 1.

Figure 3: Statistical behavior of minimum spectral efficiency. Scenario 1:
vHetNet (4 MBSs + HAPS) using JUBD Algorithm. Scenario 2: vHetNet (4 MBSs +
HAPS) beamforming optimization with max-SINR based user association. Scenario
3: Standalone terrestrial network (4 MBSs) using JUBD Algorithm. Scenario 4:
Standalone terrestrial network (5 MBSs) using JUBD Algorithm. Scenario 5:
vHetNet (4 MBSs + HAPS) using JUBD Algorithm with multiple base station
association. Figure 4: Probability distribution of UEs associated with HAPS.

In this section, we evaluate the performance of the proposed JUBD Algorithm 1
through numerical simulations. In simulations, we consider a $4$ Km by $4$ Km
square hotspot urban area with $16$ uniformly distributed UEs in it. To
facilitate comparisons, we consider five distinct scenarios. Scenario 1
involves a vHetNet comprising four MBSs and one HAPS, and employs JUBD
algorithm 1. Scenario 2 also involves a vHetNet comprising four MBSs and one
HAPS but focuses only on optimizing the beamforming weights and employs the
maximum received SINR-based scheme for user association. Scenarios 3 and 4
consider a standalone terrestrial network comprising four and five MBSs,
respectively, and employ JUBD Algorithm 1. Scenarios 1 to 4 assume that each
UE can get associated to only one base station; however, we consider Scenario
5 where Algorithm 1 is implemented in a vHetNet (one HAPS and four MBSs) with
multiple base station association, where UEs can get associated with more than
one base station. Another potential scheme to compare with our proposed JUBD
Algorithm 1 is the DQL method[10] which may not be the best for the problem
setting in this paper due to the larger scale of the network and a larger
number of variables. The results are obtained for each scenario based on
$1000$ independent and identically distributed (i.i.d.) network topologies.
All simulations are implemented in MATLAB using CVX with Mosek 9.1.9 as the
internal solver. The rest of the simulation parameters are listed in Table I.

Fig. 2 illustrates the evolution of objective function values across SCA
iterations for three randomly sampled network topologies representing UE
distributions. It can be observed that the objective function value increases
with each SCA iteration and saturates after five or six iterations. Fig. 3
illustrates the cumulative distribution function (CDF) of the minimum SE for
different scenarios. Particularly, Fig. 3(a) compares Scenarios 1, 2, 3, 4,
and 5, and from it, we can deduce four observations. First, the proposed
vHetNet architecture (i.e., Scenarios 1, 2, and 5), provides improved SE
performance as compared to those of standalone terrestrial networks (i.e.,
Scenarios 3 and 4). This improvement is due to the fact that vHetNets allow
some UEs to connect to HAPS with higher order MIMO antenna, LoS links and
better channels, resulting in improved SE performance. Fig. 4, which plots the
probability distribution of UEs associated with HAPS, corroborates this fact
that in every sample network topology, a few number of UEs (less than seven
UEs out of 16) are associated with HAPS.

Secondly, the proposed JUBD Algorithm 1 for vHetNet, used in Scenarios 1 and
5, yields better SE performance compared to the scenario where only the
beamforming weights are optimized (i.e., Scenario 2). Thirdly, the CDF curve
of vHetNet (i.e., Scenarios 1, 2, and 5) exhibited a notably steeper slope
when compared to those of standalone terrestrial networks (i.e., Scenarios 3
and 4). Steepness indicates a reduced standard deviation concerning the
minimum SE. Accordingly, in vHetNets, UEs experience almost the same quality-
of-service (QoS) regardless of their positioning within the coverage area.
This is because if a UE experiences suboptimal channel conditions within the
terrestrial network, it gets associated with HAPS with a better channel
quality due to the dominant LoS connection. The third observation reaffirms
the provision of fairer QoS to UEs in different network topologies. Lastly,
comparing Scenarios 1 and 5, it can be observed that Scenario 5, where each UE
can get associated with multiple base stations in a vHetNet, leads to a
slightly higher SE compared to Scenario 1, where UEs can get associated with
only one base station.

Fig. 3(b) plots the CDF of minimum SE for Scenario 1 with the varying number
of antenna elements on HAPS, $N_{B+1}=\\{8\times 8,~{}16\times 16,~{}32\times
32\\}$. It can be observed that the SE performance of the network increases
with the number of antenna elements on HAPS. This is because a higher order
MIMO antenna easily forms a narrower spot beam, resulting in increased gain
toward the target UE and decreased interference toward other UEs.

##  V Conclusion

In this work, we proposed the design of an optimum beamforming weight and user
association scheme for an integrated HAPS-terrestrial network (vHetNet) to
mitigate interference. We formulated an MMF optimization problem, which turned
out to be a MINLP. In general, MINLP problems are challenging to solve; thus,
we reformulated the original problem and developed a low-complexity and fast-
converging SCA-based iterative algorithm to solve it. Simulation results
showed that vHetNets with optimized beamforming and user association scheme
offer higher SE than standalone terrestrial network scenarios. This
observation demonstrates the potential of HAPS-enabled vHetNets for the 6G and
beyond wireless networks. The interference management in vHetNets, based on
machine learning approaches, opens up future research directions.

## References

  * [1] Y. Hokazono, H. Kohara, Y. Kishiyama and T. Asai, ”Extreme coverage extension in 6G: Cooperative non-terrestrial network architecture integrating terrestrial networks,” in Proc. IEEE Wireless Commun. Netw. Conf. (WCNC), Austin, TX, USA, 10–13 Apr. 2022, pp. 138-143.
  * [2] G. Karabulut Kurt, M. G. Khoshkholgh, S. Alfattani, A. Ibrahim, T. S. J. Darwish, M. S. Alam, H. Yanikomeroglu, and A. Yongacoglu, “A vision and framework for the high altitude platform station (HAPS) networks of the future,” IEEE Commun. Surveys Tuts., vol. 23, no. 2, pp. 729–779, Secondquarter 2021.
  * [3] T. Tozer and D. Grace, “High-altitude platforms for wireless communications,” Electron. Commun. Eng. J., vol. 13, no. 3, pp. 127–137, Jun. 2001.
  * [4] M. S. Alam, G. K. Kurt, H. Yanikomeroglu, P. Zhu, and N. D. Dào, ”High altitude platform station based super macro base station constellations,” IEEE Commun. Mag., vol. 59, no. 1, pp. 103-109, Jan. 2021.
  * [5] Z. Wei, L. Wang, Z. Gao, H. Wu, N. Zhang, K. Han, Z. Feng, ”Spectrum sharing between high altitude platform network and terrestrial network: Modeling and performance analysis,” IEEE Trans. Commun., vol. 71, no. 6, pp. 3736-3751, Jun. 2023.
  * [6] Y. Kawamoto, A. Matsushita, S. Verma, N. Kato, K. Kaneko, A. Sata, M. Hangai, ”HAPS-based interference suppression through null broadening with directivity control in space-air-ground integrated networks,” IEEE Trans. Veh. Technol., vol. 72, no. 12, pp. 16098-16107, Dec. 2023.
  * [7] Y. Shibata, W. Takabatake, K. Hoshino, A. Nagate and T. Ohtsuki, ”HAPS cell design method for coexistence on terrestrial mobile networks,” in Proc. IEEE Veh. Technol. Conf., Florence, Italy, 20-23 Jun. 2023, pp. 1-6.
  * [8] K. Tashiro, K. Hoshino, and A. Nagate, ”Nullforming-based precoder for spectrum sharing between HAPS and terrestrial mobile networks,” IEEE Access, vol. 10, pp. 55675-55693, 2022.
  * [9] W. Liu, X. Hou, L. Chen, Y. Hokazono and J. Zhao, ”Interference coordination method for integrated HAPS-terrestrial networks,” in Proc. IEEE Veh. Technol. Conf., Helsinki, Finland, 19-22 Jun. 2022, pp. 1-6.
  * [10] H. Khoshkbari, S. Sharifi and G. Kaddoum, ”User association in a VHetNet with delayed CSI: A deep reinforcement learning approach,” IEEE Commun. Lett., vol. 27, no. 8, pp. 2257-2261, Aug. 2023.
  * [11] S. Liu, H. Dahrouj and M.-S. Alouini, ”Joint user association and beamforming in integrated satellite-HAPS-ground networks,” IEEE Trans. Veh. Technol., 2023 (Early Access).
  * [12] A. A. Shamsabadi, A. Yadav, O. Abbasi and H. Yanikomeroglu, ”Handling interference in integrated HAPS-terrestrial networks through radio resource management,” IEEE Wireless Commun. Lett., vol. 11, no. 12, pp. 2585-2589, Dec. 2022.
  * [13] U. Ghafoor, H. Z. Khan, M. Ali, A. M. Siddiqui, M. Naeem and I. Rashid, ”Energy efficient resource allocation for H-NOMA assisted B5G HetNets,” IEEE Access, vol. 10, pp. 91699-91711, 2022.
  * [14] Z. Lian, L. Jiang, C. He and D. He, ”User grouping and beamforming for HAP massive MIMO systems based on statistical-eigenmode,” IEEE Wireless Commun. Lett., vol. 8, no. 3, pp. 961-964, Jun. 2019.
  * [15] A. Beck, A. Ben-Tal, and L. Tetruashvili, ”A sequential parametric convex approximation method with applications to nonconvex truss topology design problems,” J. Global Optim., vol. 47, no. 1, pp. 29–51, May 2010.
  * [16] Y. Mao, B. Clerckx, J. Zhang, V. O. K. Li and M. A. Arafah, ”Max-min fairness of K-user cooperative rate-splitting in MISO broadcast channel With user relaying,” IEEE Trans. Wireless Commun., vol. 19, no. 10, pp. 6362-6376, Oct. 2020.