การแปลงพื้นหลักฐานระหว่าง WGS 84 และ Indian 1975 และ การคำนวณความสูง Geoid โดยใช้ EGM96 และ EGM2008

เอกสารที่ต้องอ่าน
การสำรวจโดยการรังวัดดาวเทียม GPS ในประเทศไทย
การแปลงพื้นหลักฐานระหว่าง WGS 84 และ Indian 1975
ค่าความสูงออร์โธเมตริกจากการสํารวจด้วยดาวเทียมระบบ GPS ในประเทศไทย
สถานีโครงข่ายหลักจีพีเอสกรมโยธาธิการและผังเมือง
หมุด GPS กรมโยธาธิการและผังเมือง
งานรังวัดและทำแผนที่ของกรมที่ดิน
การรังวัดหมุดหลักฐานแผนที่โดยระบบดาวเทียม
แผนที่ทั้งหมดมาตราส่วน 1:50 000 (กรมแผนที่ทหาร)
แผนที่ทั้งหมดมาตราส่วน 1:250 000 (กรมแผนที่ทหาร)

การแปลงพื้นหลักฐานระหว่าง WGS 84 และ Indian 1975 ( Datum

Transformation between WGS 84 and Indian 1975)

เรียบเรียงโดย พ.อ.เอื้อมเกียรติ เจริญสม จาก“ รายงานผลการศึกษาวิจัย เรื่องค่าตัวแปรในการเปลี่ยนพื้นหลักฐานของ ผท.ทหาร : WGS 84 กับ อินเดียน 1975 ” โดยกองยีออเดซี่และยีออฟิสิกส์ กรมแผนที่ทหาร

1. กล่าวนำ

ตั้งแต่ปีพ.ศ.2534 กรมแผนที่ทหารได้เปลี่ยนวิธีการสำรวจขยายโครงข่ายหมุดหลักฐานทางราบ

แห่งชาติจากวิธีการสำรวจด้วยกล้อง กล่าวคือ การสำรวจสามเหลี่ยมชั้นที่ 1 และการวงรอบชั้นที่ 1 มาเป็น

การสำรวจด้วยดาวเทียม GPS ประเทศไทยยังใช้พื้นหลักฐานทางราบเป็นพื้นหลักฐาน Indian 1975 ซึ่งมีจุดศูนย์กำเนิดอยู่ที่เขาสะแกกรัง จังหวัดอุทัยธานี และมีรูปทรงรี คือ Everest ในขณะที่การสำรวจด้วย

ดาวเทียม GPS เกี่ยวข้องกับพื้นหลักฐาน WGS 84 ซึ่งมีจุดศูนย์กำเนิดหรือจุดอ้างอิงและรูปทรงรีที่

แตกต่างกัน ทำให้ค่าพิกัดของหมุดหลักฐานที่รังวัดได้ของหมุดเดียวกันแตกต่างกัน ดังนั้นจึงต้องมีการ

แปลงค่าพิกัดจากพื้นหลักฐาน WGS 84 จากการสำรวจด้วยดาวเทียม GPS มาเป็นค่าพิกัดบนพื้นหลักฐาน Indian 1975 ที่เป็นพื้นหลักฐานที่ใช้กันแพร่หลายในหน่วยงานที่เกี่ยวข้องกับการสำรวจและการทำแผนที่ ทั้งภาครัฐและเอกชน ในระยะเริ่มแรกของการสำรวจด้วยดาวเทียม GPS เพื่อการขยายโครงข่ายของหมุดหลักฐาน ทางราบจึงใช้หมุดหลักฐานที่เขาสะแกกรังเพียงหมุดเดียวเป็นหมุดแรกออก โดยใช้การแปลงพื้นหลักฐานจาก Indian 1975 มาเป็น WGS 84 ซึ่งได้รับความช่วยเหลือจาก Defense Mapping Agency (DMA) แห่งสหรัฐอเมริกา จากการใช้ข้อมูลของการรังวัดสถานี Doppler ของระบบดาวเทียม Transit ได้ค่าความสัมพันธ์ดังนี้

X75 = X84 - 206 เมตร

Y75 = Y84 - 837 เมตร

Z75 = Z84 - 295 เมตร

โดยที่ค่า X75 , Y75 , Z75 คือค่าพิกัดฉากบนพื้นหลักฐาน Indian 1975

X84 , Y84 , Z84 คือค่าพิกัดฉากบนพื้นหลักฐาน WGS 84

อย่างไรก็ดีหลังจากกรมแผนที่ทหารได้รับความช่วยเหลือจากองค์กรทางด้านยีออเดซี่จาก

ต่างประเทศหลายหน่วยงานมากขึ้น โดยเฉพาะในการสำรวจโครงข่ายหมุดหลักฐานชั้น Zero Order หรือClass AA ตามมาตราฐานของ FGCC (Federal Geodetic Control Committee ) แห่งสหรัฐอเมริกาที่มีการรังวัดเป็นแบบสัมบูรณ์ที่มีความละเอียดถูกต้องสูงในระดับที่สามารถตรวจสอบการเคลื่อนตัวของเปลือกโลกจนได้หมุดหลักฐานจำนวน 7 หมุด ซึ่งมีทั้งความละเอียดถูกต้องสูงกว่า ปริมาณมากกว่า และกระจายตัวดีกว่าในแง่ของการเป็นหมุดควบคุมในงานขยายโครงข่ายจากเดิมที่ใช้จุดศูนย์กำเนิดที่สะแกกรังเพียงแห่งเดียวเป็นหมุดควบคุม กรมแผนที่ทหารจึงได้เชื่อมโยงหมุดหลักฐานเหล่านี้เป็นหมุดควบคุมในการขยายโครงข่ายในระดับความละเอียดถูกต้องลดหลั่นลงไปอีกสองระดับอีกทั้งในปี พ.ศ. 2545 ได้รวบรวมข้อมูลการรังวัด GPS ทั้งหมดนับตั้งแต่การสำรวจด้วย GPSนำมาปรับแก้พร้อมกันจนได้โครงข่ายที่มีความถูกต้องน่าเชื่อถือและเป็นเอกภาพผลการดำเนินการดังกล่าวทำให้ค่าพิกัดของหมุดหลักฐานเปลี่ยนแปลงไปจากเดิมที่ใช้หมุดควบคุมเพียงแห่งเดียวแล้วขยายโครงข่ายและปรับแก้ออกเป็นส่วนๆ จึงได้ศึกษาวิจัยเพื่อหาความสัมพันธ์ในการการแปลงพื้นหลักฐานทั้งสองขึ้นใหม่เพื่อให้มีความละเอียดถูกต้องน่าเชื่อถือมากขึ้น

2. สมมติฐานในการแปลงพื้นหลักฐานฯ

ในหลายกรณีการแปลงพื้นหลักฐานจะใช้แบบจำลองทางคณิตศาสตร์ แบบ 7 ตัวแปร (การ

แปลงค่าพิกัดฉากแบบ 3 มิติแบบคงรูป)ที่ประกอบด้วย การเลื่อนตามแกน 3 ตัวแปร การหมุนแกน 3 ตัว

แปร และความต่างมาตราส่วนอีก 1 ตัวแปร แต่ในกรณีจะยึดถือรูปแบบของการแปลงที่ DMA ให้การ

ช่วยเหลือกล่าวคือเป็นแบบการเลื่อนตามแกน 3 ตัวแปร นอกจากนี้โดยปกติค่าการหมุนและความต่าง

มาตราส่วนมักมีค่าน้อยมากรวมทั้งขนาดพื้นที่ของประเทศไทยไม่ใหญ่มากนักเมื่อเทียบกับขนาดของโลก

จึงใช้เลือกแบบจำลองทางคณิตศาสตร์ดังกล่าวในการหาค่าตัวแปรในการแปลงพื้นหลักฐานดังนี้

X75 = X84 + Δ X เมตร

Y75 = Y84 + Δ Y เมตร

Z75 = Z84 + Δ Z เมตร

โดยที่ค่า Δ X , Δ Y และ Δ Z เป็นการเลื่อนตามแกน X , Y และZ ตามลำดับ

3. ขั้นตอนในการหาค่าตัวแปรในการแปลงพื้นหลักฐาน

3.1 ข้อมูลที่ใช้ในการหาค่าตัวแปร

การจะหาค่าตัวแปรได้หลังจากที่กำหนดแบบจำลองทางคณิตศาสตร์ที่เหมาะสมตาม

สมมติฐานแล้วขั้นตอนต่อไป คือ การรวบรวมข้อมูลที่ใช้ในการหาค่าตัวแปร ข้อมูลดังกล่าวนี้ คือ ค่าพิกัดของหมุดหลักฐานที่ทราบค่าบนพื้นหลักฐานทั้งสองที่มีค่าน่าเชื่อถือสูง และกระจายตัวอย่างดีโดยแปลงจากค่าพิกัดแบบ Geodetic Coordinates ( ϕ , λ,h ) มาเป็นค่าพิกัดฉาก 3 มิติ (X , Y , Z) ในการนี้สามารถ หาได้จำนวน 16 หมุด ตามตารางที่ 1 ปัญหาอยู่ที่ว่าค่าพิกัดของหมุดเหล่านี้บนพื้นหลักฐาน Indian 1975 ซึ่งสำรวจด้วยกล้องไม่สามารถหาค่าความสูงเหนือทรงรี(Ellipsiodal Height; h) ได้มีเพียงแต่ค่า Orthometric Height (H) เมื่อไม่มีความสูงเหนือทรงรี จึงไม่สามารถไปแปลงเป็นค่าพิกัดฉาก X , Y , Z ได้ อย่างถูกต้องดังนั้นจึงหาค่าความสูงเหนือทรงรีจากผลรวมของ และ Geoid Undulation (N) ที่ได้จาก Earth Gravity Model (EGM 96) จะได้ h = H + N ในบางหมุดรวมกับการใช้เทคนิคในการปรับแก้โครงข่ายฯ จึงได้ h บนหมุดหลักฐานทั้ง 16 หมุด

การคำนวณความสูง Geoid โดยใช้ EGM96 และ EGM2008

• ขอพักเรื่องโปรแกรมมิ่งสักตอนเพราะว่า เขียนเกี่ยวกับเรื่องโปรแกรมมิ่งต้องใช้พลังความคิดมาก นอกจาก library “GDAL” ที่ผมกำลังเขียนถึง ซึ่งมีแง่มุมให้เขียนเกี่ยวกับการใช้งานได้เป็นร้อยๆตอนเลยละครับ ถ้ามีแรงกายและใจขนาดนั้น ถ้ายังจำกันได้คือ GeographicLib ของ Charles Karney ที่เคยอ้างอิงถึงไปแล้ว มีเรื่องอีกเรื่องที่น่าสนใจคือ การคำนวณ Geoid Height ซึ่งดูโค๊ดของ Charles Karney โมดูลนี้ต้องเป็นงงเรื่องอัลกอริทึ่ม มากไปกว่านั้นมีรุ่นน้องส่งโค๊ดมาให้คำนวณเรื่องเดียวกันแต่เป็นโค๊ดของภาษา Fortran ยิ่งงงเข้าไปอีก ทั้งที่เคยร่ำเรียนตอนอยู่มหาวิทยาลัย แต่ก็ลืม syntax ไปหมด
• EGM (Earth Gravitational Model) ที่คนที่รังวัด GPS คงต้องทราบกันดี รุ่นก่อนหน้านี้คือ EGM96 ใช้งานแพร่หลายที่สุด โปรแกรมคำนวณ GPS ทั้งหลายเช่น Leica Geo Office (LGO) หรือ Trimble Geo Office (TGO) ผมยังใช้ LGO5 ยังเป็นรุ่นเก่าอยู่ไม่ทราบว่าในรุ่นใหม่ ทั้ง TGO และ LGO ปรับมาใช้ EGM2008 แล้วหรือยัง
• ที่อยากเขียนเรื่องคำนวณ Geoid Height เพราะในการใช้ GPS ในปัจจุบันการอ้างอิงความสูงจะเป็นความสูงบน Ellipsoid ก็คือ WGS84 ถ้าบอกว่าความสูงบน WGS84 เท่ากับ 35 เมตร จะเท่าไหร่เมื่อเมื่อเทียบกับ รทก. (ระดับน้ำทะเลปานกลาง หรือ MSL – Mean sea level) ที่เราคุ้นเคยกันดี การคำนวณ Geoid height จึงเป็นเรื่องที่คำนวณการทอนความสูงจาก WGS84 มาบน MSL (หรือเรียกอีกอย่างว่า Orthometric Height) แต่คนใช้ระดับชาวบ้าน บางครั้งไม่ต้องสนใจเพราะในเครื่อง GPS มือถือรุ่นใหม่ทั้งหลายได้คำนวณความสูงจากทรงรี WGS84 เป็น ให้เป็น MSL ให้เป็นที่เรียบร้อยแล้ว (คือภายในเครื่องได้บรรจุ EGM96 เข้าไปเรียบร้อย)

รูปทรง Geoid

• รูปทรงรีที่ใช้แทนสัณฐานโลก จะมีรูปทรงที่แน่นอน สามารถใช้สูตรทางคณิตศาสตร์คำนวณได้ แต่ Geoid เป็นจะเป็นทรงที่ไม่แน่นอน บุบๆบู้บี้ ดูรูปข้างล่างประกอบ

รูปทรง Geoid (ภาพจาก http://en.wikipedia.org/wiki/Geoid)

ความสัมพันธ์ระหว่าง Geoid และ ทรงรี

• เนื่องจากค่าความสูงที่เราต้องการคือต้องอยู่บน Geoid (เทียบเท่าระดับน้ำทะเลปานกลาง) แต่ความสูงที่ได้จากรังวัด GPS (ตัวอย่างเช่นการรังวัดแบบ Fast Static และ Static) จะอยู่บนทรงรี ซึ่งแต่ละสถานที่ความต่างเนื่องรูปทรงที่ไม่แน่นอนของ Geoid จึงทำให้แต่ละสถานทีค่าความสูงต่างจากทรงรีและบน Geoid จะไม่เท่ากัน

ความสัมพันธ์ Geoid และทรงรี (ภาพจาก http://principles.ou.edu/earth_figure_gravity/geoid/index.html)

• จากรูปด้านบน(ผมชอบรูปนี้มากแสดงได้ชัดเจน) ค่าความสูง H (รทก.) = h (ความสูงบนทรงรี ได้จาก GPS) – N (ค่าความสูง Geoid) จากสูตรความสูงของจุดที่เราต้องการจริงๆแล้วอยู่บนภูมิประเทศ (Topography หรือ Terrain) ถ้าเรารังวัด GPS จะได้ค่า h ออกมา ถ้าสามารถคำนวณหา N ได้จากค่าพิกัด geographic (lat/long) เราก็สามารถหาค่า H ที่เราต้องการได้ จากรูปด้านบนจะสังเกตเห็นว่าแนวคำนวณหา H จะเป็นแนวตั้งฉากกับ Geoid ซึ่งจะเป็นแนวแรงไปตามแรงโน้มถ่วงเสมอ แต่แนว h จะเป็นแนวตั้งฉากกับทรงรี ลองดูรูปอีกรูปหนึ่งด้านล่าง

ความสัมพันธ์ระหว่าง Geoid และทรงรี (ภาพจาก http://en.wikipedia.org/wiki/Geoid)

ความเป็นมาของ EGM (Earth Gravitational Model)

• EGM96 เป็นผลงานของความร่วมมือของ National Imagery and Mapping Agency (NIMA), NASA Goddard Space Flight Center (GSFC) และ Ohio State University โครงการนี้ได้รวบรวมข้อมูลความโน้มถ่วงจากหลายๆแหล่งของโลก ด้วยวิธีที่แตกต่างกันบ้าง เช่นในมหาสมุทรใช้ดาวเทียม GEOSAT และ ERS-1 ข้อมูลเหล่านี้จะนำมาหาค่าที่เรียกว่า coefficients ของ EGM96
• แรงโน้มถ่วงของโลกมีลักษณะเป็นฟังก์ชั่น Harmonic  ซึ่งฟังก์ชั่นฮาร์โมนิค จะคำนวณได้ก็ต่อเมื่อทราบค่า coefficients ที่กล่าวไว้ข้างต้นนั่นเอง โมเดล EGM96 จะประกอบไปด้วยกริดแต่ละช่องขนาด 15 ลิปดา x 15 ลิปดา และแต่ละช่องจะเก็บค่า coefficients ที่แตกต่างกันไปตามค่าพิกัด การคำนวณหาความสูง geoid (Geoid Undulation) จะรับ Input จากผู้ใช้คือค่าพิกัด แล้วนำค่าพิกัดไปดึงเอาค่า coefficient แล้วนำค่า coefficient ไปคำนวณหาค่าความสูง Geoid (N – เรียกอีกอย่างว่า Geoid separation)

EGM2008

• จัดทำโดย U.S. National Geospatial-Intelligence Agency (NGA) ชื่อเดิมก็คือ NIMA นั่นเอง โดยได้ปรับปรุงจาก EGM96 โดยเพิ่มข้ิอมูลจากดาวเทียม Gravity Recovery and Climate Experiment(GRACE) เป็น ดาวเทียมของ NASA ที่วัดสนามความโน้มถ่วงของโลกปล่อยสู่วงโคจรในปี 2002 โดยที่โมเดลนี้ได้เผยแผ่ในปี 2008 จึงเรียกว่า EGM2008 ความละเอียดของ EGM2008 ประกอบด้วยกริดที่แต่ละช่องมีขนาด 1′ x 1′  หรือประมาณบรรจุค่า coefficients ประมาณ 4 ล้านค่า

EGM96 vs. EGM2008

• EGM2008 ควรจะมี accuracy สูงกว่า EGM96 แต่บางรายงานก็ยังก้ำกึ่ง ตัวอย่างรายงานและการทดสอบที่แสดงว่า EGM2008 นั้นดีกว่า EGM96 ได้แก่ Wuhan University หรือที่ โดยศาสตราจารย์  Charles Merry, University of Cape Town ท้ายข้อสรุประบุว่ายากที่จะบอกว่า EGM2008 ละเอียดกว่า EGM96 มากเท่าไหร่ แต่ที่ทวีปอเมริกาเหนือ ทวีปยุโรป ค่า Accuracy ของ EGM2008 ดีถึง ±10 ซม.

การคำนวณความสูง Geoid บน EGM96

• หาโปรแกรมที่เป็นโปรแกรม standalone ไม่ได้ส่วนใหญ่แล้วจะแฝงอยู่ในโปรแกรมด้านการแปลงพิกัด เช่น Trimble Coordinate Calculator หรือโปรแกรมด้านการคำนวณรังวัด GPS เช่น TGO และ LGO นอกเหนือจากนั้นเป็นโปรแกรม online ซะมากกว่า ที่แรกก็คือที่ NGA โดยตรงไปที่ http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html ดูตัวอย่างด้านล่าง