Best Fit of Structural Data


We assume in all these calculations that our direction data has been converted into Cartesian Coordinate form. The code examples show the basic operations and don't contain error traps to prevent errors like division by zero and so on. Also they don't show initialization: variables should be set to zero at the beginning of any routine.

The best-fit calculations involve simple multivariate calculus. The fundamental principle in multivariate calculus is to treat one variable at a time, regard everything else as constant, and integrate or differentiate exactly as in ordinary calculus.

Best Fit of a Line to a Set of Direction Data

The best fit of a line to a set of directions is easy. If we have direction data (x1, y1, z1), (x2, y2, z2), etc. then calculate x as the average of x1, x2 .... , y as the average of y1, y2 .... and z as the average of z1, z2 .... Normalize x, y, and z by calculating s = x2 + y2 + z2 and dividing x, y, and z by s.

If we think of all the direction data as a set of vectors, the best fit vector is just the average of all the direction vectors.

Best Fit of a Plane (Great Circle) to a Set of Direction Data

This problem has many applications in geology. The best-known example is the fit of a fold axis to a collection of bedding-plane attitudes. Another application is finding the pole of rotation for the opening of an ocean basin, given matching points on opposite sides. We have a set of direction cosines, C1, C2, and C3 for each line. The data could be input by the user, or READ from a set of DATA statements. We want to find either a plane that is as nearly parallel to all the direction cosines as possible or, what is really the same thing, a line as nearly perpendicular to all the direction cosines as possible. In either case, the condition we want to satisfy as best we can is:

A * C1 + B * C2 + C3 = 0

In this equation, the direction numbers of the desired plane, or its pole, are A, B, and 1. Since the equation sums to zero, we lose no generality by assuming one of the direction numbers to be equal to 1 and it simplifies the math considerably by eliminating one variable.

In reality, the expression above will not sum exactly to zero but to some small value E. We want the sum of the squares of all the values of E to be a minimum. We have:

To evaluate the sum of the errors squared we will need the sums of C1*C1, C1*C2, and so on for all the directions in our data set. Let us define a small array S such that S(1,2) equals the sum of C1 * C2, and so on. We can write a short routine as follows:

DIM S(3,3)
WHILE C1 < 9999
S(1,1) = S(1,1) + C1 * C1
S(1,2) = S(1,2) + C1 * C2
S(1,3) = S(1,3) + C1 * C3
S(2,2) = S(2,2) + C2 * C2
S(2,3) = S(2,3) + C2 * C3
S(3,3) = S(3,3) + C3 * C3

In this example, note the use of a dummy variable to exit the input routine. The dummy variable was chosen because no real direction cosine would ever be greater than 1, and the variable number of data inputs precludes a loop of fixed length.

The condition we need to attain is the following:

= A*A*S(1,1) + B*B*S(2,2) + S(3,3) + 2*A*B*S(1,2) + 2*A*S(1,3) + 2*B*S(2,3)

In this case, we don't know A and B, so we treat them as the variables and differentiate the above formula with respect to A and B. The minimum for each derivative, of course, occurs when the derivative is zero. We obtain the two conditions:

We can then solve these two equations for A and B. The following routine will do the job.

A = S(1,2)*S(2,3) - S(1,3)*S(2,2)
A = A/(S(1,1)*S(2,2) - S(1,2)*S(1,2))
B = S(1,2)*S(1,3) - S(2,3)*S(1,1)
B = B/(S(1,1)*S(2,2) - S(1,2)*S(1,2))
C = 1/SQR(1 + A*A + B*B)
A = A*C: B = B*C

We could have solved for the third direction cosine, C, when we derived the least-squares fit, but then we would have had a messy job of solving three simultaneous equations. This approach is a lot simpler. The next-to-last step normalizes A and B so that the squares of A, B and C sum to 1

Best Fit of a Cone (Small Circle) to a Set of Direction Data

A circular cone is a useful approximation to a variety of geologic phenomena. Many folds are conical rather than cylindrical. Shatter cones, or fracture surfaces probably produced by shock waves during meteor impacts, have surface striations which define the surface of a cone. Conical distributions, when plotted on a stereonet, lie on a small circle. An application to plate tectonics might be finding the best fit pole of rotation between two plates given a transform fault (which is ideally an arc of a small circle).

Assume the cone axis has direction cosines A1, A2, and A3 and a given line we wish to fit to the cone has direction cosines C1, C2, and C3. The angle between the cone axis and any line on the surface of the cone is constant and given by:

A1 * C1 + A2 * C2 + A3 * C3 = K

where K is the cosine of half the apical angle of the cone. We can simplify the equation by letting A = A1/A3, B = A2/A3, and C = -K/A3, and the equation becomes

A * C1 + B * C2 + C3 - C = 0

As before, the formula will generally not be exactly equal to zero but will have a small value E, and we find the best fit by reducing the sum of the squares of all E to a minimum. We have

We will need the sums of C1*C1, C1*C2, and so on, but in addition we will need the sums of C1, C2 and C3. The routine we used in the previous section will work just as well here, but with modifications because we need additional sums.

DIM S(3,3)
WHILE IF C1 < 9999
S1 = S1 + C1: REM SUM OF C1
S2 = S2 + C2: REM SUM OF C2
S3 = S3 + C3: REM SUM OF C3

We follow basically the same approach as before. However, we have three independent variables and the math is correspondingly more complex. We want to minimize

SUM OF (A*C1 + B*C2 + C3 + C) * (A*C1 + B*C2 + C3 + C) =

Again, we differentiate with respect to each coefficient A, B and C, and obtain:

We now have three equations in three unknowns, and the easiest way to program the solution is through determinants:

D = S(1,1) * (S(2,2)*N - S2*S2)
D = D + S(1,2) * (S1*S2 - N*S(1,2))
D = D + S1 * (S(1,2)*N - S1*S2)

D1 = -S(1,3) * (S(2,2)*N - S2*S2)
D1 = D1 - S(2,3) * (S1*S2 - N*S(1,2))
D1 = D1 - S3 * (S(1,2)*N - S1*S2)

D2 = S(1,1) * (-S(2,3)*N + S2*S3)
D2 = D2 + S(1,2) * (S(1,3)*N - S1*S3)
D2 = D2 + S1 * (S(2,3)*S1 - S2*S(1,3))

D3 = S(1,1) * (S(2,3)*S2 - S(2,2)*S3)
D3 = D3 + S(1,2) * (S(1,2)*S3 -S(1,3)*S2)
D3 = D3 + S1 * (S(2,2)*S(1,3)-S(1,2)*S(2,3))

A = D1/D: B = D2/D: C = D3/D
A3 = 1/SQR(1 + A*A + B*B)
A1 = A*A3: A2 = B*A3: K = -C*A3

Return to Course Syllabus
Return to Course Notes Index
Return to Professor Dutch's Home Page

Created 1 April 1999, Last Update 30 January 2012

Not an official UW Green Bay site