Suppose the system $Ax=b$ has a unique solution $c \in \mathbb{Q}_{\geq 0}^n$, where $A \in \mathbb{Z}_{\geq 0}^{n \times n}$ and $b \in \mathbb{Z}_{\geq 0}^n$.  Note that the equation $x_1=c_1$ can be written as a linear combination of equations which appear in $Ax=b$.  Thus, there exist $q_1, \dots, q_n \in \mathbb{Q}$ such that $\sum_{i=1}^n q_iA_i=e_1$ and $\sum_{i=1}^n q_i b_i=c_1$, where $A_i$ is the $i$th row of $A$ and $e_1$ is the first standard basis vector.  Let $I$ be the set of indices $i$ such that $q_i \geq 0$.  Then for all $j \notin I$, the coefficients of $\sum_{i \in I} q_iA_i$ dominate the coefficients of $q_jA_j$, and $\sum_{i \in I} q_ib_i \geq q_jb_j$.  Thus, we may subtract these equations from $\sum_{i \in I} q_iA_i$ to derive $x_1=c_1$.  Then just repeat for each variable. 

Note that this method only needs to retain at most $n+1$ equations throughout. I am not sure if keeping only $n$ equations is possible.