在上一篇文章中，我们介绍了椭球模型下的一系列基础的坐标操作。本节，介绍观测与遮挡问题。
观测主要用于从观察点A观测大地标准点B，用来解决观测的仰角、方位角与大地坐标系之间的关系。
在没有GPS卫星的时代，为了测量一个位置的坐标，往往会设置多个采样点，不断测角、测距、测气压，“跑断腿”。在现代，这种基于方位俯仰的测量技术已经用的不多了，但方位俯仰的计算还是很有用的。
遮挡问题和观测是同一个问题，解决的是从A点能不能看到B的问题。有时候计算A,B的共视问题时用的很多，尤其是其中一方为高山或者飞机时。
1 根据AB两点位置求取方位角、俯仰角
在 A 点观测B点的方位俯仰，是在A点的ENU坐标下完成的。
接前文，A的ECEF坐标为
A
⃗
\vec A
A, B的ECEF坐标为
B
⃗
\vec B
B, 在A的ENU坐标系下的坐标为
B
E
N
U
⃗
\vec{B_{ENU}}
BENU，他们的关系：
要计算站在t位置下，查看任意一点 M 的 ENU 坐标，只要执行运算
B E N U ⃗ T = [ E ⃗ N ⃗ U ⃗ ] × ( B E C E F ⃗ − t ⃗ ) T = R × ( B E C E F ⃗ − t ⃗ ) T \vec {B_{ENU}}^T=\begin{bmatrix} {\vec {E} \\ \vec N \\ \vec U}\end{bmatrix} \times \left ( \vec {B_{ECEF}} - \vec t \right )^T =R \times \left ( \vec {B_{ECEF}} - \vec t \right )^T BENUT=[ENU]×(BECEF−t)T=R×(BECEF−t)T
为了行文方便，重复一下R矩阵的特点。上式中坐标都是行向量,转置后是列向量。R的左逆、右逆都是自身的转置,故而两类坐标的转换很方便。如果是求取单位向量，则把 t ⃗ \vec t t设置为0即可。
R = [ − sin θ cos θ 0 − sin φ c o s θ − sin φ sin θ cos φ cos φ c o s θ cos φ sin θ sin φ ] R = \begin{bmatrix} {-\sin \theta} & {\cos \theta} & 0\\ {-\sin \varphi cos\theta} & {-\sin \varphi \sin \theta} & {\cos \varphi} \\ \cos \varphi cos\theta & \cos \varphi \sin \theta & \sin \varphi \end{bmatrix} R= −sinθ−sinφcosθcosφcosθcosθ−sinφsinθcosφsinθ0cosφsinφ
R × R T = R T × R = I R\times R^T=R^T \times R=I R×RT=RT×R=I
(1) 俯仰El
E l = sin − 1 U ∣ B E N U ⃗ ∣ El = \sin^{-1} \frac{U}{\left| \vec {B_{ENU}}\right|} El=sin−1 BENU U
(2) 方位Az
A z = tan − 1 E N Az = \tan^{-1} \frac{E}{N} Az=tan−1NE
(3) 基于大气折射的俯仰修正
对于大气层外的B位置，还需要考虑电磁波的折射问题。我们参考 http://www.zeptomoby.com/satellites/ 中的代码，其修正原理如下：
#ifdef WANT_ATMOSPHERIC_CORRECTION
double saveEl = el;
// Elevation correction for atmospheric refraction.
// Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104
// Note: Correction is meaningless when apparent elevation is below horizon
el += deg2rad((1.02 /
tan(deg2rad(rad2deg(el) + 10.3 /
(rad2deg(el) + 5.11)))) / 60.0);
if (el < 0.0)
{
// Reset to true elevation
el = saveEl;
}
if (el > (PI / 2))
{
el = (PI / 2);
}
#endif
主要接口：
/*!
* \brief observe_az_el_range 计算方位角、俯仰角
* \param llaA_in A点(观察者)
* \param llaB_in B点(被观察者)
* \param paz 方位角
* \param pel 俯仰角
* \param prange 长度
* \param ATMOSPHERIC_CORRECTION 大气折射矫正开关
* \param rad 经纬度采用弧度
* \param latfirst 纬度在第一维度[0]
*/
inline void observe_az_el_range(
const double llaA_in[/*3*/],
const double llaB_in[/*3*/],
double * paz,
double * pel,
double * prange,
bool ATMOSPHERIC_CORRECTION = false,
const bool rad = false,
const bool latfirst = true
)
2 可视判断问题
上一节的问题，可以引申出一个命题。如果知道A\B的位置，如何判断A与B之间是不是发生遮挡？同学们有时候会想当然认为只要仰角小于0，就会遮挡。这其实只存在于A位于海面的情况。真实情况要复杂的多，和二者所处的地形有很大关系。
(1) 基于DEM的精确测算方法
一般认为：在海面上，只要A、B两点的连线与当前椭球面没有交点，则可视；在陆地上，只要A、B两点的连线与当前地形没有交点，则可视。在可视关系发生时，仰角完全可以小于0.
此问题可等效为线段A-B与地表数字高程DEM平面有无交点的问题。我们构造参数方程描述线段AB上的任意一点 E，游标
ϱ
\varrho
ϱ是介于0~1的实数。
E ⃗ = A ⃗ + ( B ⃗ − A ⃗ ) ϱ \vec E=\vec A + (\vec B-\vec A)\varrho E=A+(B−A)ϱ
设函数
h
=
F
d
e
m
(
E
⃗
)
h=F_{dem}(\vec E)
h=Fdem(E)能够返回ECEF坐标
E
⃗
\vec E
E对应的椭球切投影处的地面高程h（单位），而
H
=
L
L
A
H
(
E
⃗
)
H=LLA_H({\vec E})
H=LLAH(E)返回的是当前位置下的海拔，则可视的条件是：
∀
ϱ
∈
[
0
,
1
]
:
F
d
e
m
(
E
⃗
)
<
L
L
A
H
(
E
⃗
)
\forall {\varrho \in \left [0,1\right ]}: F_{dem}(\vec E)<LLA_{H}(\vec E)
∀ϱ∈[0,1]:Fdem(E)<LLAH(E)
通俗的说，就是AB连线上的每一点都高于地形。
这种判断，需要有精确的数字高程数据的支持。函数F查询DEM阵列，根据经纬度网格，进行插值，获取高程。此类计算极为耗时，且与DEM数据的品类高度相关，接口上，采用functional对象进行泛化处理。
(2) 基准海平面算法
当采用标准海平面作为表面时，DEM数据恒为0. 该命题变为线段AB与椭球面有无交点的问题。可以使用斜率折半查找的方法，找到距离地面最近的参数 ϱ \varrho ϱ。
(3) 引入仰角约束
如果还有附带条件，如试图排除“擦边球”场景，避免线段上某个位置与DEM曲面恰好相切，则可以使用仰角约束。可以要求在路径上任意一点观测A、B两点的仰角均大于 E l El El，如5度或者10度。
要注意的是，引入的约束越多，则计算越慢。
主要接口：
/*!
* \brief observe_vis_judge 假设地球没有地形，是光滑的椭球，求取AB是不是可视
* \param llaA_in 观测点的经纬度
* \param llaB_in 被观测点经纬度
* \param minEl 仰角约束，路径上所有地表位置观察A、B的仰角均要高于minEl
* \param rad 经纬度采用弧度
* \param latfirst 纬度在第一维度[0]
* \return true 是AB可视，false是不可视
*/
inline bool observe_vis_judge(
const double llaA_in[/*3*/],
const double llaB_in[/*3*/],
const double minEl = 0,
const bool rad = false,
const bool latfirst = true
)；
/*!
* \brief observe_dem_judge 根据数字高程DEM判断可视关系
* \param llaA_in 观测点的经纬度
* \param llaB_in 被观测点经纬度
* \param dem 数字高程回调，给入纬度（度）、经度（度），返回高度（米）
* \param dem_in_meters dem数据的分辨率（米）,将使用高于1/2分辨率的精度进行搜索测试
* \param minEl 仰角约束，路径上所有地表位置观察A、B的仰角均要高于minEl
* \param rad 经纬度采用弧度
* \param latfirst 纬度在第一维度[0]
* \return true 是AB可视，false是不可视
*/
inline bool observe_dem_judge(
const double llaA_in[/*3*/],
const double llaB_in[/*3*/],
std::function<double (double lat,double lon)> dem,
const double dem_in_meters,
const double minEl = 0,
const bool rad = false,
const bool latfirst = true
)；
