Vehicular traveling direction measuring system
BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     This invention relates generally to a system for measuring a traveling direction of a vehicle using a geomagnetic direction sensor and, more specifically, to a vehicular traveling direction measuring system which is capable of providing an accurate measurement of a traveling direction of the vehicle even when a vehicle body is magnetized. 
     2. Description of the Background Art 
     One of previously proposed systems for measuring a vehicle traveling direction is exemplified by, such as, Japanese First Patent Publication No. 59-100812 published on June 11, 1984. 
     In the system of this publication, a geomagnetic direction sensor is used for measuring the traveling direction of the vehicle. The geomagnetic direction sensor has a pair of windings which are intersected in perpendicular to each other and arranged on an annular core in a horizontal posture. The geomagnetic direction sensor outputs detected voltages (output values) corresponding to geomagnetic directional components interlinked to the respective windings. 
     When the vehicle makes one turn under the uniform geomagnetism, an output circle is described on a coordinate plane using coordinate positions each defined by the voltages detected at the respective windings. Accordingly, during a normal traveling of the vehicle, a direction from a center of the output circle toward a coordinate position defined by the detected voltages is derived as a traveling direction of the vehicle. 
     On the other hand, when the vehicle body is magnetized due to the disturbance of the environmental geomagnetism, the center of the output circle is caused to vary or move, thus resulting in an error in measurement of the traveling direction of the vehicle. To correct this error, one turn of the vehicle is required to derive a new center of a corresponding output circle. Specifically, the output values of the geomagnetic direction sensor are sampled during one turn of the vehicle and then averaged to derive a new or corrected center of the output circle. 
     In the foregoing background art, however, because a small number of the output values are sampled to be averaged even when the magnetic environment is not preferable, it is substantially impossible to derive a corrected center coordinate position of a corresponding output circle with high accuracy. Further, one turn of the vehicle has to be performed by a driver for the correction. 
     Still further, although it is desirable to make one turn of the vehicle as soon as possible after the vehicle body is magnetized, it is practically difficult to seek a proper place for making the vehicle turn. Accordingly, the direction measurement error remains during the traveling of the vehicle until finding the proper place for making the vehicle turn. 
     Further, when the foregoing correction of the center of the output circle is performed in the condition of the large disturbance of the environmental geomagnetism, such as, on a road under a high level road or among high buildings, the correction accuracy becomes less in comparison with the correction accuracy performed in the condition of smaller disturbance of the geomagnetism. Accordingly, the corrected coordinate position of the center of the output circle tends to be unreliable so that the correction itself becomes meaningless. 
     SUMMARY OF THE INVENTION 
     Therefore, it is an object of the present invention to provide a system for measuring a traveling direction of a vehicle that can eliminate the above-noted defects inherent in the background art. 
     It is another object of the present invention to provide a system for measuring a traveling direction of a vehicle than can provide an accurate measurement of a traveling direction of the vehicle by correcting a coordinate position of a center of an output circle with high accuracy, but without requiring a vehicle driver to make any turn of the vehicle, irrespective of the condition of the disturbance of the environmental geomagnetism around the vehicle. 
     To accomplish the above-mentioned and other objects, according to one aspect of the present invention, a system comprises a geomagnetic sensor for sequentially deriving two directional component data of the geomagnetism around a vehicle, the two directions being orthogonal to each other on a horizontal plane, first means for deriving a traveling direction of the vehicle based on a direction toward a coordinate position defined by the two directional component data from a first coordinate position of a center of a corresponding output circle. 
     The system further comprises second means for deciding whether a magnetization level on a vehicle body is shifted from a first level to a second level based on the two directional component data, third means for storing a first group of the two directional component data, fourth means for storing a second group of the two directional component data, fifth means for deriving a second coordinate position as a preliminarily corrected value of the first coordinate position based on the first group of the data, the second coordinate position derived for the first magnetization level in response to a negative decision of the second means and sixth means for deriving a third coordinate position as a preliminarily corrected value of the first coordinate position based on the second group of the data, the third coordinate position derived for the second magnetization level in response to a positive decision of the second means. 
     The system further comprises seventh means for deriving a first value indicative of accuracy of the second or third coordinate position, eighth means for deriving a second value which is variable depending on a variation in the magnetization level on the vehicle body, and ninth means for calculating a fourth coordinate position as a finally corrected value of the first coordinate position based on the second or third coordinate position and a latest value of the fourth coordinate position, the latest value derived in a prior calculation executed by the ninth means, by changing a rate of displacement of the fourth coordinate position between the latest value and the second or third coordinate position, the rate determining a point of the fourth coordinate position between the latest value and the second or third coordinate position. 
     According to another aspect of the present invention, a system comprises a geomagnetic sensor for sequentially deriving two directional component data of the geomagnetism around a vehicle, the two directions being orthogonal to each other on a horizontal plane, first means for deriving a first traveling direction of the vehicle based on a direction toward a coordinate position defined by the two directional component data from a first coordinate position of a center of a corresponding output circle, second means for deriving a variation in a traveling direction of the vehicle based on data output from a non-geomagnetic second sensor, the second sensor outputting the data without being affected by the geomagnetism around the vehicle, and third means for deriving a second traveling direction of the vehicle based on the variation derived by the second means. 
     The system further comprises fourth means for deciding whether a magnetization level on a vehicle body is shifted from a first level to a second level based on the two directional component data, fifth means for storing a first group of the two directional component data, sixth means for storing a second group of the two directional component data, seventh means for storing a third group of the two directional component data and corresponding data of the second traveling direction of the vehicle, eighth means for deriving a second coordinate position as a preliminarily corrected value of the first coordinate position based on the first group of the data, the second coordinate position derived for the first magnetization level in response to a negative decision of the fourth means, ninth means for deriving a third coordinate position as a preliminarily corrected value of the first coordinate position based on the data stored in the seventh means, the third coordinate position derived for the second magnetization level in response to a positive decision of the fourth means, tenth means for deriving a first value indicative of accuracy of the second coordinate position, eleventh means for deriving a second value indicative of accuracy of the third coordinate position, twelfth means for deriving a fourth coordinate position as a preliminarily corrected value of the first coordinate position based on the second group of the data, the fourth coordinate position derived for the second magnetization level when the second value is less than a first predetermined level, thirteenth means for deriving a third value indicative of accuracy of the fourth coordinate position, and fourteenth means for deriving a fourth value which is variable depending on a variation in the magnetization level on the vehicle body. 
     The system further comprises fifteenth means for calculating a fifth coordinate position as a finally corrected value of the first coordinate position based on the second or third or fourth coordinate position and a latest value of the fifth coordinate position, the latest value derived in a prior calculation executed by the fifteenth means, by changing a rate of displacement of the fifth coordinate position between the latest value and the second or third or fourth coordinate position, the rate determining a point of the fifth coordinate position between the latest value and the second or third or fourth coordinate position. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The present invention will be understood more fully from the detailed description given hereinbelow and from the accompanying drawings of the preferred embodiment of the invention, which are given by way of example only, and are not intended to be limitative of the present invention. 
     In the drawings: 
     FIG. 1 is a schematic block diagram showing a preferred embodiment of the present invention; 
     FIG. 2 is a schematic circuit wiring diagram of a geomagnetic direction sensor used in the preferred embodiment of FIG. 1; 
     FIG. 3 is a graph of magnetic flux density B versus magnetic field H showing characteristics of energization of the geomagnetic direction sensor of FIG. 2; 
     FIG. 4 is a graph of time versus magnetic flux showing a variation in magnetic flux on a permalloy core of the geomagnetic direction sensor under no application of the magnetic field; 
     FIG. 5 is a graph of time versus magnetic flux showing a variation in magnetic flux on a permalloy core of the geomagnetic direction sensor under application of the magnetic field; 
     FIG. 6 is a graph of time versus detected voltage showing characteristics of the detected voltages; 
     FIG. 7 is a schematic plan view of a vehicle incorporated with the geomagnetic direction sensor of FIG. 2 to explain a traveling direction of the vehicle; 
     FIG. 8 is an explanatory view of an output circle derived by outputs of the geomagnetic direction sensor; 
     FIG. 9 is an explanatory view showing the geomagnetic direction sensor, wherein a magnetic field other than the geomagnetic field is applied to the geomagnetic direction sensor in addition to the geomagnetic field; 
     FIG. 10 is an explanatory view showing a displacement of the output circle due to magnetization on the vehicle body; 
     FIG. 11 is a general flowchart to be executed by an arithmetic processing unit in the preferred embodiment of FIG. 1 for deriving a new center position of a corresponding output circle; 
     FIG. 12 is a flowchart showing an interrupt routine to be executed per a predetermined time by the arithmetic processing unit in the preferred embodiment of FIG. 1; 
     FIG. 13 is a flowchart showing an interrupt routine to be executed per a predetermined distance by the arithmetic processing unit in the preferred embodiment of FIG. 1; 
     FIG. 14 is a flowchart showing an interrupt routine of a normal mode process to be executed by the arithmetic processing unit per a predetermined timing when new data enough for calculating a preliminary center position are collected; 
     FIG. 15 is a flowchart showing a subroutine of a normal mode process and a magnetization mode process for deriving a preliminary center position and its calculation accuracy; 
     FIG. 16 is a flowchart showing an interrupt routine for deriving a new center position of a corresponding output circle; 
     FIG. 17 is a flowchart showing a subroutine of initialization for a magnetization mode process; 
     FIG. 18(a) and FIG. 18(b) are a flowchart showing an interrupt routine of a magnetization mode process to be executed by a predetermined timing when new data enough for calculating a preliminary center position are collected; 
     FIG. 19 is a flowchart showing a subroutine of the magnetization mode process for deriving a preliminary center position and its calculation accuracy; 
     FIG. 20(a) is a graph showing a V x  -V y  coordinate plane for explaining how to derive a new preliminary center position of a corresponding output circle in the normal mode process; 
     FIG. 20(b) is a graph showing a V x  -V y  coordinate plane for explaining how to derive a new preliminary center position of a corresponding output circle in the magnetization mode process; 
     FIG. 21 is an explanatory view for explaining a matrix buffer; 
     FIG. 22 is an explanatory view for explaining how to derive a new center position using an interior division rate; and 
     FIG. 23 is a graph showing a V x  -V y  coordinate plane for explaining concepts for deriving a new preliminary center position using the method of isosceles triangle. 
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     Referring now to the drawings, a preferred embodiment of a vehicular traveling direction measuring system according to the present invention will be described with reference to FIGS. 1 to 23. 
     FIG. 1 shows a block diagram of the vehicular traveling direction measuring system of the preferred embodiment. In FIG. 1, the geomagnetic direction sensor 1, preferably of a flux gate type, resolves the geomagnetism into two directional components, the two directions being perpendicular to each other on a horizontal plane, and outputs the two directional component data in the form of electrical signals (VX, VY) indicative of a coordinate position on a X-Y coordinate plane. A vehicle speed sensor 2 outputs a count signal per a predetermined unit distance ΔD travelled by the vehicle. The gyro sensor 3 outputs a signal indicative of an angular velocity ω G  in a yawing direction of the vehicle. The gyro sensor 3 is replaceable by, such as, a steering angle sensor and a speed differential sensor which derives a difference in speed between right and left side wheels, as long as its output is not affected by the geomagnetic environment. 
     The output signals from the geomagnetic direction sensor 1, the speed sensor 2 and the gyro sensor 3 are respectively fed to an arithmetic processing unit 5 through an interface circuit 4. The processing unit 5 executes several calculation routines based on the input data to derive a traveling direction of the vehicle, which will be described later in detail. 
     FIG. 2 shows a structure of the geomagnetic direction sensor 1 used in the preferred embodiment. The geomagnetic direction sensor 1 includes an annular magnetic core 6 made of a permalloy and a pair of windings 7X, 7Y wound on the magnetic core 6. The windings 7X, 7Y are arranged in perpendicular to each other. Further, a winding 8 is wound on the magnetic core 6. The winding 8 is energized by a power supply 9 so as to energize the magnetic core 6 immediately before saturation as shown in FIG. 3. 
     When the geomagnetic direction sensor 1 is placed in a non-magnetic field with the power supplied to the winding 8, a magnetic flux φ 1  and a magnetic flux φ 2  passing through positions S 1  and S 2  of the permalloy core 6, respectively, have the same magnitudes but run in opposite directions to each other (180° out of phase with each other), as shown in FIG. 4. Accordingly, as the magnetic flux interlinked to the winding 7X becomes zero, its detecting voltage expressed by V x  -N dφ/dt (N denotes the number of turns of the winding), also becomes zero, and similarly, a detecting voltage V y  of the winding 7Y also becomes zero. 
     Further, when the geomagnetic field He is applied in perpendicular to the winding 7X as shown in FIG. 2, a bias of the geomagnetic flux density expressed as Be=μHe (μ denotes a permeability of the permalloy core) is given to the magnetic flux in the magnetic core 6. Accordingly, the magnetic flux φ 1  and φ 2  becomes asymmetrical as shown in FIG. 5, so that the winding 7X derives the detecting voltage V x  of a waveform as shown in FIG. 6. 
     On the other hand, because the geomagnetic field He is parallel to the winding 7Y in FIG. 2, no geomagnetic field is interlinked with the winding 7Y so that no voltage V y  is generated in the winding 7Y. 
     The geomagnetic direction sensor 1 is mounted on the vehicle in a horizontal posture as shown in FIG. 7. For example, when the geomagnetic field He is applied to both of the windings 7X, 7Y, the detected voltages V x  and V y  (output values) are respectively derived at the windings 7X, 7Y according to the geomagnetic field, as shown in FIG. 8. 
     The detected voltages V x  V y  are expressed by the following equations (1) and (2) provided that a value K is a winding constant and a value B is a horizontal component of the geomagnetism He. 
     
         V.sub.x =KB cos φ                                      (1) 
    
     
         V.sub.y =KB sin θ                                    (2) 
    
     Accordingly, an angle θ of a traveling direction of the vehicle with respect to a width direction of the vehicle body is expressed by the following equation (3). 
     
         θ=tan.sup.-1 (V.sub.x /V.sub.y)                      (3) 
    
     As appreciated from the equations (1) and (2), when the vehicle makes one turn under the uniform geomagnetic field He, an output circle is described on the X-Y coordinate plane based on a coordinate position defined by the detected voltages V x  and V y , as shown in FIG. 8. The output circle is expressed by the following equation (4). 
     
         V.sub.x.sup.2 +V.sub.y.sup.2 =(KB).sup.2                   (4) 
    
     Accordingly, because the coordinate position defined by the detected voltages V x , V y  falls on the output circle, a direction from the center 0 of the output circle toward the coordinate position (an output point shown in FIG. 8) is derived by the arithmetic processing unit 5 as a traveling direction of the vehicle. 
     Assuming that the vehicle is magnetized, an additional magnetic field G is generated and interlinked with the windings 7X, 7Y along with the geomagnetism He, as shown in FIG. 9. Accordingly, the output circle is moved from a broken line position to a solid line position, as shown in FIG. 10. Consequently, the traveling direction of the vehicle derived by the processing unit 5 inevitably includes an error because a coordinate position derived after the magnetization on the vehicle body does not fall on the output circle derived before the vehicle magnetization. 
     In order to correct the above-noted measurement error, the processing unit 5 executes a below-described process for deriving a corrected coordinate position of the center of a corresponding output circle. 
     A general operation of the preferred embodiment will be described hereinbelow with reference to FIG. 11 which shows a general flowchart to be executed by the processing unit 5 for deriving a corrected center coordinate position of a corresponding output circle. 
     In a first step 100, it is determined whether the vehicle is magnetized. This determination is executed based on the input data in the form of the output values V x  and V y  fed from the geomagnetic sensor 1 and stored in an annular buffer A which will be described later. Specifically, when the disturbance of the environmental geomagnetism large enough to magnetize the vehicle occurs, the output values V x , V y  of the geomagnetic direction sensor 1 abruptly varies. Accordingly, by calculating the variation between the output values V x , V y  and comparing it with a predetermined threshold level, it is determined whether the vehicle is magnetized or whether the vehicle is shifted to a new magnetization level. In other words, the processing unit 5 determines whether a magnetization level including zero level is shifted from a current level to a new magnetization level. In order to derive an accurate determination whether the shift of the magnetization level occurs, the processing unit 5, for example, first calculates an average value of a predetermined number (for example, 20) of V x  stored in the annular buffer A and an average value of a predetermined number (for example, 20) of V y  stored in the annular buffer A. Then, the processing unit 5 calculates a distance between a coordinate position defined by newly input V x  and V y  and a coordinate position defined by the average values. When the calculated distance is above a predetermined level, the processing unit 5 determines that the vehicle is magnetized to a new magnetization level. 
     The details of the process for yielding the above-noted determination are disclosed in Japanese First Patent Publication No. 64-46612 published on February 21, 1989 and in W. German Patent Application No. P 40 00 345.0 filed on January 8, 1990, the contents of which are incorporated by reference for the sake of disclosure. 
     If the decision at the step 100 is NO, i.e. there is no possibility of the substantial shift of the magnetization level, the routine goes to a step 102 where a condition of a mode flag is checked. If the mode flag is not set at 1, i.e. being set at 0, then the routine goes to a step 104 where a normal mode process is executed. On the other hand, if the mode flag is set at 1, then the routine goes to a step 106 where a magnetization mode process is executed. 
     In the normal mode process at the step 104, a preliminary coordinate position of a center of a corresponding output circle as well as its calculation accuracy are derived, which will be described later in detail. At a step 108, it is checked whether the above-noted preliminary coordinate position and its calculation accuracy are derived at the step 104. If YES, then the routine goes to a step 110 where the mode flag is cleared to zero in the case of the normal mode process, because the mode flag is set at 0 at the step 102, the mode flag is remained 0 at the step 116 and further goes to a step 112 where a finally corrected value indicative of the center coordinate position is derived based on the above-noted preliminary center position and the calculation accuracy, which will be also described later in detail. 
     On the other hand, if the decision at the step 100 is YES, i.e. there is a possibility of the substantial shift of the magnetization level, the routine goes to a step 114 where an initialize routine is executed, which will be described later in detail, and then goes to a step 116 where the mode flag is set to 1. Subsequently, the routine returns to the step 100. 
     If the decision at the step 100 is NO, then the routine goes to the step 102. This time, because the mode flag has been set to 1 at the step 116, the routine goes to the step 106 where a preliminary coordinate position of a center of a corresponding output circle as well as its calculation accuracy are derived, which will be described later in detail. At the step 108, it is checked, similarly for the normal mode process, whether the above-noted preliminary center position and its calculation accuracy are derived at the step 106. If YES, then the routine goes to the step 110 where the mode flag is cleared to zero, and further goes to the step 112 where a finally corrected value indicative of the center position is derived based on the above-noted preliminary center position and the calculation accuracy, which will be also described later in detail. 
     On the other hand, if the decision at the step 100 is again YES, then the steps 114 and 116 are repeated until the decision at the step 100 becomes NO. This is likely to happen when the vehicle successively passes through a plurality of railroad crossings. 
     FIG. 12 shows a flowchart of an interrupt routine which is executed per a predetermined unit time ΔT, for example, every 100 msec. 
     At a step 200, the output values V x , V y  derived by the geomagnetic sensor 1 and an angular velocity ω G  in a yawing direction of the vehicle derived by the gyro sensor 3 are read out every AT. At a step 202, the output values V x , V y  read out at the step 200 are stored in the annular buffer A. The annular buffer A sequentially stores a predetermined number (for example, 20) of a pair of the output values BV x , V y . When a new pair of the output values is stored, the oldest pair of the output values is discarded. Accordingly, the annular buffer A always stores the predetermined number of the newest pairs of the output values V x , V y . The geomagnetic data stored in the annular buffer A are used at the step 100 in FIG. 11 for determining whether the magnetization level is shifted, as described before. 
     Before going to a step 204, an interrupt routine for counting up a speed counter is explained with reference to FIG. 13. This interrupt routine is executed per a predetermined distance ΔD travelled by the vehicle. Specifically, if the speed sensor 2 outputs 24 pulses per one rotation of a vehicle tire, then ΔD is, for example, 6 to 7 centimeters, which however is variable according to an outer diameter of the tire, and the interrupt routine is executed 24 times per one rotation of the tire. At a step 300, the speed counter is counted up by 1 per every execution of the interrupt routine and is used at the step 204 in FIG. 12 for determining whether the vehicle is stopped. Specifically, when a value of the speed counter is zero, then it is determined that the vehicle is stopped. 
     Referring back to the step 204 in FIG. 12, when the answer at the step 204 is YES, i.e. the vehicle is stopped, then subsequent processes are not executed in order to avoid collections of the data while the vehicle is stopped. This is because if the data are collected during the vehicle being stopped, weight of these data increases when calculating the corrected center position so that the correction of the center position can not be effectively performed. 
     On the other hand, if the answer at the step 204 is NO, i.e. the vehicle is running, the routine goes to a step 206 where the speed counter is reset to zero for a subsequent execution of this interrupt routine. At a step 208, a gyro direction θ G  is derived by the following equation (5). 
     
         θ.sub.G =θ.sub.GOLD +ω.sub.G ×ΔT (5) 
    
     wherein, θ GOLD  is θ G  which has been obtained during the prior execution of the interrupt routine based on the outputs from the gyro sensor and stored in an annular buffer B at a step 224. 
     Next, the routine goes to a step 210 where a reliability of the data V x , V y  is derived. Specifically, the output data V x , V y  includes not only the regular geomagnetism components but also irregular components generated due to the disturbance of the geomagnetism around the vehicle. Naturally, these irregular components become an obstacle against the precise measurement of the vehicle traveling direction. Accordingly, by estimating an amount of the disturbance, the data reliability can be derived such that as the amount of the disturbance is estimated less, the data reliability is assumed higher. 
     The estimation of the disturbance amount is performed for example, by using a combination of the geomagnetic direction sensor 1 and the gyro sensor 3. The geomagnetic direction sensor 1 is capable of deriving an absolute traveling direction of the vehicle while its output is affected by the magnetism environment. The gyro sensor 3 is capable of deriving a variation in the traveling direction of the vehicle by integration of the angular velocity and accumulation thereof without being affected by the magnetic environment. The followings are examples how to derive the estimation of the disturbance amount. 
     (1) A variation in a direction derived by outputs from the gyro sensor 3 and a variation in a direction derived by outputs from the geomagnetic direction sensor 1 are compared to estimate the amount of the disturbance larger as a difference between two variations is larger. This is based on the fact that though the output of the gyro sensor 3 includes a drift error, it is considered to be accurate for a short time, i.e. for a small variation. 
     (2) A first distance travelled by the vehicle for ΔT between time points t 1  and t 2  is derived by Δθ G  ×R, wherein Δθ G  is a variation in a vehicle traveling direction for ΔT derived by the outputs from the gyro sensor 3 and R is a radius of a corresponding output circle. Next, a second distance travelled by the vehicle for AT between the time points t 1  and t 2  is derived by the following equation (6). 
     
         {(V.sub.x1 -V.sub.x2).sup.2 +(V.sub.y1 -V.sub.y2).sup.2 }.sup.1/2(6) 
    
     wherein, (v x1 , V y1 ) is an output value at the time point t 1  derived by the geomagnetic direction sensor, and (V x2 , V y2 ) is an output value at the time point t 2  derived by the geomagnetic directional sensor. 
     The first and second distances are compared to estimate the amount of the disturbance larger as a difference between the first and second distances becomes larger. 
     Further details are disclosed in Japanese Patent Application No. 63-211967 contents of which are hereby incorporated by reference for the sake of disclosure. 
     Then, at a step 212, the data reliability is compared with a predetermined value α 1 . If an answer at the step 212 is YES, i.e. the data reliability is less than the predetermined value, the subsequent process is not executed. If the answer at the step 212 is NO, i.e. the data reliability is no less than the predetermined value, the routine goes to a step 214 where a matrix buffer B is renewed. At a step 216, it is checked whether the mode flag is set at 1. If an answer is NO, a matrix buffer A is also renewed at a step 218. If the answer is YES, the matrix buffer A is not renewed. 
     Before going to a step 220, explanation will be given hereinbelow for the matrix buffers A and B with reference to FIG. 21. 
     Each matrix buffer has a V x  -V y  coordinate plane which is partitioned into a plurality of sections similar to a lattice pattern, and each section has a counter. Assuming that an X coordinate value V x  derived by the geomagnetic directional sensor 1 falls within a range of x n+1  to x n+2  and a corresponding Y coordinate value V y  falls within a range of y n  to y n+1 , the counter for the hatched section is counted up according to the foregoing data reliability. Accordingly, a value of the counter becomes larger with a higher reliability of the data and with a larger number of times of the corresponding section to be matched with the input data. 
     The matrix buffer A is used for the normal mode process and the matrix buffer B is used for the magnetization mode process, which will be described later. 
     Referring back to FIG. 12, at the step 220, a traveling direction θ of the vehicle is calculated by the following equation (7). ##EQU1## wherein, C x  and C y  are X and Y coordinate values of a center of an output circle derived by a prior execution of a center position correction routine as shown in FIG. 16, which will be described later. 
     The traveling direction θ derived at the step 220 is used at a step 400 in FIG. 14 for determining a timing of starting the execution of the normal mode process. 
     Next, at a step 222, the data reliability derived at the step 210 is compared with another predetermined value α 2  which is set larger than the predetermined value α 1 . If a decision at the step 222 is NO, i.e. the data reliability is above the predetermined value α 2 , the output value V x , V y  and the corresponding θ G  derived at the step 208 are stored in the annular buffer B at a step 224. Similar to the annular buffer A, the annular buffer B stores a predetermined number of the sequentially derived data sets. The data stored in the annular buffer B is used in the magnetization mode process for deriving the preliminary center coordinate position and its calculation accuracy, which will be described later in detail. 
     Now, the normal mode process will be described hereinbelow with reference to FIG. 14. 
     At the first step 400, a current traveling direction θ derived at the step 220 in FIG. 12 is compared with an old traveling direction θ OLD1  stored at a step 430 in a prior execution of the normal mode process. If a difference between 25 the directions θ and θ OLD1  is larger than a given threshold level THD 1  which is, for example, THD 1  =100°, a preliminary coordinate position PC of the corresponding output circle is calculated at a step 410, and further, a calculation accuracy AC of the preliminary center position PC derived at the step 410 is derived at a step 420, based on the data stored in the matrix buffer A. The matrix buffer A stores all the data during traveling of the vehicle except for those data to be erased in the magnetization mode process, which will be described later. Then, at the step 430, the current traveling direction θ is stored as θ OLD1  for a subsequent execution of the normal mode process. 
     On the other hand, if the answer at the step 400 is NO, i.e. the difference between the directions θ and θ OLD1  is less than the given threshold level THD 1 , the subsequent steps as mentioned above are not executed in order to avoid the calculation of the preliminary center value without new data enough for the calculation thereof. 
     It is to be noted that when power is first applied to the arithmetic processing unit 5, the old traveling direction θ OLD1  is initialized to a traveling direction derived at the step 220 at the time of the initialization. 
     Now, the details of the steps 410 and 420 will be described hereinbelow with reference to FIG. 15. 
     At a first step 500, V xo  and V yo  are respectively initialized to C x  and C y . C x  and C y  are virtual X and Y coordinate values of a center of an output circle derived by the data stored in the matrix buffer A based on a premise that a circle is described by the data stored in the matrix buffer A. 
     For example, C x  and C y  are derived by the following equation. ##EQU2## wherein, i denotes the number of X and Y coordinate values, respectively, derived by the geomagnetic sensor 1. 
     At the first step 500, a virtual radius R o  of the output circle having its center value of V xo  and V yo , is also initialized to 300 mG (milligauss) which is an average output value of the geomagnetism. 
     Through subsequent steps 510 to 560, the preliminary center coordinate position is calculated using the method of least squares. 
     Specifically, a distance Ri between the virtual center coordinate position (V xo , V yo ) and a coordinate position (V x1 , V yi ) derived by the output of the geomagnetic sensor 1 is calculated by the following equation (8). 
     
         Ri=[{V.sub.xi -V.sub.xo }.sup.2 +{V.sub.yi -V.sub.yo }.sup.2 ].sup.1/2(8) 
    
     Accordingly, a value J according to the method of least squares is derived by the following equation (9) at the step 510. ##EQU3## 
     According to the method of least squares, by deriving values of V xo , V yo  and R o  which make the value J the least, the corrected values of C x  and C y  can be obtained. 
     At the subsequent step 520, partial derivatives of J with respect to V xo , V yo  and R o  are respectively derived. 
     In order to derive the values of V xo , V yo  and R o  which render the value J the least, the following condition should be satisfied. ##EQU4## 
     By deriving values of V xo , V yo  and R o  which satisfy the equations (10), (11) and (12), a coordinate position defined by the derived V xo  and V yo  becomes the corrected preliminary center value PC. 
     However, it is impossible to derive the values of V xo , V yo  and R o  which satisfy the equations (10), (11) and (12). In view of this, the present preferred embodiment utilizes Newton-Raphson&#39;s method at the step 530. At the step 530, correction values h, m and f are derived by the following equations (13), (14), (15) and (16). ##EQU5## wherein, giz=δgi/δz. 
     At the subsequent step 540, the correction values are compared with predetermined reference values R 1 , R 2  and R 3 , respectively. If a decision at the step 540 is YES, i.e. all the correction values are no more than the corresponding reference values, then the X and Y coordinate values V xo  and V yo  derived at the step 500 are set as the corrected preliminary center position PC at the step 560. 
     On the other hand, if any one of the three correction values is decided above the corresponding reference value at the step 540, the correction value h is added to v xo , the correction value l is added to v yo , and the correction value m is added to R o  at the step 550. In this case, the steps 520 to 540 are repeated until each of the correction values becomes no more than the corresponding reference value. 
     The details of deriving the preliminary center position is disclosed in Japanese First Utility Model Publication No. 1-144814 published on Oct. 4, 1989, the contents of which are incorporated hereby by reference for the sake of disclosure. 
     Subsequently, at a step 570, the accuracy of calculation AC of the preliminary center coordinate position derived at the step 560 is calculated. This calculation of the accuracy is derived by using a function f(J) with respect to the value J. A value of the function f(J) gets less as the correction values h, m, l get larger, and is given in a range of 0 to 1. 
     Now, the correction routine for deriving a new or finally corrected center position referred to at the step 112 in FIG. 11 will be described in detail with reference to FIG. 16. 
     At a first step 500, an interior division rate DR is derived by the following equation (17). 
     
         DR=AC/IVC                                                  (17) 
    
     wherein, Ac denotes the calculation accuracy derived at the step 570 in FIG. 15, and IVC denotes a coefficient which is variable based on a variation in the magnetization level on the vehicle body, which will be described later in detail. 
     The interior division rate DR is a coefficient for determining to what point between the preliminarily center position PC and a latest center position OC derived by a previous execution of the center position correction routine, a new center position NC should be approached. In other words, the interior division rate DR determines a point of the new center position between the preliminarily center position PC and the latest center position OC. 
     It is to be noted that the calculation accuracy is given, as mentioned above, within the range of 0 to 1, and the coefficient IVC is given no less than a value 1, which will be described later. Accordingly, the interior division rate DR is given in a range 0 to 1. 
     Subsequently, at a step 610, the new center position NC is derived by the following equation (18). 
     
         NC=OC+DR×|PC-OC|                   (18) 
    
     wherein, OC denotes the latest cent position derived by a previous execution of the center position correction routine, DR denotes the interior division rate, and PC denotes the preliminary center position derived at the step 560 in FIG. 15. 
     As understood from the equation (18), the new center position NC approaches closer to the preliminary center position PC from the latest center position OC with a larger value of the interior division rate DR, as clearly shown in FIG. 22. 
     Subsequently, the coefficient IvC is increased when conditions of steps 620 and 630 are satisfied. Specifically, if the calculation accuracy AC derived at the step 570 in FIG. 15 is decided to be larger than a predetermined threshold level TH 2  at the step 620, and an absolute difference between the preliminary center position PC and the new center position NC is decided to be less than a predetermined threshold level THL 1 , a new coefficient NIVC is derived by the following equation (19). 
     
         NIVC=OIVC+K.sub.2                                          (19) 
    
     wherein, K 2  is a set constant, and OIVC denotes latest coefficient IVC derived at the step 640 in a prior execution of the center position correction routine and used at the foregoing step 600 in the current execution of this correction routine, or a latest coefficient IVC used at the step 600 both in the prior and current execution of the center position correction routine in which the increment of the coefficient IVC at the step 640 was not executed in the prior execution the center position correction routine. 
     Accordingly, when reliable preliminary center positions are continuously obtained through sequential executions of the normal mode process, the coefficient IVC gets considerably large through sequential executions of the center position correction routine. Therefore, even if the vehicle travels across the bad geomagnetic environment, i.e. under the geomagnetic disturbance not large enough to magnetize the vehicle, after having travelled long in the good geomagnetic environment, so that the preliminary center position calculated deviates from the latest center position derived through a prior execution of the center position correction routine, the new center position does not move substantially because the coefficient IVC at the step 600 is considerably large to make the interior division rate at the step 610 considerably small. Accordingly, a reliable and stable measurement of the vehicle traveling direction is attained. 
     It is to be appreciated that when the power is first applied to the arithmetic processing unit 5, the coefficient IVC is initialized to a given proper value. 
     Now, a center position correction system through the magnetization mode process will be described hereinbelow with reference to FIGS. 16 to 19. 
     FIG. 17 shows the initialize routine for the magnetization mode process as mentioned at step 114 in FIG. 11. 
     As described above, in the magnetization mode process, the calculation of the preliminary center position is executed based on only those data derived after the positive decision at the step 100 that the magnetization level is shifted. Accordingly, in this initialize routine, the data stored in the matrix buffer B and the annular buffer B are erased through respective steps 700 and 710. Subsequently, a traveling direction θ G  derived at the step 208 in FIG. 12 at the time of this initialization is stored as θ OLD2  at a step 720 and as θ OLD3  at a step 730. Accordingly, θ OLD2  equals θ OLD3  at the time of the initialization. 
     After the initialization is executed, the mode flag is set to 1 at the step 116 in FIG. 11 and then the routine goes to the magnetization mode process if the answer at the step 100 is NO, i.e. the possibility of occurrence of a further magnetization level shift is determined negative. 
     FIGS. 18(a); and 18(b) show the magnetization mode process referred to at the step 106 in FIG. 11. 
     At a first step 800, a current traveling direction θ G  derived at the step 208 in FIG. 12 is read out. Then, at a step 802, a difference between the current traveling direction θ G  and the traveling direction θ OLD2  stored at the step 720 in FIG. 17 is compared with a predetermined threshold level THD 2  which is set, for example, to 10°. If an answer at the step 802 is YES, i.e. the difference is larger than the threshold level THD 2 . then the routine goes to a subsequent process. On the other hand, if the answer at the step 802 is NO, then the routine returns to the step 800 to repeat the steps 800 and 802 in order to avoid the subsequent calculations without new data not enough for the calculations thereof, similar to the step 400 in FIG. 14. 
     It is to be noted that, at the step 802, a variation in the traveling direction is derived by the output data from the gyro sensor, while, in the normal process mode, the variation in the traveling direction is derived by the output data from the geomagnetic direction sensor 1. This is because, in the magnetization process mode, the center coordinate position is actually moved from the center coordinate position derived by the prior execution of the center position correction routine before the current magnetization level shift occurs and the traveling direction θ derived at the step 220 in the magnetization mode process is calculated using the center coordinate position derived before the current magnetization level shift occurs. Accordingly, the variation in the traveling direction can not be measured precisely based on the output data from the geomagnetic direction sensor 1. 
     It is also to be noted that although the traveling direction obtained by the output data of the gyro sensor includes a substantial error on a long-term basis due to a drift, it is considered to be accurate on a short-term basis in comparison with that derived by the output data from the geomagnetic direction sensor in the magnetization mode process. 
     Subsequently, at a step 804, a preliminary center position is calculated based on the data stored in the annular buffer B using a so-called method of isosceles triangle which will be described hereinbelow with reference to FIG. 19. 
     At a first step 900, the current traveling direction θ G  read out at the step 800 in FIG. 18(a) is stored as θ G1 . Subsequently, at a step 910, a traveling direction θ G  which has a maximum angular difference from the current traveling direction θ G1  among the traveling directions θ G  stored in the annular buffer B, is read out as θ G2 . At a next step 920, a directional difference Δθ G  between the traveling directions θ G1  and θ G2  is derived. 
     The preliminary center position can be derived in the following way. 
     As shown in FIG. 23, a coordinate position A defined by the stored values V x1  and V y1  corresponding to θ G1  and a coordinate position B defined by the stored values V x2  and V y2  corresponding to θ G2  are set on a VX-VY coordinate plane. Then, the preliminary center position C (C xc , C yc ) is derived such that the preliminary center position C is distanced by R1 from the respective coordinate positions A and B and an angle between two lines connecting the preliminary center position C and the respective coordinate positions A and B is Δθ G  (∠ACB=Δθ G ). 
     This is equivalent to deriving a preliminary center position by the following equation which is executed at a subsequent step 930. ##EQU6## 
     Subsequent to the step 930 where the preliminary center position is calculated as mentioned above, its calculation accuracy AC is derived at a step 940 using a function f(Δθ G ) with respect to the directional variation Δθ G . As seen from a graph of directional variation Δθ G  versus calculation accuracy AC shown at the step 940, a value of the function f(Δθ G ) gets larger as the directional variation θ G  gets larger, and is given in a range of 0 to 1. 
     The details of the so-called method of isosceles triangle described above is disclosed in Japanese First Utility Model Publication No. 1-117712 published on August 9, 1989, which is incorporated by reference for the sake of disclosure. 
     Referring back to FIG. 18(a), subsequent to a step where the calculation accuracy AC is derived as mentioned above, the routine goes to a step 808 where the current traveling direction θ G  is stored as θ OLD2  for a subsequent execution of the magnetization mode process if an answer at a step 812 is negative, which will be described later. 
     Subsequently, at a step 810, the derived calculation accuracy AC is compared with a predetermined threshold level TH 1 . If an answer at the step 810 is NO, i.e. the calculation accuracy AC is no less than the threshold level TH 1 , the routine goes to a step 818 bypassing subsequent steps 812 to 816. On the other hand, if the answer at the step 810 is YES, i.e. the calculation accuracy is less than the threshold level TH 1 , the routine goes to the step 812. This is because, since the calculation accuracy is less in the foregoing method of isosceles triangle than in the foregoing method of least squares, the routine skips the steps 812 to 816 only when the derived calculation accuracy AC is no less than the threshold level TH 1 . 
     At the step 812, a directional difference between a new traveling direction θ G  and the stored direction θ OLD3  is compared with a predetermined threshold level THD 3  which set, for example, to 90° in this embodiment. As appreciated, the threshold level THD 3  is considerably larger than the threshold level THD 2  which is set at 10°. If an answer at the step 812 is NO, i.e. the directional variation is less than the threshold level THD 3 , then the routine returns to the step 800 to repeat the foregoing process. On the other hand, if the answer at the step 812 is YES, i.e. the directional difference is no less than the threshold level THD 3 , then the routine goes to the steps 814 and 816 where the preliminary center position PC and its calculation accuracy AC are derived based on the data stored in the matrix buffer B using the foregoing method of least squares as in the same manner as described with reference to FIG. 15. 
     As appreciated, in the magnetization mode process, the preliminary center position is derived by the method of isosceles triangle per every small directional variation until the calculation accuracy AC gets no less than the threshold level TH 1 . If the calculation accuracy AC gets no less than the threshold level TH 1 , it is considered that the reliable preliminary center position PC is derived and the routine directly goes to the step 818. On the other hand, if it is continued that the calculation accuracy is less than the threshold level TH 1  and if an accumulated directional variation commencing from the initialization at the step 730 in FIG. 17 becomes no less than the threshold level THD 3 , the preliminary center position PC and its calculation accuracy AC are derived based on the data stored in the matrix buffer B using the method of least squares. 
     Referring now to the step 818 in FIG. 18(b), a coefficient VC indicative of a variation in the magnetization level shift on the vehicle body is derived by the following equation (21). 
     
         VC=K.sub.1 ×|PC-OC|                (21) 
    
     wherein, K 1  a set constant, PC denotes the preliminary center position, and OC is the latest center position derived through a prior execution of the center position correction routine as shown in FIG. 16. The coefficient VC is given in a range of 0 to 1 and is in proportion to a difference between the preliminary center position PC and the latest center position OC. Accordingly, the coefficient VC gets larger as the variation in the magnetization level gets larger. 
     Subsequently, at a step 820, the coefficient IVC, as mentioned before with reference to FIG. 16, is renewed. Specifically, a new coefficient NIvC is derived by the following equation (22). 
     
         NIVC=OIVC×(1-VC)                                     (22) 
    
     wherein, OIVC is a latest coefficient IVC derived at the step 640 in a prior execution of the center position correction routine of FIG. 16, or a latest coefficient IVC used at the step 600 in the prior execution of the center position correction routine in which the step 640 was not executed in the prior execution thereof. The coefficient IVC is given no less than a value 1. 
     As appreciated from the equation (22), the new coefficient NIVC is set in a less value as the variation in the magnetization level gets larger, i.e. as the coefficient VC gets larger. 
     Subsequently, at a step 822, the data stored in the matrix buffer A are erased according to a rate determined by the coefficient VC. For example, each counter of the matrix buffer A is multiplied by (1-VC). Accordingly, if the coefficient VC is at its maximum value 1, then all the data stored in the matrix buffer A are erased at the step 822. 
     Subsequently, at a step 824, all the data stored in the matrix buffer B is added to the matrix buffer A. 
     The execution of the steps 822 and 824 is necessary for a subsequent execution of the normal mode process under the magnetization level substantially the same as the magnetization level under which the current magnetization mode process as mentioned above has been executed. 
     Accordingly, as shown in FIG. 20(a), in the normal mode process, a preliminary center position PC is derived based on all the data vo stored in the matrix buffer A. On the other hand, as shown in FIG. 20(b), in the magnetization mode process, an output circle C 1  which is derived by the data V 1  collected before a shift of the magnetization level moves to an output circle C 2  which is derived by the data V 2  collected after the shift of the magnetization level. Accordingly, the preliminary center position in the magnetization mode process should be derived based on only those data collected after the shift of the magnetization level. In FIG. 20(b), O 3  denotes a virtual center position derived based on mixture of the data V 1  and V 2 . 
     Now, the center position correction routine of FIG. 16 after the execution of the magnetization mode process will be described hereinbelow. It is to be appreciated that the center position correction routine itself is the same as that after the execution of the normal mode process. 
     At the 600, the interior division rate DR is derived by the following equation (17). 
     
         DR=AC/IVC                                                  (17) 
    
     wherein, AC denotes the calculation accuracy derived at the step 806 or 816 in FIG. 18(a), and IVC denotes the new coefficient NIVC derived at the step 820 in FIG. 18(a). 
     As described above, the interior division rate DR is a coefficient for determining to what point between the preliminary center position derived at the steps 804 or 814 and a latest center position derived by a previous execution of the center position correction routine, a new center position should be approached. In other words, the interior division rate DR determines a point of the new center position between the preliminary center position and the latest center position. 
     The calculation accuracy is given, as mentioned above, within the range of 0 to 1, and the coefficient IVC is given no less than a value 1. At the step 820, it is set that the minimum value of the new coefficient NIVC is 1. Accordingly, the interior division rate DR is given in a range 0 to 1. 
     Subsequently, at the step 610, the new center position NC is derived by the following equation (18). 
     
         NC=OC+DR×|PC-OC|                   (18) 
    
     wherein, OC denotes the latest center position derived by a previous execution of the center position correction routine, DR denotes the interior division rate, and PC denotes the preliminary center position derived at the step 804 or 814 in FIG. 18 (a); 
     As understood from the equation (18), the new center position approaches closer to the preliminary center position from the latest center position with a larger value of the interior division rate, as clearly shown in FIG. 22. Since the coefficient IVC is set small at the step 820 in FIG. 18(b), the interior division rate DR is set large at the step 600 so that the displacement of the center position gets larger toward the preliminary center position PC in the center position correction routine after the execution of the magnetization mode process. 
     Subsequently, the coefficient IvC is increased when conditions of the steps 620 and 630 are satisfied. Specifically, if the calculation accuracy AC derived at the step 806 or 816 in FIG. 18 (a) is decided to be larger than the predetermined threshold level TH 2  at the step 620, and the absolute difference between the preliminary center position PC and the new center position NC is decided to be less than a predetermined threshold level THL 1 , a new coefficient NIVC is derived by the following equation (19). 
     
         NIVC=OIVC+K.sub.2                                          (19) 
    
     wherein, K 2  is a set constant and OIVC den latest coefficient IVC derived at the step 820 as the new coefficient NIVC in FIG. 18(b) and used at the foregoing step 600 in the current execution of this correction routine. 
     Generally, since the calculation accuracy is not so high in the magnetization mode process in comparison with that derived in the normal mode process when the vehicle is traveling across the good geomagnetic environment and since the threshold level TH 2  is set quite high, the step 640 is unlikely to be executed through the execution of the magnetization mode process. 
     It is to be appreciated that, according to the foregoing preferred embodiment of the present invention, when reliable preliminary center positions are continuously obtained through sequential executions of the normal mode process, the coefficient IvC gets considerably large through sequential executions of the center position correction routine. Therefore, even if the vehicle travels across the bad geomagnetic environment, i.e. under the geomagnetic disturbance not large enough to magnetize the vehicle, after having travelled long in the good geomagnetic environment, so that the preliminary center position calculated deviates from the latest center position derived through a prior execution of the center position correction routine, the new center position does not move substantially because the coefficient IVC at the step 600 is considerably large to make the interior division rate at the step 610 considerably small. Accordingly, a reliable and stable measurement of the vehicle traveling direction is attained. On the other hand, immediately after the magnetization level on the vehicle body has shifted to a new magnetization level, i.e. in the magnetization mode process, the center position should be approached to a new center position, i.e. toward the preliminary center position without delay. This is accomplished by reducing the coefficient IVC drastically at the step 820 in FIG. 18(b) so as to derive a large interior division rate DR at the step 600 in FIG. 16. 
     Further, as appreciated from the foregoing description of the preferred embodiment, in the center position correction system through the normal mode process, not so substantial variation in the magnetization level, i.e. the variation less than the predetermined level, which is otherwise accumulated to the substantial variation no less than the predetermined level, can be gradually compensated each time the new center position is derived, while, in the center position correction system through the magnetization mode process, the substantial variation in the magnetization level no less than the predetermined level is drastically compensated. 
     It is to be understood that this invention is not to be limited to the embodiments described above, and that various changes and modifications may be made without departing from the spirit and scope of the invention as defined in the appended claims.