## 2 System Model and Problem Formulation  
### System Model  
As shown in Figure 1, an uplink communication scenario for the terrestrial-satellite network is investigated, in which \(M\) BSs and one satellite (The proposed algorithm can be easily extended to multiple satellite settings. The additional satellites can be viewed as additional APs.) are deployed to serve \(K\) terrestrial UEs. The UEs' data is uploaded to the BSs or satellite and then transmitted to the gateway of the core network via the backhaul links. Define the sets of UEs and BSs as \(\mathcal{K}=\{1,2,\ldots,K\}\) and \(\mathcal{M}=\{1,2,\ldots,M\}\), respectively. The satellite is denoted as \(M+1\). And all the APs (i.e., BSs and satellite) constitute a set denoted by \(\mathcal{R}=\mathcal{M}\cup M+1\). Assume that the \(m\)-th AP is equipped with \(N_{m}\), \(m=1,\ldots,M+1\) antennas, while the terrestrial UEs are all deployed with single antennas. Denote \(\mathcal{U}_{m}\) as the UEs associated to AP \(m\), satisfying \(\mathcal{K}=\{\mathcal{U}_{1}\cup\mathcal{U}_{2}\cdots\cup\mathcal{U}_{M+1}\}\), then the received signal of AP \(m\) from UE \(k\) is denoted by  
\[\mathbf{y}_{m,k}=\mathbf{h}_{m,k}\sqrt{p_{m,k}}s_{m,k}+\sum_{j\in\mathcal{U}_{ m}\backslash k}\mathbf{h}_{m,j}\sqrt{p_{m,j}}s_{m,j}+\sum_{i=1,i,j\
eq m}^{M+1} \sum_{j\in\mathcal{U}_{i}}\mathbf{h}_{m,j}\sqrt{p_{i,j}}s_{i,j}+\mathbf{n}_{m}, \tag{1}\]  
where \(p_{m,k}\) and \(s_{m,k}\) denote the transmit power and symbol of UE \(k\) served by AP \(m\). \(\mathbf{h}_{m,k}\in\mathbb{C}^{N_{m}\times 1}\) denotes the channel vector of UE \(k\) towards AP \(m\). \(\mathbf{n}_{m}\sim\mathcal{CN}(0,\sigma_{m}^{2}\mathbf{I})\) is the additive white Gaussian noise (AWGN) vector. Similar to most of the related works [4; 18], we assume that all the wireless channels undergo slow fading, and the perfect channel state information of all channels are available at the gateway through feedback/training sent from the user terminals via the backhaul channel. This mechanism has already been adopted in DVB-S2 [19]. For the uplink NOMA communication scenario, the UEs with better channel conditions have the higher priority to be decoded through SIC [2]. Assume that channel coefficients of UEs towards AP \(m\) are sorted as  
\[Q_{m}\triangleq\left|\mathbf{h}_{m,\pi_{m}(1)}\right|^{2}\geq\left|\mathbf{h }_{m,\pi_{m}(2)}\right|^{2}\geq\cdots\geq\left|\mathbf{h}_{m,\pi_{m}(K)} \right|^{2},\forall m\in\mathcal{R}. \tag{2}\]Let \(\mathcal{S}_{m,k}=\{\pi_{m}(1),\pi_{m}(2),\cdots,\pi_{m}(k)\}\); thus, the received SINR of UE \(k\) at AP \(m\) is calculated by  
\[SINR_{m,k}=\frac{\left|\mathbf{v}_{m,k}\mathbf{h}_{m,k}\right|^{2}p_{m,k}}{ \sum\limits_{j\in\mathcal{U}_{m}\backslash\mathcal{S}_{m,k}}\left|\mathbf{v}_ {m,k}\mathbf{h}_{m,j}\right|^{2}p_{m,j}+\sum\limits_{i=1,j+m}^{M+1}\sum \limits_{j\in\mathcal{U}_{i}}\left|\mathbf{v}_{m,k}\mathbf{h}_{m,j}\right|^{2}p _{i,j}+\left\|\mathbf{v}_{m,k}\right\|^{2}\sigma_{m}^{2}}, \tag{3}\]  
where \(\mathbf{v}_{m,k}\in\mathbb{C}^{1\times N_{m}}\) is the receive beamformer for AP \(m\) to detect the signal from UE \(k\). Thus, the achievable rate of UE \(k\) at AP \(m\) is given by  
\[R_{m,k}=\log(1+SINR_{m,k}) \tag{4}\]  
Define \(\alpha_{m,k}\) as the UE association indicator, where \(\alpha_{m,k}=1\) denotes that UE \(k\) is served by AP \(m\), and \(\alpha_{m,k}=0\), otherwise. Therefore, the achievable transmission rate of UE \(k\) is written as  
\[R_{k}=\sum\limits_{m\in\mathcal{R}}\alpha_{m,k}R_{m,k} \tag{5}\]  
And the achievable sum rate for AP \(m\) is given by  
\[R_{m}=\sum\limits_{k\in\mathcal{K}}\alpha_{m,k}R_{m,k} \tag{6}\]