Controlled Flight Into Terrain (CFIT)

I want to build a collision avoidance system

How can I check for whether an aircraft with a given heading will crash into terrain?

How can I check whether terrain is on the right hand side of the aircraft?

The terrain data is given in a Matrix.

[Rik]

Original question:

Controlled Flight Into Terrain (CFIT)

I want to build a collision avoidance system

How can I check for whether an aircraft with a given heading will crash into terrain?

How can I check whether terrain is on the right hand side of the aircraft?

The terrain data is given in a Matrix.

[Rik]

@Manuel Regarding your flag ("I have rewritten the question how it should have been formulated. The accepted answer is way too detailed and already demonstrates the whole program - leaving nothing to be done besides the actual thing I asked for. For people that want to also do the work, they will always be judged by the answer given. Sam's answer is what I was looking for and William did too much. I would like to delete my question."):

Deleting your question is no longer possible after an answer has been posted (unless and until all answers are deleted). You got the help, why shouldn't future readers with a similar question?

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[William Rose]

@Manuel,

Suppose you have an array whose values are the terrain elevations. For example:

Elev=1000*membrane(1,300);     %sample terrain matrix (height in m)
Elev=Elev-min(min(Elev));      %adjust terrain height so all points are >=0

Suppose you have vectors for the X coordinates (columns) and Y coordinates (rows) corresponding to each grid point.

% Assume grid is square with 0.5 nautical mile long sides.
% 1 degree of latitude = 60 n.m.; 1 degree of longitude = 60*cos(Lat) n.m.
[ny,nx]=size(Elev);
Lat0=40; Long0=-75;                     %coordinates of array element(1,1)
Lat=Lat0+(0:ny-1)/120;                  %latitude in degrees of the terrain array
Long=Long0+(0:nx-1)/(120*cosd(Lat0));   %longitude in degrees

Supose you have the aircraft current position, heading, and ground speed. Assume level flight.

p0=[-75.0,45.0,500];    %initial position (Long, Lat) & altitude (m above MSL)
gs=300;                 %ground speed (knots=n.m./hour)
hdg=135;                %compass heading (true)

Compute the predicted path.

theta=90-hdg;           %heading (degrees) measured the usual mathematical way
v=[gs*cosd(theta), gs*sind(theta), 0];  %velocity vector (knots)
t=[0:1/60:1.5]';        %time vector (hours), 1 minute spacing
vd=v.*[1/(60*cosd(Lat0)),1/60,1]; %convert velocity from nm/hr to [degrees Long/hr,degreesLat/hr,]
ppred=p0+t*vd;          %predicted position at 1 minute intervals

Plot the terrain and the predicted path.

%% Plot terrain and current and predicted position
surf(Long,Lat,Elev,'EdgeColor','none');         %terrain
xlabel('Longitude'); ylabel('Latitude'); zlabel('Elevation (m)'); 
grid on; hold on
plot3(p0(1),p0(2),p0(3),'r*');                  %current position
plot3(ppred(:,1),ppred(:,2),ppred(:,3),'-r');   %predicted position

That does not look good.

Since this is probably a homework project, the rest will be just ideas, and you can write the code. Now that you have all the information organized in a nice way, you can think of a method to determine if the predicted path intersects the terrain, and when and where that occurs.

Write some code. If it has an error or gives an unexpected results, post it on this site, with a description of the method you are tyring to use, and what the error message or unexpected result is, and what steps you have tried to fix the error.

You can calculate the slope of the terrain perpendicular to the flight direction, at the impact point (numerical estimate of the directional derivative). If the sign of the slope is one way or the other, then you know the terrain is higher on the right, or higher on the left, at the impact point. Alternatively, if there is a collision predicted, you could try a heading that is 1 degree to the left of the original heading, and try a heading 1 degree to the right of the original heading. Determine the time to impact for these alternate headings. Depending on whether the impact occurs sooner, later, or not at all with the alternate headings, you can conclude that the terrain is on the right, or on the left, or straight ahead.

Good luck!

[Sam Chak]

Wow, @William Rose, this is good stuff. 👍

Assuming that the aircraft system has the ability to generate the contour of the 'Membrane' Hill, can we use the info from the contour line for path planning, like detouring around the 'safe' contour line?

[William Rose]

@Sam Chak, thank you very much!

Great idea. The code below adds contour lines to the plot, and it returns the coordinates of the contour at the current flight level in array M. One could devise a route that avoids the area enclosed by M. You could modify the line that returns M, so that it returns the contour for (flight level-X), where X is a vertical margin of safety.

Elev=1000*(membrane(1,300));     %sample terrain matrix (height in m)
Elev=Elev-min(min(Elev));      %adjust terrain height so all points are >=0
% Assume grid is square with 0.5 nautical mile long sides.
% 1 degree of latitude = 60 n.m.; 1 degree of longitude = 60*cos(Lat) n.m.
[ny,nx]=size(Elev);
Lat0=40; Long0=-75;                     %coordinates of array element(1,1)
Lat=Lat0+(0:ny-1)/120;                  %latitude in degrees of the terrain array
Long=Long0+(0:nx-1)/(120*cosd(Lat0));   %longitude in degrees

%% Aircraft current position, heading, and ground speed. Assume level flight. 
p0=[-75.0,45.0,500];    %initial position (Long, Lat) & altitude (m above MSL)
gs=300;                 %ground speed (knots=n.m./hour)
hdg=135;                %compass heading (true)

%% Compute the predicted path.
theta=90-hdg;           %heading (degrees) measured the usual mathematical way
v=[gs*cosd(theta), gs*sind(theta), 0];  %velocity vector (knots)
t=[0:1/60:1.5]';        %time vector (hours), 1 minute spacing
vd=v.*[1/(60*cosd(Lat0)),1/60,1]; %convert velocity from nm/hr to [degrees Long/hr,degreesLat/hr,]
ppred=p0+t*vd;          %predicted position at 1 minute intervals

%% Plot terrain and current and predicted position
surf(Long,Lat,Elev,'EdgeColor','none');         %terrain
grid on; hold on
contour3(Long,Lat,Elev, '-k')                   %add contour lines
M=contour3(Long,Lat,Elev,[p0(3) p0(3)],'--r');  %contour at current flight level
plot3(p0(1),p0(2),p0(3),'r*');                  %current position
plot3(ppred(:,1),ppred(:,2),ppred(:,3),'-r');   %predicted position
xlabel('Longitude'); ylabel('Latitude'); zlabel('Elevation (m)')

@Manuel, maybe you can implement Sam's great idea.

[Manuel Dullnig]

This will work. thanks for your help

I will let you know how it goes :)

[William Rose]

@Manuel, you're welcome. Good luck!

[Sam Chak]

@Manuel, Good luck and don't forget to acknowledge or credit @William Rose's contribution in your project.

[William Rose]

Thank you, @Sam Chak.

[William Rose]

@Manuel, that is a good quesitn and it is tue tht it is not a trivial question.

If you do the evaluation as I suggested, which is to compare the (plane) altitude to the (terrain) elevaiton at each minute, then you should have time i, when altitude>elevation, and time i+1 when altitude<elevation. And you should have the exact locaiton of hte aircraft at i and i+1 (where i corresponds to time in minutes). Then you do linear interpolation to estimate when altitude=elevation.

Example:

time altitude elevation

t1=95 a1=800 e1=790

t2=96 a2=800 e2=815

Then t_CFIT = t1+ (t2-t1)*(a2-e1)/(e2-e1)

Once you know the exact time of CFIT, you can use the initial position and the velcocity vector (bot of which you know) to compute the location at the CFIT time.