0

我有一个问题一直在尝试解决，但我无法给出答案。我在 Matlab 中编写了一个函数，给定两个纬度/经度点和两个方位角，将返回两个大圆的交点。

但是，我真正需要的是沿两个初始航向的第一个大圆交点。即，如果两架飞机从纬度/经度点 1 和 2 开始，初始方位为轴承 1 和轴承 2，那么两个大圆交点中的哪一个是它们遇到的第一个？有很多解决方案（使用haversine）可以让我更接近两点，但我实际上并不关心哪个更接近，我关心的是在给定特定起点和标题的情况下我会首先遇到哪个。在很多情况下，两个路口中距离较近的实际上是遇到的第二个路口。

现在，我意识到我可以用很多条件语句来处理不同的情况，但我认为必须有一种方法来处理我采取所有交叉产品的顺序（下面给出的功能代码），但是我根本想不出正确的解决方案！我还应该提到这个函数将用于大型计算密集型模型，所以我认为这个问题的解决方案需要相当优雅/快速。谁能帮我解决这个问题？

以下不是我的函数代码（我无法在此处列出），但它是我的函数所基于的伪代码：

 %Given inputs of lat1,lon1,Bearing1,lat2,lon2,Bearing2:
%Calculate arbitrary secondary point along same initial bearing from first
%point
dAngle = 45;
lat3 = asind( sind(lat1)*cosd(dAngle) + cosd(lat1)*sind(dAngle)*cosd(Bearing1));
lon3 = lon1 + atan2( sind(Bearing1)*sind(dAngle)*cosd(lat1), cosd(dAngle)-sind(lat1)*sind(lat3) )*180/pi;
lat4 = asind( sind(lat2)*cosd(dAngle) + cosd(lat2)*sind(dAngle)*cosd(Bearing2));
lon4 = lon2 + atan2( sind(Bearing2)*sind(dAngle)*cosd(lat2), cosd(dAngle)-sind(lat2)*sind(lat4) )*180/pi;


%% Calculate unit vectors
% We now have two points defining each of the two great circles.  We need
% to calculate unit vectors from the center of the Earth to each of these
% points
[Uvec1(1),Uvec1(2),Uvec1(3)] = sph2cart(lon1*pi/180,lat1*pi/180,1);
[Uvec2(1),Uvec2(2),Uvec2(3)] = sph2cart(lon2*pi/180,lat2*pi/180,1);
[Uvec3(1),Uvec3(2),Uvec3(3)] = sph2cart(lon3*pi/180,lat3*pi/180,1);
[Uvec4(1),Uvec4(2),Uvec4(3)] = sph2cart(lon4*pi/180,lat4*pi/180,1);


%% Cross product
%Need to calculate the the "plane normals" for each of the two great
%circles
N1 = cross(Uvec1,Uvec3);
N2 = cross(Uvec2,Uvec4);


%% Plane of intersecting line
%With two plane normals, the cross prodcut defines their intersecting line
L = cross(N1,N2);
L = L./norm(L);
L2 = -L;
L2 = L2./norm(L2);


%% Convert normalized intersection line to geodetic coordinates
[lonRad,latRad,~]=cart2sph(L(1),L(2),L(3));
lonDeg = lonRad*180/pi;
latDeg = latRad*180/pi;
[lonRad,latRad,~]=cart2sph(L2(1),L2(2),L2(3));
lonDeg2 = lonRad*180/pi;
latDeg2 = latRad*180/pi;

更新： Mathworks 论坛上的一位用户指出了这一点：

实际上，他们可能每个人都达到了不同的点。我可能误解了这个问题，但你措辞的方式表明两条轨迹将收敛到同一点，这是不正确的。如果你想象第一个点在一个交叉点之后一点点，而第二个点在另一个交叉点点点之后，你就会遇到这样一种情况，每个“平面”都会朝着一开始最靠近另一个平面的交叉点移动。

我没有考虑到这一点。实际上很有可能每个平面首先与不同的大圆相交。这只是让事情变得更加复杂......

4

1 回答 1

1

大概你有一个平面方向的速度矢量（你想要寻找你的第一个交点的方向 - “地球表面上的方位矢量”）。您实际上并没有指定这一点。让我们称之为v。

您还有两点的笛卡尔坐标P1和P2，以及平面的初始位置P0。假设每个都已经在单位球面上（长度 = 1）。

现在我们想知道哪个是更短的距离 - 到 P1 或 P2 - 所以我们需要知道“从法线向量的方向看”的角度。为此，我们需要 thesin和cos- 然后我们可以使用 找到角度atan2。我们sin从叉积（记住所有向量都被归一化）和cos点积中获得。为了得到sin“从正确的方向看”，我们用法向量取点积。

这是一段将所有内容放在一起的代码-我对点和方向矢量使用了非常简单的坐标，因此我可以在脑海中确认这给出了正确的答案：

% sample points of P0...P2 and v
P0 = [1 0 0];
P1 = [0 1 0];
P2 = [0 -1 0];
v = [0 1 0];
% compute the "start direction normal":
n0 = cross(P0, v);
n0 = n0 / norm( n0 );  % unit vector 

% compute cross and dot products:
cr01 = cross(P0, P1);
cr02 = cross(P0, P2);
cos01 = dot(P0, P1);
cos02 = dot(P0, P2);

% to get sin with correct sign, take dot product with direction normal:
sin01 = dot(cr01, n0);
sin02 = dot(cr02, n0);
% note - since P0 P1 and P2 are all in the same plane
% cr02 and cr02 are either pointing in the same direction as n0
% or the opposite direction. In the latter case we get a sign flip for the sin
% in the former case this does nothing

% Now get the angle, mapped from 0 to 2 pi:
ang1 = mod(atan2(sin01, cos01), 2*pi);
ang2 = mod(atan2(sin02, cos02), 2*pi);

if( ang1 < ang2 ) 
  fprintf(1,'point 1 is reached first\n');
  % point 1 is the first one reached
else
  fprintf(1,'point 2 is reached first\n');
  % point 2 is the first
end

当您更改速度矢量的方向（指向 P2 而不是 P1）时，程序会正确地告诉您“首先到达点 2”。

让我知道这是否适合您！

于 2013-09-27T23:16:17.880 回答