Recently, many papers have been published on the earthquake-related ionospheric anomalies. In particular, event studies and statistical studies on Ionospheric electrons using Ionosonde and GPS TEC (total electron content) indicate the pre-seismic anomalies. In Japan region, the previous statistical study during 1998-2010 shows that there is an increase anomaly 1-5 days before the earthquake with M>=6 D<=40 km. In this study, we performed statistical analysis on ionosonde data observed at Kokubunji (35.71N, 139.49E), Japan over 60 years of 1958 ~ 2017 in order to investigate the earthquake-related ionospheric anomalies. The ionosonde is operated by the National Institute of Information and Communication Technology, Japan. The studied parameters are NmF2, which indicates the maximum electron density in the F2 layer, hmF2, which shows the NmF2 altitude, and foEs which is the critical frequency of the Es layer. We select earthquakes occurred within 1000 km and 350 km radius from the ionosonde station, and M >= 6 and D <= 40 km. There are 352 and 148 earthquakes. Ionospheric anomalies are also induced by geomagnetic storm. Therefore, to remove the geomagnetic storm effects, we adopt the Dst index. We divided the values into some classes by onset time and intensity. For each class, we performed statistical processes to calculate duration of ionospheric anomalies triggered by a geomagnetic storm. We found these anomalies last about 1-3 days, and we remove the data based on results. To investigate the correlation between ionospheric anomalies and earthquakes, we perform the Superposed Epoch Analysis (SEA). We define the anomaly as 15 days backward median value ±1.5 IQR and an anomalous day has the anomaly with more than 10 hours. Then, we make up the binary time series data on the earthquake-related anomaly. We extract dataset 30 days before and after the earthquake day for the all selected earthquakes and perform SEA. To evaluate the statistical significance, we compute the random test. SEA results show that there is a significant positive anomalies 1-10 days before the earthquakes in NmF2 data. For hmF2 data, we can find there are no clear correlation. As for foEs, we need further investigation. For NmF2, we checked Molchan’s Error Diagram (MED). From the results of MED, NmF2 may be precursor for earthquakes. The detailed will be given in the presentation.