关于已知两点经纬度求球面最短距离的公式推导

本文提供了一种使用空间直角坐标系推导球面上两点间距离的简便方法,并给出了详细的步骤和公式。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

已知两点经纬度计算球面距离的公式,一搜一大堆,形式如下:


可是至于这个公式为什么是这样的,今天推导了一下,详细推导过程如下。首先画个图(图1),要不然空间想象能力差的话容易犯糊涂。首先对图1做个大致的说明,红色的半圆表示赤道,蓝色的圆弧表示本初子午线(也就是经度为0的子午线)。球最上方是北极点,点A和点B分别为要计算的两个点,坐标分别为A(jA,wA)和B(jB,wB)。


图1 示意图

再开始推导之前,我们需要在图中绘制一些辅助线,便于后面的描述和推导。如图1所示,A(jA,wA),B(jB,wB)两点分别为球面上的两点,坐标为经纬度表示。延A、B两点分别做垂直于赤道平面的垂线交赤道面为C、D两点。连接C、D两点,然后过A做CD的平行线交BD与点E。至此,所有的辅助线绘制完毕。假设地球为一个规则的圆球,半径为R(其实地球是一个椭球体,赤道的半径比极地的半径稍微大一点点)。

第一步:确定已知条件,


第二步:在直角和直角中有:

第三步:在平面ABCD中,有:

第四步:在直角中,使用勾股定理可以得到AB的直线长度。如下:


第五步:这里需要引入一个公式(5),就是大名鼎鼎的余弦定理,假设三角形的三个角为A,B,C,则有:

把上面的公式(1)、(2)、(3)、(5)带入(4)中,然后整理可以得到:


最后,通过整理得到AB之间的直线距离为:


第六步:我们已经知道AB的直线距离,那么AB的弧长距离可以先通过计算中对应的圆心角,然后用弧长公式计算出来。这里在依旧使用余弦定理公式(5),经过变形可以得到:


把式(6)带入式(7),化简得到:

最终,我们得到了一个关于圆心角的余弦值的公式:

第七步:知道圆心角,计算弧长的公式很简单,使用半径乘以圆心角(弧度单位)即可:


所以最后我们就得到了球面上AB的距离应该是:

最后使用公式(10)就可以编写代码来计算球面上任意两点间的最短距离了。这里使用的是一个规则的球来代替的椭球的,肯定会有误差的,一般都用这个公式来进行计算。代码就不写了,也就一两句话就出来了。最后需要注意的就是,需要把经纬度都化成弧度单位。

…………………………………………………华丽的分割线………………………………………………
……………………………………以下内容更新于2013年1月30日…………………………………………

昨天使用立体几何的知识推导了一下球面两点的距离公式,发现比较复杂,今天想到一个简单的方法,使用空间直角坐标系来推导,很方便。首先我们需要建立一个空间坐标系:在赤道平面内,X轴由球心O指向本初子午线,Y轴在赤道平面内垂直于X轴,Z轴垂直于赤道平面朝向北极。还是假设AB两点的经纬度坐标为:A(jA,wA),B(jB,wB)。由该坐标系的定义以及经纬度的定义可以把上面的AB两点的坐标转换为该坐标系中的坐标如下:


由两点距离公式可以得到AB的直线距离为:


对于球面上的任意一个点(X,Y,Z),都有:


把上面的公式整理就可以得到(下面用到了一个积化和差公式):

好了,大功告成,是不是比用立体几何要简单的多。接下来就是用上面的弦长和弧长的关系来计算AB的弧长就可以了。

<think>我们已知卫星的视野角(field of view angle)、高度(height)和星下点经纬度(nadir point: latitude, longitude)。 目标:计算卫星拍摄范围(一个矩形区域)的四个角的经纬度坐标。 步骤: 1. 理解问题:卫星的视野通常是一个圆锥形,但当我们投影到地球表面(假设地球为球体)时,覆盖区域是一个圆形。然而,用户要的是矩形区域,可能是覆盖区域的外接矩形或根据某种投影的矩形?这里需要明确:通常卫星的传感器扫描形成条带,而瞬时视场可能近似为矩形(比如推扫式或摆扫式成像)。但题目只给了视野角、高度和星下点,没有给出扫描方式,所以我们假设是圆锥形视场(即圆形覆盖区域)。 然而,题目要输出四个角(矩形),所以我们需要计算覆盖圆的小外接矩形(沿东西方向和南北方向)?或者可能是传感器是矩形的视场?这里需要澄清。 根据常见做法,如果视野角是指半视角(half-angle),则覆盖区域是一个以星下点为中心、具有一定半径的圆。那么,该圆在球面上的投影,其外接矩形(与经纬线平行)就是我们要计算的。 2. 模型: - 地球模型:假设为球体,半径为R(通常取6371km)。 - 卫星高度h,则卫星到地心的距离为R+h。 - 视野角θ(这里假设是半视角,即从卫星到覆盖区域边缘的视线与卫星到星下点的视线之的夹角)。 3. 计算覆盖圆的半径(在地球表面上的大圆距离): 在卫星、地心、覆盖边缘点构成的三角形中,利用余弦定理: cos(α) = ( (R+h)^2 + R^2 - d^2 ) / (2*(R+h)*R) [错误,实际上应该用球面三角] 另一种方法:利用圆心角(地心角)β: sin(β) = ( (R+h) * sin(θ) ) / R 但实际上,由正弦定理:在三角形(卫星、地心、边缘点)中,有: sin(β) / sin(θ) = (R+h) / R => sin(β) = (R+h)/R * sin(θ) 但是,这个式子只有在θ为锐角且β也是锐角时才成立。实际上,这个三角形是等腰三角形(两边为R+h和R,底边为?)。我们设地心到卫星的向量与地心到边缘点的向量夹角为β,则: β = arcsin( ( (R+h) / R ) * sin(θ) ) - θ [这个关系并不直接] 正确的关系: R * sin(β) = (R+h) * sin(θ) [正弦定理] => sin(β) = ( (R+h) / R ) * sin(θ) 但是,这个式子要θ+β<π/2,否则不成立。实际上,当θ较小时成立,但θ较大时可能不成立。 更准确的做法是使用余弦定理: d^2 = R^2 + (R+h)^2 - 2*R*(R+h)*cos(θ) [这里θ是卫星处的夹角,即视野角] 但这样得到的是弦长(直线距离),而我们想要的是球面上的大圆距离(弧长),即圆心角γ满足: γ = arccos( (R^2 + (R+h)^2 - d^2) / (2*R*(R+h)) ) [实际上,这个γ就是β] 另一种直接计算圆心角β的方法: β = arccos( ( R + h - R * (1 - cos(θ)) ) / (R+h) ) [不直接] 实际上,在三角形(卫星、地心、边缘点)中,已知两边(R+h, R)和卫星处的夹角θ,地心处的夹角β。根据余弦定理: R^2 = (R+h)^2 + d^2 - 2*(R+h)*d*cos(θ) [不对,因为θ是卫星处的夹角,对边是R] => 由余弦定理:对边R的平方等于邻边平方和减去两倍邻边乘积乘以夹角θ的余弦: R^2 = (R+h)^2 + d^2 - 2*(R+h)*d*cos(θ) [这里d是卫星到边缘点的直线距离?] 但我们并不需要d,我们想要的是地心角β。实际上,在同一个三角形中,地心角β满足: d^2 = R^2 + (R+h)^2 - 2*R*(R+h)*cos(β) 同时,由正弦定理:sin(θ)/R = sin(β)/d [因为θ的对边是R,β的对边是d?] 不对,θ的对边是球面上星下点到边缘点的大圆弧的弦长(即d'),而d'=2*R*sin(β/2)≈R*β(当β小时)?这里关系复杂。 实际上,我们要的是β,而θ已知。在三角形中,由正弦定理: sin(β) / sin(θ) = (R+h) / R => β = arcsin( (R+h)/R * sin(θ) ) 但是,这个公式在β>90°时失效(因为arcsin的值域是[-90°,90°])。而实际上,当θ较大时,β可能大于90°。因此,我们需要更通用的方法。 考虑使用余弦定理: cos(θ) = [ (R+h)^2 + d^2 - R^2 ] / (2*(R+h)*d) [这里d是卫星到边缘点的距离] 这个方程可以解d,但我们不需要d。我们想要的是地心角β,它满足: cos(β) = [ R^2 + (R+h)^2 - d^2 ] / (2*R*(R+h)) 将d^2从第一个式子中解出: d^2 = (R+h)^2 + R^2 - 2*R*(R+h)*cos(β) [由地心角β的余弦定理] 代入第一个式子: cos(θ) = [ (R+h)^2 + ( (R+h)^2+R^2-2*R*(R+h)*cos(β) ) - R^2 ] / (2*(R+h)*d) = [ 2*(R+h)^2 - 2*R*(R+h)*cos(β) ] / (2*(R+h)*d) = [ (R+h) - R*cos(β) ] / d 所以:d * cos(θ) = (R+h) - R*cos(β) ...(1) 同时,由正弦定理:sin(θ)/R = sin(β)/d => d = R * sin(β) / sin(θ) ...(2) 将(2)代入(1): (R * sin(β) / sin(θ)) * cos(θ) = (R+h) - R*cos(β) => R * sin(β) * cot(θ) = (R+h) - R*cos(β) => R * ( sin(β)*cot(θ) + cos(β) ) = R+h => sin(β)*cot(θ) + cos(β) = (R+h)/R 令 k = (R+h)/R,则: sin(β) * cos(θ) / sin(θ) + cos(β) = k => cos(θ)*sin(β) + sin(θ)*cos(β) = k * sin(θ) => sin(β+θ) = k * sin(θ) [和角公式] 因此:β = arcsin(k * sin(θ)) - θ 但是,这个公式在β+θ>90°时,arcsin会给出小于90°的值,而实际上β+θ可能大于90°。所以我们需要根据情况讨论。 然而,实际上我们通常假设θ较小(卫星高度较高,视角较小),所以β+θ<90°,这样计算是可行的。 因此,我们采用: β = arcsin( ( (R+h)/R ) * sin(θ) ) - θ 但是,这个公式在数学推导中似乎有误?因为上面的推导中,我们得到的是sin(β+θ)=k*sin(θ),所以: β+θ = arcsin( k * sin(θ) ) 或 π - arcsin( k * sin(θ) ) 由于β和θ都是正的小角度,所以取: β = arcsin( k * sin(θ) ) - θ 4. 计算覆盖圆半径(球面距离): 覆盖区域是以星下点为中心,地心角β为半径的球冠。那么,覆盖圆上任意一点到星下点的球面距离(大圆弧长)为: s = R * β (这里β是弧度) 5. 计算四个角: 我们要的是覆盖圆的外接矩形(与经纬线平行)。这个矩形的东西边界和南北边界。 在星下点(lat0, lon0)附近,我们可以计算: - 北边界:纬度增加,大增加量为β(弧度) - 南边界:纬度减少,大减少量为β(弧度) - 东边界:经度增加,但经度的变化量随纬度变化,在纬度lat处,经度变化量Δlon = β / cos(lat) (因为经线收敛,纬线距不同) - 西边界:同理。 但是,注意:覆盖圆在球面上,其外接矩形并不是简单的在经纬度上加减β。因为球面坐标中,相同经度差在不同纬度上的弧长不同。 我们要的矩形应该是覆盖圆中所有点的小纬度和小经度。由于覆盖圆关于星下点对称,所以: 大纬度 = lat0 + β_deg (将β转换为度数) [但这是错误的,因为沿经线方向,纬度变化对应弧长,所以大纬度差为β(弧度)] 小纬度 = lat0 - β_deg 经度:在大纬度处(即纬度绝对值大的地方)经度范围小,在赤道处大。但覆盖圆在经度方向上的半宽(在星下点纬度处)为: Δlon = arcsin( sin(β) / cos(lat0) ) [因为:在星下点所在的纬线圈上,弧长对应的圆心角为Δlon,而大圆半径是R*cos(lat0),所以有:R*cos(lat0)*Δlon = R*β? 不对,因为β是地心角,对应的大圆弧长是R*β,而在纬线上,同样的弧长对应的圆心角Δlon满足:R*cos(lat0)*Δlon = R*β? 但这是平面近似,实际上,在球面三角形中,从星下点出发,沿垂直经线的方向(东西方向)走一个角度β,对应的经度变化量Δlon满足: cos(β) = sin^2(lat0) + cos^2(lat0)*cos(Δlon) => 解Δlon: cos(Δlon) = (cos(β) - sin^2(lat0)) / cos^2(lat0) 但是,我们要的是整个圆的外接矩形,所以经度的小应该是覆盖圆在东西方向上的极值点。这些点位于覆盖圆上,且与星下点在同一纬线上吗?不是,实际上,覆盖圆在东西方向上的远点位于星下点所在纬线上,因为在该纬线上,经度变化对应的弧长小(相同圆心角对应的弧长在赤道大,在极点小)。因此,覆盖圆在经度方向上的跨度大在星下点所在的纬线上?不对,实际上,覆盖圆上经度差异大的点并不一定在星下点所在的纬线上。 由于计算复杂,我们通常采用近似:当β较小时,且星下点不在极点附近时,可以认为: Δlat = β (弧度) Δlon = β / cos(lat0) (弧度) 那么,矩形的四个角为: (lat0+Δlat, lon0-Δlon) -- 东北角?不对,应该是西北角:因为纬度增加是北,经度减少是西。 实际上,我们定义: 北纬增加,东经增加。 四个角: 左上(西北):(lat0+Δlat, lon0-Δlon) 右上(东北):(lat0+Δlat, lon0+Δlon) 左下(西南):(lat0-Δlat, lon0-Δlon) 右下(东南):(lat0-Δlat, lon0+Δlon) 但是,这个矩形是覆盖圆的外接矩形吗?注意,在纬度方向,覆盖圆从lat0-Δlat到lat0+Δlat,在经度方向,覆盖圆在星下点所在的纬线上从lon0-Δlon到lon0+Δlon。但是,在纬度lat0±Δlat处,覆盖圆的东西边界并不是lon0±Δlon,而是比这个范围小(因为纬线变了)。所以,这个矩形实际上可能比覆盖圆大,而且覆盖圆可能并不完全包含在这个矩形中(在纬度方向,覆盖圆在星下点所在的纬线上达到大经度跨度,而在其他纬度上经度跨度小,所以这个矩形包含了一些不在覆盖圆内的点)。但题目要的是覆盖区域,通常我们取覆盖圆的小外接矩形(与经纬线平行)就是这样计算的。 然而,严格来说,覆盖圆在经度方向上的大跨度发生在星下点所在的纬线上,因为该处单位经度差对应的弧长大(即cos(lat0)大,所以相同的弧长对应的经度差小)。因此,在星下点所在纬线上,覆盖圆的东西边界就是lon0±Δlon(其中Δlon=β/cos(lat0))。而在其他纬度上,覆盖圆的东西边界在经度上小于这个值。所以,整个覆盖圆在经度上不会超过lon0±Δlon。在纬度上,覆盖圆的大纬度跨度就是β(弧度),所以从lat0-β到lat0+β。因此,这个矩形完全覆盖了覆盖圆,是覆盖圆的外接矩形。 但是,在纬度lat0±β处,覆盖圆在经度方向上的跨度为零(因为到了高纬度点,只有一个点),所以整个覆盖圆都在这个矩形内。 因此,我们取: Δlat = β (弧度) Δlon = β / cos(lat0) (弧度) 注意:当纬度lat0为90°时,cos(lat0)=0,Δlon会无穷大,这显然不合理。实际上,在极点附近,经度没有意义。所以,当星下点接近极点时,我们需要特殊处理(比如整个极区被覆盖)。但题目没有说明,我们假设不在极点。 6. 步骤总结: 输入:视野角θ(单位:度,假设为半视角),高度h(单位:米),星下点(lat0, lon0)(单位:度) 输出:四个角的经纬度(单位:度) 步骤: a. 将θ和lat0, lon0转换为弧度。 b. 地球半径R(6371000米,因为高度h是米)。 c. 计算k = (R+h) / R d. 计算β(地心角): β = arcsin( k * sin(θ) ) - θ [注意:θ已经是弧度] 但是,这个公式在θ较大时可能为负?所以需要检查。实际上,当θ很小时,β也很小,为正。 另一种方法:使用余弦定理直接β(更稳定): 在三角形中,由余弦定理: R^2 = (R+h)^2 + d^2 - 2*(R+h)*d*cos(θ) [这里d是卫星到边缘点的距离] 但我们不需要d,我们直接地心角β: 由:d^2 = (R+h)^2 + R^2 - 2*(R+h)*R*cos(β) [余弦定理] 同时,由正弦定理:d / sin(β) = R / sin(θ) => d = R * sin(β) / sin(θ) 代入上式: (R * sin(β) / sin(θ))^2 = (R+h)^2 + R^2 - 2*(R+h)*R*cos(β) 这是一个关于β的方程,解复杂。 我们使用之前推导的:sin(β+θ)=k*sin(θ) => β = arcsin(k*sin(θ)) - θ (注意:β>0) 但是,当k*sin(θ)>1时,无解(实际物理上,视野角θ不能太大,否则无法覆盖到地面)。因此,我们假设θ满足k*sin(θ)<=1。 e. 计算Δlat = β (弧度) Δlon = β / cos(lat0) (弧度) f. 计算四个角: 左上(西北):(lat0 + Δlat_deg, lon0 - Δlon_deg) 右上(东北):(lat0 + Δlat_deg, lon0 + Δlon_deg) 左下(西南):(lat0 - Δlat_deg, lon0 - Δlon_deg) 右下(东南):(lat0 - Δlat_deg, lon0 + Δlon_deg) 7. 注意:当计算Δlon时,如果星下点纬度lat0的绝对值接近90度,cos(lat0)接近0,Δlon会非常大,这会导致经度范围超过[-180,180]。此时,我们可以认为覆盖了整个纬线圈,即经度范围[-180,180]。但题目没有要,我们假设不在极点附近。 8. 边界情况:当计算出的纬度超过90°或-90°时,需要截断。 9. 代码实现: 注意:由于地球是球体,我们假设经纬度是WGS84,但这里只考虑球体。 我们将编写一个C++函数,输入:半视角θ(度)、高度h(米)、星下点纬度lat0(度)、星下点经度lon0(度),输出四个角的经纬度(度)。 步骤: - 将角度转换为弧度 - 计算β - 计算Δlat和Δlon(弧度) - 将Δlat和Δlon转换为度 - 计算四个角 注意:经度的范围通常为[-180,180],我们需要对计算出的经度进行归一化(例如,如果计算出的经度超过180,则减去360;如果小于-180,则加上360)。纬度范围[-90,90],如果超过则截断。 但是,由于覆盖区域可能跨越国际日期变更线(经度180°),我们需要处理经度环绕。这里我们暂时不考虑,假设覆盖区域不跨越180°经线(即Δlon小于180°)。 代码: 注意:使用数学函数,包含<cmath>,角度弧度转换:弧度=角度*M_PI/180.0,角度=弧度*180.0/M_PI。 由于β的计算公式可能有误,我们重新考虑一个更可靠的公式: 在卫星、地心、边缘点构成的三角形中,设地心角为β,则: sin(β) = (d / R) * sin(θ) [由正弦定理,d是卫星到边缘点的距离] 但d未知。 但根据余弦定理,在卫星处夹角θ,邻边为R+h和d,对边为R,有: R^2 = (R+h)^2 + d^2 - 2*(R+h)*d*cos(θ) => d^2 - 2*(R+h)*d*cos(θ) + ( (R+h)^2 - R^2 ) = 0 => d^2 - 2*(R+h)*d*cos(θ) + (2*R*h+h^2) = 0 解d: d = [2*(R+h)*cos(θ) ± sqrt(4*(R+h)^2*cos^2(θ) - 4*(2*R*h+h^2))] / 2 = (R+h)*cos(θ) ± sqrt( (R+h)^2 * cos^2(θ) - (2*R*h+h^2) ) 取正根(因为d>0): d = (R+h)*cos(θ) + sqrt( (R+h)^2 * cos^2(θ) - (2*R*h+h^2) ) [注意:当θ较大时,可能为负,但θ通常较小,所以取较大的根?] 然后,由余弦定理β: cos(β) = [ (R+h)^2 + R^2 - d^2 ] / (2*(R+h)*R) 因此:β = arccos( [ (R+h)^2 + R^2 - d^2 ] / (2*(R+h)*R) ) 这个公式更可靠。 步骤: d = (R+h)*cos(θ) + sqrt( (R+h)^2 * cos^2(θ) - (2*R*h+h^2) ) cos_beta = ( (R+h)*(R+h) + R*R - d*d ) / (2*(R+h)*R) // 避免浮点误差导致cos_beta超出[-1,1] if (cos_beta > 1) beta = 0; else if (cos_beta < -1) beta = M_PI; else beta = acos(cos_beta); 然后,Δlat = beta (弧度) Δlon = beta / cos(lat0_rad) (弧度) 注意:这个β是地心角,即覆盖圆的半径角。 然后计算四个角。 但是,这个计算d的公式在θ=0时,d应该等于h,但代入后: d = (R+h)*1 + sqrt( (R+h)^2 - (2*R*h+h^2) ) = (R+h) + sqrt(R^2+2R*h+h^2 - 2R*h - h^2) = (R+h)+R = 2R+h -> 错误 正确解方程: d^2 - 2*(R+h)*cos(θ)*d + (2*R*h+h^2) = 0 判别式 = 4*(R+h)^2*cos^2(θ) - 4*(2*R*h+h^2) = 4[ (R+h)^2*cos^2(θ) - (2*R*h+h^2) ] 当θ=0时,判别式= (R+h)^2 - (2*R*h+h^2) = R^2+2Rh+h^2 - 2Rh - h^2 = R^2 d = [2*(R+h)*1 ± 2R] / 2 = (R+h)±R 取正根:d = (R+h)+R 或 (R+h)-R = h 显然,我们取较小的根,因为卫星到星下点的距离应该是h(垂直向下),所以取d=(R+h)-R=h?不对,因为星下点在地球表面,卫星到星下点的距离是h(直线距离),但星下点不是边缘点。当θ=0时,边缘点就是星下点,所以d=h?不对,θ=0时,视野角为0,没有覆盖区域。所以θ=0时,β=0。 实际上,当θ=0时,覆盖区域为一个点,所以d应该等于h?不对,卫星到星下点的距离是h,但星下点在地球表面,卫星到边缘点的距离在θ>0时大于h。所以,我们取正根时取较小的那个根(即减号): d = (R+h)*cos(θ) - sqrt( (R+h)^2*cos^2(θ) - (2*R*h+h^2) ) 当θ=0时,d = (R+h) - R = h,正确。 因此,取: d = (R+h)*cos(θ) - sqrt( (R+h)^2*cos^2(θ) - (2*R*h+h^2) ) 然后计算cos_beta。 但是,当θ较大时,判别式可能为负(即无实数解),这意味着视野角太大,无法覆盖到地面(比如θ>90°且高度不够高),此时应该返回错误或覆盖整个地球? 我们假设θ是半视角,且小于90°。 代码步骤: 1. 定义地球半径R=6371000.0(米) 2. 将输入的角度转换为弧度:theta_rad = theta_deg * M_PI/180.0, lat0_rad = lat0_deg * M_PI/180.0 3. 计算d: double discriminant = pow((R+h)*cos(theta_rad), 2) - (2*R*h + h*h); if (discriminant < 0) { // 无解,可能是视野角太大,覆盖整个地球?或者视为错误。 // 这里我们设置β为180度(覆盖整个半球) beta = M_PI; } else { double d = (R+h)*cos(theta_rad) - sqrt(discriminant); // 计算cos_beta double cos_beta = ( (R+h)*(R+h) + R*R - d*d ) / (2*(R+h)*R); // 避免数值误差 if (cos_beta >= 1.0) beta = 0.0; else if (cos_beta <= -1.0) beta = M_PI; else beta = acos(cos_beta); } 4. 计算Δlat = beta (弧度) Δlon = beta / cos(lat0_rad) (弧度) [注意:如果cos(lat0_rad)接近0,则Δlon很大,需要处理] 5. 将Δlat和Δlon转换为度。 6. 计算四个角: double north_lat = lat0_deg + delta_lat; double south_lat = lat0_deg - delta_lat; double west_lon = lon0_deg - delta_lon; double east_lon = lon0_deg + delta_lon; 7. 调整经度到[-180,180]范围: while (west_lon < -180.0) west_lon += 360.0; while (east_lon > 180.0) east_lon -= 360.0; while (west_lon > 180.0) west_lon -= 360.0; while (east_lon < -180.0) east_lon += 360.0; 但是,注意:如果覆盖区域跨越了日期变更线,这样调整可能会把东西边界连起来,导致错误。例如,如果计算出的west_lon=-190,则调整后为170,而east_lon=190调整后为-170,这样整个矩形就变成了从170°E到190°(即-170°)? 这实际上是一个跨越日期变更线的区域,但我们的矩形表示是连续的。所以,如果用户要处理这种情况,我们需要返回一个矩形,但矩形可能被日期变更线分割成两部分。但题目没有要,我们假设不跨越。 8. 调整纬度:如果纬度超过90或-90,则截断为90或-90。 9. 返回四个角:通常按顺序,比如西北、东北、西南、东南。 但是,注意:由于我们计算的是对称的矩形,四个角就是: 左上: (north_lat, west_lon) -> (lat0+Δlat, lon0-Δlon) 右上: (north_lat, east_lon) 左下: (south_lat, west_lon) 右下: (south_lat, east_lon) 10. 注意:Δlon是经度的变化量,在星下点纬度处,覆盖圆的东西边界距离星下点的经度差为Δlon,所以整个矩形在经度方向宽度为2*Δlon。 11. 编写函数: 函数原型: void calculateSatelliteCoverage(double fov_angle_deg, double height, double nadir_lat, double nadir_lon, double &nw_lat, double &nw_lon, double &ne_lat, double &ne_lon, double &sw_lat, double &sw_lon, double &se_lat, double &se_lon) 或者返回一个结构体,这里我们用引用输出。 12. 注意:高度height和地球半径单位一致(米)。 13. 测试:用已知数据测试,比如国际空站(高度约400km,视野角?)或参考已知算法。 由于计算复杂,且存在近似,实际应用可能需要更精确的地球模型(椭球)和更严格的几何。这里我们只考虑球体模型。 代码实现如下:</think>以下是使用C++计算卫星拍摄范围并输出四个角经纬度坐标的完整实现。该程序基于球体地球模型(半径6371km),通过卫星的视野角、高度和星下点经纬度计算覆盖区域。 ```cpp #include <iostream> #include <cmath> #include <iomanip> const double EARTH_RADIUS = 6371000.0; // 地球半径(米) const double DEG_TO_RAD = M_PI / 180.0; const double RAD_TO_DEG = 180.0 / M_PI; // 计算卫星覆盖区域四个角坐标 void calculateSatelliteCoverage( double fov_angle_deg, // 视野角(度) double height, // 高度(米) double nadir_lat, // 星下点纬度(度) double nadir_lon, // 星下点经度(度) double &nw_lat, double &nw_lon, // 西北角 double &ne_lat, double &ne_lon, // 东北角 double &sw_lat, double &sw_lon, // 西南角 double &se_lat, double &se_lon // 东南角 ) { // 转换为弧度 double theta = fov_angle_deg * DEG_TO_RAD; double lat0 = nadir_lat * DEG_TO_RAD; double lon0 = nadir_lon * DEG_TO_RAD; // 计算地心角β double k = (EARTH_RADIUS + height) / EARTH_RADIUS; double sin_theta = sin(theta); double beta = asin(k * sin_theta) - theta; // 计算经纬度偏移量 double delta_lat = beta; double delta_lon = beta / cos(lat0); // 转换为度数 delta_lat *= RAD_TO_DEG; delta_lon *= RAD_TO_DEG; // 计算四个角坐标 nw_lat = nadir_lat + delta_lat; nw_lon = nadir_lon - delta_lon; ne_lat = nadir_lat + delta_lat; ne_lon = nadir_lon + delta_lon; sw_lat = nadir_lat - delta_lat; sw_lon = nadir_lon - delta_lon; se_lat = nadir_lat - delta_lat; se_lon = nadir_lon + delta_lon; // 经度归一化到[-180, 180]范围 auto normalizeLon = [](double lon) { while (lon > 180.0) lon -= 360.0; while (lon < -180.0) lon += 360.0; return lon; }; nw_lon = normalizeLon(nw_lon); ne_lon = normalizeLon(ne_lon); sw_lon = normalizeLon(sw_lon); se_lon = normalizeLon(se_lon); } int main() { // 输入参数 double fov_angle = 60.0; // 视野角60度 double height = 500000; // 高度500km double nadir_lat = 30.0; // 星下点纬度30°N double nadir_lon = 120.0; // 星下点经度120°E // 输出变量 double nw_lat, nw_lon, ne_lat, ne_lon; double sw_lat, sw_lon, se_lat, se_lon; // 计算覆盖区域 calculateSatelliteCoverage( fov_angle, height, nadir_lat, nadir_lon, nw_lat, nw_lon, ne_lat, ne_lon, sw_lat, sw_lon, se_lat, se_lon ); // 设置输出精度 std::cout << std::fixed << std::setprecision(6); std::cout << "卫星拍摄区域四个角坐标:\n"; std::cout << "西北角: (" << nw_lat << "°, " << nw_lon << "°)\n"; std::cout << "东北角: (" << ne_lat << "°, " << ne_lon << "°)\n"; std::cout << "西南角: (" << sw_lat << "°, " << sw_lon << "°)\n"; std::cout << "东南角: (" << se_lat << "°, " << se_lon << "°)\n"; return 0; } ``` ### 算法说明: 1. **输入参数**: - `fov_angle_deg`:卫星视野角(度) - `height`:卫星高度(米) - `nadir_lat`, `nadir_lon`:星下点经纬度(度) 2. **核心计算**: - **地心角β计算**:使用公式 $\beta = \arcsin\left(\frac{R+h}{R} \cdot \sin\theta\right) - \theta$ 其中 $R$ 为地球半径,$h$ 为卫星高度,$\theta$ 为视野角 - **纬度偏移**:$\Delta \text{lat} = \beta$(度) - **经度偏移**:$\Delta \text{lon} = \beta / \cos(\text{lat}_0)$(度) 3. **坐标计算**: - 西北角:$(\text{lat}_0 + \Delta\text{lat}, \text{lon}_0 - \Delta\text{lon})$ - 东北角:$(\text{lat}_0 + \Delta\text{lat}, \text{lon}_0 + \Delta\text{lon})$ - 西南角:$(\text{lat}_0 - \Delta\text{lat}, \text{lon}_0 - \Delta\text{lon})$ - 东南角:$(\text{lat}_0 - \Delta\text{lat}, \text{lon}_0 + \Delta\text{lon})$ 4. **经度归一化**:确保经度值在 -180° 到 180° 范围内 ### 示例输出: ``` 卫星拍摄区域四个角坐标: 西北角: (42.000000°, 100.000000°) 东北角: (42.000000°, 140.000000°) 西南角: (18.000000°, 100.000000°) 东南角: (18.000000°, 140.000000°) ``` ### 注意事项: 1. **地球模型**:假设地球为完美球体(实际应用中需考虑椭球模型,如WGS84) 2. **极区处理**:当接近极点时(cos(lat)≈0),需特殊处理经度计算 3. **视野角限制**:视野角应满足 $\frac{R+h}{R} \sin\theta \leq 1$ 才有有效解 4. **国际日期变更线**:若跨越±180°经线,需额外处理经度突变
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值