Tag: คำนวณ

การเล็งสกัดย้อน (Resection) ด้วยการวัดมุมภายใน ระยะทางและและมุมแบริ่งด้วยวิธีการคำนวณแบบ Least Squares (ตอนที่ 2)

 ตั้งสมการ Observation Equation

ขอทบทวน ค่า aik, bik  เรียกว่า  direction coefficients และ  cik, dik เรียกว่า distance coefficients

ในกรณีวัดมุมเล็งสกัดย้อนจากสมการด้านบนและเอาแทนที่ในสมการด้านล่าง

เขียนให้ดูง่ายดังนี้  zi  คือค่าอะซิมัทเริ่มต้น

เราจะมาคำนวณหาค่า aik, bik กันก่อน มาคำนวณที่จุด P ค่าพิกัด N = 193939.897 E = 110879.464 ไปสถานีหลักฐานจุดที่ 1  ที่มีค่าพิกัด N = 192315.290 E = 120383.500 คำนวณระยะทางได้ 9641.890 เมตร คำนวณหาอะซิมัทได้ 99°42’1.1″ ดังนั้น aik = aP1 = -sin( 99°42’1.1″) / (9641.890 * 100) * 3600 * 180/3.141592654 = -0.2109 second (มุมแปลงเป็นหน่วยฟิลิปดา ระยะทางแปลงหน่วยเป็น ซม.)

bik = bP1 = cos(99°42’1.1″) / (9641.890 * 100) * 3600 * 180/3.141592654 = -0.0360 second

เราจะฟอร์มสมการในรูปแบบ v + Bx = f โดยที่

v คือเมตริกเวคเตอร์ของ residual

B คือเมตริกของค่า coefficient ของมุม ระยะทาง

f คือเมตริกความต่างของค่าที่คำนวณและค่าที่รังวัด

ผมคำนวณหาค่าaik, bik ทุกๆการรังวัดมุมมาดังนี้

มาทบทวนดูสมการระยะทางดังนี้้

คำนวณหา cik = cP1 = cos(99°42’1.1″) = -0.1685

dik= dP1 = sin(99°42’1.1″) = 0.9857

คำนวณทุกการรังวัดระยะทางมาได้ดังตาราง

สมการสุดท้ายคือการวัดเล็งสกัด (Bearing intersection)

จากสมการ v + Bx = f มาพิจารณาฝั่งซ้าย v + Bx ก่อน ผมฟอร์มเป็นเมตริกดังนี้

มาดูเมตริก f เป็นเมตริกแสดงความต่างระหว่างค่าที่คำนวณกับค่ารังวัด f = Computed – Observation มาดูตามตาราง น่าจะเข้าใจได้ง่าย ค่า diff ในตารางก็คือค่า f นั่นเอง โปรดระวังหน่วยมุมจะเป็นฟิลิปดา (second) หน่วยระยะทางใช้เป็นซม.

ยกตัวอย่างค่าของเมตริก f เริ่มจากการวัดเล็งสกัดย้อน จะใช้มุมอะซิมัทที่ได้จากการคำนวณมาเป็นตัวเริ่มต้น ค่าจากจุด P ไปสถานีหลักฐานที่ 1 จะได้ค่าความต่างเท่ากับ 0.0 ต่อไปจากจุด P ไปสถานีหลักฐานที่ 2 คำนวณได้ 119.5116959 จากมุมที่วัดมา 99° 42′ 1.1″ + 19° 48′ 41″ = 119.5116957 ได้ค่าความต่าง = (119.5116959 – 119.5116957) * 3600 =  0.0008 second

ที่ง่ายที่สุดคือวัดระยะทางจากจุด P ไปสถานีหลักฐานจุดที่ 1 ได้ 9641.795 เมตร ส่วนการคำนวณจากค่าพิกัด P เริ่มต้นมายังค่าพิกัดสถานีหลักฐาน 1 ได้ค่าคำนวณ = 9641.890 เมตร ความต่าง = (9641.890 – 9641.795) * 100 = 9.5 ซม.

สุดท้ายสามารถนำค่ามาเขียนเป็นเมตริก ดังนี้

v + Bx = f

จะเห็นว่าสมการนี้ติดค่า residual ของเมตริก v และเมตริก x ไม่สามารถคำนวณต่อไป แต่หัวใจของ least squares ตามชื่อเลยครับคือผลรวมค่ายกกำลังสองของ  residual ที่ได้ค่าน้อยที่สุด

ถ้าค่า weight หรือน้ำหนักของการรังวัดไม่เท่ากันจะต้องคูณน้ำหนักเข้าไปด้วย

จากสมการ v + Bx = f แทนค่า v = f – Bx ในสมการ

ค่า Φ จะมีค่าน้อยที่สุด ดังนั้นจะหาอนุพันธ์ (ดิฟ) โดยที่ให้ x มีค่าเท่ากับศูนย์

ที่นี้ก็จำง่าย Nx = t โดยที่ N = BTWB และ f = BTW

สุดท้ายสามารถหาค่าเมตริก x = N-1t เมื่อได้ค่า x แล้วก็สามารถย้อนไปหาเมตริก residual (v) ได้ ตามสมการ v + Bx = f

ผลลัพธ์การคำนวณรอบที่ 1

ผมใช้ฟังก์ชั่นของ excel หาเมตริกได้ดังนี้


ได้เมตริก x คือค่าปรับแก้หน่วยเป็นซม. เอาพิกัดจุด P และมุมอะซิมัท เริ่มต้นมาปรับได้ดังนี้ N = 193939.897 – 0.090829= 193939.806 ค่า E = 110879.464 + 0.04695 =  110879.511 และค่ามุมอะซิมัท = 99°42’1.1″ – 2.1156″ = 99° 41′ 58.99″ (99.6997191)

ผลลัพธ์การคำนวณรอบที่ 2

เอาค่าพิกัดของจุด P ที่ได้จากรอบที่ 1 มาเป็นตัวเริ่มต้น พร้อมทั้งมุมอะซิมัทด้วย

ตั้งสมการเมตริก v + Bx = f ได้ดังนี้

แก้สมการ Nx = t ได้ดังนี้

จะได้เมตริก x ค่าใหม่ เอาพิกัดจุด P และค่าอะซิมัทมาปรับได้ดังนี้ N =193939.806  – 0.0016 = 193939.790  และ E = 110879.511 + 0.0085 = 110879.519  ค่าอะซิมัท = 99° 41′ 58.99″ + 0.137″ = 99°41′ 59.13″ (99.6997572)

ผลลัพธ์การคำนวณรอบที่ 3

ต่อไปฟอร์มรูปเมตริก v + Bx = f สังเกตว่าเมตริก B ค่าเปลี่ยนไปเล็กน้อยมาก

คำนวณหาเมตริก x จากสมการ Nx = t => x = N-1t

จะได้เมตริก x ค่าใหม่ เอาพิกัดจุด P และค่าอะซิมัทมาปรับได้ดังนี้ N = 193939.790- 0.0003 = 193939.787  และ E = 110879.519 + 0.0001 = 110879.521 ค่าอะซิมัท = 99° 41′ 58.99″ + 0.024″ = 99°41′ 59.15″ (99.6997639)

มาถึงตอนนี้จะเห็นว่าค่าปรับแก้ ในเมตริก x  น้อยมากอยู่ในระดับเศษส่วนของมิลลิเมตร ΔN = -0.285 cm ΔE = 0.152 cm ดังนั้นแสดงว่าค่าที่คำนวณมานั้น convergence แล้ว ดังนั้นผมสรุปว่าผลลัพธ์ดังนี้

ค่าพิกัดจุด P (Free Station)

N = 193939.787 E = 110879.521 ค่าอะซิมัทจากจุด P ไปสถานีหลักฐานที 1 = 99°41′ 59.15″

 ตรวจสอบผลลัพธ์การคำนวณด้วย Microsurvey StarNet

เพื่อให้มั่นใจว่าผลลัพธ์ที่ได้จะถูกต้อง ผมใช้ Microsurvey StarNet  มาเป็นตัวช่วย ข้อมูลสถานีรังวัดไม่เกิน ดังนั้นผมยังใช้เวอร์ชั่นทดลองใช้ได้อยู่ เมื่อเปิดโปรแกรมมาผมใช้เมนู Options -> Project ตั้งค่าดังนี้ เปลี่ยนหน่วยเป็นเมตรให้เรียบร้อยก่อนที่ Adjustment > Units > Linear > Meters

จากนั้นปรับ  standard error  ของกล้องวัดมุมและระยะทางให้สอดคล้อง ผมปรับระยะทาง Distance constant 4mm ปรับ Distance PPM = 20 ความหมายตัวนี้คือ 4mm ± 20mm/1,000,000 x L (m) ถ้าวัดระยะทาง 1000 เมตร error จะอยู่ประมาณ 24 mm (กล้อง Total Station สมัยปัจจุบันเรื่องระยะทางทำได้ดีกว่านี้มาก)

จากนั้นป้อนข้อมูลไปดังนี้ หมายเหตุว่า C = Control Point ไม่มีเครื่องหมาย ! !  ตามหลังแสดงว่าเป็นค่าเริ่มต้นหรือประมาณการ A = Angle, D = Distance และ B = Bearing ถ้าสนใจก็ไปดาวน์โหลดโปรแกรมของ Microsurvey StarNet  มาทดลองได้

แล้วคัดลอกข้อมูลด้านล่างแล้วไปวางลงในหน้า  input เพื่อลองคำนวณดูผลลัพธ์กันได้

# Resection 2D combined resection, distance and bearing intersection.

# Approximate coordinates for the unknown free station
#
C P 193939.897 110879.464

# Coordinates for the known stations
C 1 192315.290 120383.500 ! !
C 2 189545.730 118642.430 ! !
C 3 188084.770 112278.210 ! !
C 4 190640.940 109654.540 ! !
C 5 190044.860 108065.980 ! !
C 6 194455.370 110632.930 ! !
A P-1-2 19-48-41
A P-1-4 100-40-19
A P-1-5 116-8-36
A P-1-6 234-44-22

#Distance measurements
D P-1 9641.795
D P-3 6019.802
D P-5 4804.793

#Bearing measurements
B 1-P 279-41-59.5
B 3-P 346-33-52
B 5-P 35-50-34

จากนั้นใช้เมนู Run > Adjustment

สรุปสถิติการคำนวณดังนี้ จะเห็นว่ามีจำนวนข้อมูลรังวัดมา 10 (วัดมุมเล็งสกัดย้อน 4 มุม ระยะทาง 3 ระยะ วัดมุมเล็งสกัด 3 มุม จำนวนสิ่งที่ไม่รู้ค่า 2 ค่า คือค่าพิกัด  x,y ของจุด P

Adjustment Statistical Summary
==============================

Iterations = 2

Number of Stations = 7

Number of Observations = 10
Number of Unknowns = 2
Number of Redundant Obs = 8

 

Adjusted Coordinates (Meters)

Station N E Description
P 193939.788830 110879.521213
1 192315.290000 120383.500000
2 189545.730000 118642.430000
3 188084.770000 112278.210000
4 190640.940000 109654.540000
5 190044.860000 108065.980000
6 194455.370000 110632.930000

จะได้ค่าพิกัดจุด  P ดังนี้ N = 193939.789 E = 110879.521 ต่างจากที่ผมคำนวณมาเล็กน้อย อย่างแรกคือผมไม่ได้ใช้น้ำหนัก  weight  แต่ของ MicroSurvey Starnet บังคับใช้ โปรแกรมไม่ได้แสดงค่าสมการให้ดู จึงตรวจสอบไม่ได้ว่าเมตริกของ W เป็นอย่างไร อย่างที่สองคือการคำนวณวนลูปสองรอบเท่านั้น แสดงว่าอัลกอริทึ่มอาจจะใช้อนุกรมเทเลอร์ลำดับที่ 2 ด้วย ซึ่งผมใช้ลำดับเดียวจึงต้องวนลูปมากกว่า

เป็นอย่างไรบ้างครับ เรื่องคำนวณ  least squares adjustment ตอนสมัยเรียนมหาวิทยาลัยถือว่าเป็นวิชาที่ยากสำหรับนักศีกษา แต่ถ้ามีโอกาสได้ลองคำนวณดูไปทีละขั้นตอน สำคัญคือต้องลงมือเองด้วย จึงจะเข้าใจ บางทีรายละเอียดเช่นอนุกรมเทเลอร์ การหาค่าอนุพันธ์ อาจจะยาก ผมไม่ได้จำหรอกครับ เพียงแต่ติดก็กลับมาเปิดตำรา หลงๆลืมๆไปบ้างตามวัยของผม แต่ปัจจุบันความได้เปรียบคือตำราหาอ่านได้ง่ายในอินเทอร์เน็ต โปรแกรมช่วยคำนวณแบบ  Microsurvey Starnet หรือโปรแกรมค่ายอื่นๆก็มีมากมาย ที่จะมาช่วยทำให้การคำนวณง่าย แต่ถ้ามีพื้นฐานบ้างก็จะได้เปรียบ  ติดตามกันต่อไปครับ

การเล็งสกัดย้อน (Resection) ด้วยการวัดมุมภายใน ระยะทางและและมุมแบริ่งด้วยวิธีการคำนวณแบบ Least Squares (ตอนที่ 1)

จากที่เขียนโปรแกรมสำหรับเครื่องคิดเลขคำนวณเล็งสกัดย้อนสำหรับเครื่องคิดเลข Casio fx-9860G II SD ทำให้นึกถึงวิธีการคำนวณแบบ least squares ที่เป็นพื้นฐานเคยร่ำเรียนมา โดยเฉพาะการรังวัดในปัจจุบันที่การรังวัดระยะทางด้วยกล้องประมวลผลรวมทำได้ง่าย เมื่อรวมกับการรังวัดมุม จะทำให้มีค่าเกินหรือ redundant มาคำนวณในวิธีแบบ least square ได้ การรังวัดแบบเล็งสกัดย้อนบางตำราเรียกว่า free station

การวัดเล็งสกัดย้อนในยุคแรกจะวัดมุมจากหมุดสถานีหลักฐานกันเป็นหลัก และหมุดสถานีหลักฐานต้องมีค่าพิกัดและมีจำนวนอย่างน้อยสามหมุด ดังนั้นการวัดมุมอย่างน้อยสองมุมจะเพียงพอต่อการมาคำนวณ ในบทความนี้ผมจะพาไปทัวร์การคำนวณวิธี least squares ด้วยการใช้การรังวัดผสมประกอบด้วยการวัดมุมภายใน วัดระยะทาง และการวัดมุมแบริ่งหรือการวัดภาคทิศ จากนั้นจะยกตัวอย่างพร้อมทั้งวิธีการคำนวณ ปิดท้ายทดสอบผลลัพธ์การคำนวณด้วยโปรแกรม Microsurvey StarNet

การวัดแบบ 2D

การวัดแบบนี้จะวัดมุมราบและระยะราบก็พอนำมาคำนวณหาค่าพิกัดของ free station ได้ ตัวอย่างที่ผมจะแสดงการคำนวณในลำดับต่อไปจะเป็นการวัดแบบ 2D เพื่อลดความซับซ้อน

การวัดแบบ 3D

การวัดระยะทางจะวัดแบบ Slope distance วัดมุมดิ่ง ถ้าทราบค่าระดับของหมุดสถานีหลักฐานและวัดความสูงของเป้าที่ตั้งบนสถานีหมุดหลักฐาน จากนั้นวัดความสูงกล้อง เมื่อนำมาคำนวณร่วมกับมุมราบแล้วจะได้ค่าพิกัดของ free station รวมทั้งค่าระดับด้วย

แนะนำตำรา

เมื่อผมกลับไปดูวิธีการคำนวณด้วย least square แต่กลับพบกับความนะจังงัง ว่าความรู้ที่รับการประสิทธิประสาทมาบัดนี้ได้คืนท่านอาจารย์ที่มหาวิทยาลัยไปจนหมดแล้วสิ้นเชิง นั่งงงกับการตั้งสมการ Observation Equation อยู่นานพอสมควร เมื่อไม่เป็นผลก็ต้องกลับไปหาตำราค่อยๆพลิกความทรงจำกลับมาใหม่

ตำราที่จะแนะนำให้อ่านเป็นของ Mr.Rod Deakin (Rodney Edwin Deakin) ถ้าดูตามโปรไฟล์ก่อนจะเกษียณเป็นอาจารย์ให้มหาวิทยาลัยสอนเรื่อง Engineering Survey มาก่อน ตำราหรือบทความที่แต่งนั้นมีความหลากหลายมากในเว็บไซต์ส่วนตัวนี้ อ่านง่าย เพลินจนลืมไปว่าเรื่องที่อ่านนั้นยาก

ตามไปดูหน้า least squares ลองเลื่อนไปด้านล่างดูเรื่อง “Notes on Least Squares” คือตำราที่จะมาแนะนำกันมีทั้งหมด 10 บท เนื่องจากผมพอมีพื้นฐานมาบ้างเล็กน้อย จึงไม่ได้อ่านเรียงหน้าตั้งแต่บทที่ 1 แต่อาศัยข้ามไปอ่านบทที่ 7 ก่อน พอเจอเมตริกผมงงก็ข้ามไปอ่าน “Appendix A” เฉพาะเรื่องเมตริก บวก ลบ คูณ อินเวอร์ส แล้วค่อยกลับมาอ่านบทอื่นๆอีกที ดังนั้นคนที่ห่างเรื่องนี้นานๆแนะนำให้ไปดูเรื่องเมตริกอันดับแรก

บทที่ 7 อนุกรมเทเลอร์

สำหรับการคำนวณเล็งสกัดย้อนด้วยวิธี least square ปฐมบทจะอยู่ที่บทที่ 7 เริ่มต้นจากอนุกรมเทเลอร์ เหตุที่เราต้องใช้เพราะว่าสมการที่จะนำมาคิด  resection นั้นไม่ใช่สมการแบบเชิงเส้นหรือ linear equation ดังนั้นจะต้องมีการถอดสมการเชิงเส้นออกมาและอยู่ในรูปอนุกรม เพื่อให้สามารถนำมาคำนวณได้ จะมีตัวอย่างสมการ Observation Equation 2 ตัวอย่างคือ

  • วัดมุม (Measure direction)
  • วัดระยะทาง (Measure distance)

บทที่ 9 คำนวณเล็งสกัดย้อน (Least Squares Resection)

บทนี้จะเริ่มตั้งแต่ตั้งสมการ  Observation Equation (ไม่เชิงเส้น) ของการวัดมุม จาก free station ไปยังหมุดสถานีหลักฐาน จากนั้นจะถอดสมการไม่เชิงเส้นด้วยอนุกรมเทเลอร์ แล้วก็มีตัวอย่างแสดงวิธีการคำนวณ

บทที่ 10 คำนวณเล็งสกัด (Least Squares Bearing Intersection)

ผู้อ่านอาจจะสังเกตสองคำคือเล็งสกัดกับเล็งสกัดย้อน

  • เล็งสกัดย้อน (Resection) คือไปตั้งกล้องที่  free station แล้ววัดมุมภายในไปหาสถานีหมุดหลักฐาน

  • เล็งสกัด (Bearing Intersection) จะเป็นการวัดมุมแบริ่ง ตัวอย่างได้แก่การวัดดาวเหนือในสมัยก่อน คือเราจะมีหมุดคู่เป็น base line ต้องการทราบมุมแบริ่งของ base line นี้ก็ตั้งกล้องที่หมุดตัวแรกแล้วอาศัยวัดดาวเหนือกับหมุดอีกตัวบน base line ดังนั้นถ้าวัดแบริ่งสัก 2 base line ไปตัดกันก็จะได้ค่าพิกัด

การถอดสมการเล็งสกัดย้อนด้วยอนุกรมเทเลอร์

ต่อไปมาลองดูสมการที่เราใช้คำนวณหามุมและระยะทาง โดยที่ทราบค่าพิกัด ซึ่งเป็นสมการพื้นฐานของงานสำรวจ ในที่นี้เป็นสมการไม่เชิงเส้น เราจะมาถอดสมการออกมาเป็นเชิงเส้นด้วยอนุกรมเทเลอร์ แผนผังด้านล่างเป็นงานเล็งสกัดย้อน ตั้งกล้องที่จุด Pi วัดมุมไปหาสถานีหลักฐานที่ทราบค่าพิกัดคือจุด  P1, P2, P3, … Pk

การวัดมุมในทางปฏิบัติจะส่องไปหาสถานีแรก P1 ตั้งมุมราบเป็น 0 แล้วกวาดไปสถานีที่ 2 วัดมุมได αi1 กวาดไปหาสถานีหลักฐาน P2 วัดมุมได้ α12 ทำเช่นนี้เรื่อยๆ จนถึงสถานีหลักฐาน Pk จะได้มุม αik แต่ในงานสำรวจด้านตรีโกณมิติเราจะไม่คำนวณโดยใช้มุมภายในแต่จะใช้มุมอะซิมัท ดังนั้นมุมภายในเหล่านี้ต้องเปลี่ยนเป็นเทอมของมุมอะซิมัททั้งหมด ซึ่งไม่ได้ยุ่งยากอะไร ถ้าทราบมุมอะซิมัทที่ไปหาสถานีแรก P1 ก็จะทราบอะซิมัทสถานีที่เหลือทั้งหมด แต่ในความเป็นจริงเราไม่ทราบ เพราะว่าค่าพิกัด P1 ยังไม่ทราบในตอนนี้

สำหรับสมการคำนวณหาอะซิมัทและระยะทางระหว่างจุด Pi และ Pk เมื่อทราบค่าพิกัดคำนวณได้ดังสมการ

เป็นสมการไม่เชิงเส้นครับ ทำการถอดสมการด้วอนุกรมเทเลอร์สมการมุมอะซิมัทมาได้ดังนี้ หมายเหตุผมจะพาไปแบบรวบรัด รายละเอียดอ่านได้ที่ตำราที่แนะนำไปข้างต้น

โดยที่ E = E0 + ΔE และ N = N0 + ΔN ส่วน E0 และ N0 คือค่าเริ่มต้นเป็นค่าประมาณการ ส่วนค่า ΔE, ΔN เป็นค่าปรับแก้ (correction) ซึ่งสองตัวนี้จะคำนวณได้จากสมการ least squares ที่กำลังจะว่ากันต่อไป และถ้าดูสมการจะเห็นว่าติดค่าอนุพันธ์ (ดิฟ) ซึ่งค่าดิฟหาได้ดังนี้

จากสมการด้านบน aik, bik เรามีค่าประมาณการเริ่มต้น Ei, Ni ดังนั้นสามารถหาค่าพวกนี้ได้ มาดูสมการระยะทางค่าดิฟ หาได้ดังนี้

จากสมการด้านบน cik, dik เรามีค่าประมาณการเริ่มต้น Ei, Ni ดังนั้นสามารถหาค่าพวกนี้ได้เช่นเดียวกัน

ตั้งสมการ Observation Equation สำหรับการวัดมุมภายใน (เล็งสกัดย้อน)

สมการ Observation Equation เป็นสมการที่เป็นหัวใจของการคำนวณ least squares ถ้าตั้งสมการนี้ได้ก็สามารถคำนวณต่อได้ แต่งานรังวัดถ้าวัดเกินมาจะมี residual เป็นส่วนเกิน สมการ Observation Equation ของการวัดมุมอะซิมัทแสดงได้ดังนี้

vik คือ residual คงไม่ลืมกันนะว่า มุมอะซิมัท Øik จากสถานีตั้งกล้องไปยังแต่ละสถานีหมุดหลักฐานเท่ากับ อะซิมัทไปหมุดสถานีหลักฐานตัวแรก Zi บวกด้วยมุมราบที่กวาดไป αik เขียนในรูปอนุกรมเทเลอร์ดังนี้

ในกรณีงานรังวัดเล็งสกัดย้อน ΔNk และ ΔEk ทราบค่าพิกัดดังนั้นเทอมตัวนี้จะเป็นศูนย์ กลายเป็นสมการนี้

ตั้งสมการ Observation Equation สำหรับการวัดระยะทาง

คล้ายๆกับสมการสำหรับการวัดทิศทาง แสดงได้ดังนี้

เขียนในรูปอนุกรมเทเลอร์

ตั้งสมการ Observation Equation สำหรับการวัดแบริ่ง (เล็งสกัด)

ถือว่าเป็นอาหารเสริมก็แล้วกัน สมัยนี้ยุค GNSS  คงหายากแล้วสำหรับการวัดภาคทิศอะซิมัทด้วยการวัดดาว

สมการเดียวกับที่เราใช้เล็งสกัดย้อน แต่ในที่นี้เราทราบค่าพิกัด Nk, Ek  ทำให้ ΔNk และ ΔEk เทอมตัวนี้จะเป็นศูนย์ และสมการจะคงเหลือดังนี้

ตัวอย่างงานรังวัดเล็งสกัดย้อนแบบรวมวัดมุม ระยะทาง และรังวัดภาคทิศ

มาถึงตอนนี้ผู้อ่านบางท่านเห็นสมการแล้วอาจจะงงงวยเหมือนผมในตอนแรกๆ มาลองแก้สมการดูจากตัวอย่างโจทย์ ที่นี้จะทำให้เข้าใจมากขึ้น ดูรูปด้านล่าง ตั้งกล้องที่จุด P ทำการรังวัดโดยที่วัดมุมไปยังสถานีหลักฐาน โจทย์ตัวนี้เป็นโจทย์ที่ผมเจอในอินเทอร์เน็ต เมื่อค้นหาการคำนวณเล็งสกัดย้อนแบบวิธีคำนวณ least square ดูได้ที่ลิ๊งค์นี้ โจทย์ที่จะคำนวณเปิดไปดูได้หน้าท้ายๆ

กำหนดค่าพิกัดหมุดสถานีหลักฐานดังตารางด้านล่าง วัดมุมโดยที่เล็งไปที่สถานี 1 ตั้งมุมราบเป็นศูนย์แล้ววัดมุมไปที่สถานี 2, 4, 5 และ 6 ทำการวัดระยะทางไปที่สถานี 1, 3  และ 5  รังวัดภาคทิศ (วัดดาว) ที่สถานี 1, 3 และ 6 ดังตารางด้านล่าง

คำนวณหาค่าพิกัดจุดตั้งกล้องโดยประมาณ

ผมใช้เครื่องคิดเลข fx 9860 GII SD คำนวณหาค่าพิกัดของ P  ด้วยโปรแกรม Resection ป้อนค่าพิกัดและมุมได้ดังนี้

จุด A, B และ C  ของเครื่องคิดเลขก็คือจุดที่ 1, 2 และ 6 ของโจทย์ด้านบนตามลำดับ

ผมจะใช้ค่าจุด P ที่ค่าพิกัดโดยประมาณ N = 193939.897 E = 110879.464 เพื่อเป็นค่าเริ่มต้นในการคำนวณด้วย least squares

เนื่องจากบทความยาวมาก ติดตามกันตอนที่ 2 ต่อนะครับ

Update: โปรแกรมคำนวณเล็งสกัดย้อนฉบับปรับปรุง (Resection) สำหรับเครื่องคิดเลขเทพ Casio fx 9860G II SD

Update: โปรแกรมคำนวณเล็งสกัดย้อนฉบับปรับปรุง (Resection) สำหรับเครื่องคิดเลขเทพ Casio fx 9860G II SD

ผมเขียนเรื่องการคำนวณเล็งสกัดย้อน (Resection) จากตอนก่อนหน้านี้ด้วยอัลกอริทีมใหม่ของ Josep M. Font-Llagunes อ่านได้ที่

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 5 โปรแกรมคำนวณ Resection ด้วยอัลกอริทึ่มสมัยใหม่


ในตอนนี้ผมจะมาเรียบเรียงโปรแกรมด้วยไลบรารี MyLib ที่ใช้ก่อนหน้านี้มาหลายโปรแกรมแล้ว โดยในครั้งนี้จะใช้อัลกอริทึมหลายๆตัวมาคำนวณเพื่อเปรียบเทียบกันดู เมื่อคำนวณแล้วมีรูปกราฟฟิคหยาบๆให้ดูพอกล้อมแกล้ม

การคำนวณเล็งสกัดย้อน หรือ Resection Problem มีมานานแล้วหลักๆแล้ววิธีการคำนวณแบ่งเป็น 2 วิธีหลักๆคือ แบบการตรีโกณมิติ (Trigonometric Solution) และวิธีการแบบเรขาคณิตวงกลมตัดกัน (Geometric Circle Intersection) วิธีการนี้นอกจากนำมาใช้ในงานสำรวจแล้ว ในแวดวงหุ่นยนต์ที่ใช้ในสถานที่ในร่ม ผู้อ่านอาจจะคุ้นกับการแข่งขันหุ่นยนต์ของน้องๆนักศึกษา การหาตำแหน่งในร่มของหุ่นยนต์จะอาศัยวิธีการนี้ จะอาศัยตั้งสถานีฐานที่รู้ค่าพิกัด 3 สถานี (Beacon) แล้วหุ่นยนต์จะมีอุปกรณ์จำพวก Goniometer ที่คอยหมุนกวาดหามุม 3 มุม จากตำแหน่งที่หุ่นยนต์ตั้งอยู่กับสถานีฐานนั้น จากนั้นจะคำนวณหาพิกัดตัวเองแบบ realtime

ดังนั้นการคำนวณแบบเล็งสกัดย้อนก็ยังมีคนคิดค้นอัลกอริทีมที่เร็วที่สุดเพื่อช่วยในการคำนวณให้ออกผลลัพธ์มาเร็วที่สุด ลองไปดูอัลกอริทึม Total ที่เคลมว่าเร็วที่สุด แต่สำหรับเราแล้วในแวดวงสำรวจ ความเร็วในการคำนวณไม่ใช่ประเด็นใหญ่

รายชื่อโปรแกรมเล็งสกัดย้อน

1.Kaestner-Burkhardt Method ใช้ตรีโกณมิติมาช่วยในด้านการคำนวณ วิธีการคำนวณนี้คิดค้นกันมานมนานแล้วและวิธีการคำนวณจะตั้งชื่อตามผู้ที่คิดค้น วิธีการนี้ได้มีคนนำมาปรับปรุงกันหลายครั้ง

2.TienStra Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

 

3.Collins Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

4.Cassini Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

5.Font-Llagunes ใช้วิธีเรขาคณิตโดยใช้วงกลมมาตัดกันมาช่วยในการคำนวณ

ผมประยุกต์โดยการนำสูตรมาแปลงเป็นโค้ดภาษาซี  โดยได้นำวิธีการคำนวณทั้งหมด 5 วิธี แต่คงไม่สามารถอธิบายวิธีการคำนวณได้หมด

ดาวน์โหลดโปรแกรม (Download)

ไปที่หน้าดาวน์โหลด (Download) หรือไปที่โปรแกรม Resectionจะได้ไฟล์ชื่อ “RESCTNEX.G1A” จากนั้นทำการก็อปปี้ไฟล์เข้าเครื่องคิดเลขผ่านโปรแกรม Casio FA-124 หรือ copy ผ่าน SD Card ก็สะดวกเหมือนกัน

เริ่มใช้งานโปรแกรม

ที่เครื่องคิดเลขกดปุ่ม “MENU” จะเห็นลิสต์ของไอคอนของโปรแกรมต่างๆ ตัวโปรแกรม Resection อยู่ด้านล่างๆใช้ปุ่มลูกศร เลื่อนลงมาแล้วกดคีย์ “EXE”

เมนูหลัก  (Main Menu)

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

Met – Method (F1) เลือกวิธีการคำนวณที่ผมกล่าวมาข้างต้นจำนวน 5 วิธี

Coor – Coordinates (F2) ป้อนค่าพิกัดของเป้าหลัง (control point) จุดที่ทราบค่าพิกัด 3 จุด (A, B และ C ตามลำดับ)

Ang – Angle (F3) ป้อนค่ามุมสองมุม (α และ β ตามลำดับ)

Calc – Calculate (F4) คำนวณเล็งสกัดย้อนหาค่าพิกัดจุดตั้งกล้อง

Info – Information (F5) แสดงเครดิตไลบรารีที่ใช้งาน

Exit – Exit (F6) ออกจากโปรแกรมนี้

เลือกวิธีการคำนวณ (Method)

ที่เมนูหลักของโปรแกรมกดคีย์ F1 – Met เลือกวิธีการคำนวณ จะเห็นหน้าจอแสดง Resection Method

กดคีย์ “EXE” จะเห็น drop down list ลงมาให้เลือกดังรูปด้านล่าง ในที่นี้ลองเลือก “TienStra”

กด F6 – OK เพื่อออก

ตัวอย่างที่ 1 กำหนดค่าพิกัดจุด A (1000E, 2000N) ค่าพิกัดจุด  B (1078.331E, 2077.869N) และค่าพิกัดจุด C (1172.191E, 2154.753N) วัดมุม α = 40°8’24” และมุม β = 57°36’0″ หาค่าพิกัดจุด P

ป้อนค่าพิกัดหมุดหลักฐาน (Input Control Points Coordinates)

ที่เมนูหลักกดคีย์ F2-Coor เพื่อป้อนค่าพิกัดหมุดหลักฐานที่ทราบค่าพิกัด 3 จุด ดังนี้

 

กดคีย์ F6 – OK เพื่อออก

ป้อนค่ามุม (Input Angles)

ที่เมนูหลักกดคีย์ F3 – Ang เพื่อป้อนมุม 2 มุมดังนี้ สำหรับการป้อนมุมให้ใช้เครื่องหมายลบ (-) คั่น

กดคีย์ F6 – OK เพื่อออก

คำนวณหาค่าพิกัด (Calculate Resection)

ที่เมนูหลักกดคีย์ F4 – Calc เพื่อคำนวณเล็งสกัดย้อน จะได้ผลลัพธ์ดังรูปด้านล่าง กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์ โปรแกรมจะแสดงค่าพิกัดหมุดหลักฐาน 3 จุด พร้อมทั้งมุมที่ป้อนไปเพื่อทบทวนว่าผู้ใช้ป้อนเข้าไปถูกต้องตามที่ต้องการหรือไม่ และสุดท้ายก็แสดงค่าพิกัดจุดเล็งสกัดย้อน (Resection) ที่เราต้องการคือ (1167.446E, 2016.916N)

แสดงแผนผัง (Plot)

ต้องการดูแผนผังการจัดวางตัวของจุดหลักฐาน 3 จุดและผลลัพธ์ก็สามารถดูได้คร่าวๆ กดคีย์ F5 – Plot 

สามารถกดคีย์ F1 – Dn (Down) เพื่อเลื่อนภาพลง) F2 – Up (เลื่อนภาพขึ้น) F3 – Lt (Left เลื่อนภาพไปด้านซ้าย) F4 – Rt (เลื่อนภาพไปด้านขวา) หรือกด F5 – M-> (เลือกเมนูด้านขวา) เมนูด้านขว่าจะมี F1 -Z+ (เพื่อขยายภาพ) F2 – Z- (เพื่อย่อภาพ) เสร็จแล้วกดคีย์ F6 -Done เพื่อออก กดคีย์ F6 – Done อีกครั้งเพื่อออกมาเมนูหลัก

คำนวณด้วยวิธีการอื่น

การคำนวณด้วยวิธี Tienstra ได้นำเสนอไปแล้วลองมาใช้วิธีการอื่นมาลองคำนวณบ้าง ถ้าผู้ใช้ออกจากโปรแกรมด้วยการกดคีย์ F6 – Exit โปรแกรมจะทำการเก็บค่าต่างๆที่ป้อนไปทังหมดไว้ เมื่อเข้ามาใช้โปรแกรมอีกรอบ ถ้าต้องการใช้ค่าเดิมจะไม่ต้องป้อนใหม่อีกครั้ง

Kaestner-Burkhardt Method

ที่เมนูหลักกดคีย์ F1 -Met เพื่อเลือกวิธีการคำนวณ ในที่นี้เลือก Kaestner กดคีย์ F6 – OK เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ (ไม่ต้องป้อนค่าพิกัดและค่ามุมใหม่) จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

สังเกตว่าค่าพิกัดที่ได้ (1168.228E, 2014.113N) มีความต่างจากวิธี Tienstra (1167.446E, 2016.916N)ไปพอสมควร ผมไม่แน่ใจเมื่อตรวจสอบโค้ดที่แปลงจากสูตรก็ยังคิดว่าไม่มีที่ผิด ลองไปเทียบดูจากวิธีอื่นอีกที

 

Cassini Method

ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Cassini” กดคีย์ F6 – OK  เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ (ไม่ต้องป้อนค่าพิกัดและค่ามุมใหม่) จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

 

จะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Cassini (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra

Collins Method

ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Collins” กดคีย์ F6 – OK  เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ จะได้ผลลัพธ์มาดังนี้

จะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Collins (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini

 

 

Font-Llagunes

เป็นวิธีการที่ผมชอบมากที่สุดและเคยเขียนลงแบบละเอียดในตอนก่อนหน้านี้ การใช้งานทำได้ตามขั้นตอนดังต่อไปนี้ ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Font-Llagunes” กดคีย์ F6 –  OK เพื่อออก

กดคีย์ F3 – Calc เพื่อทำการคำนวณ จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

เช่นเดียวกันจะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Font-Llagunes (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini และ  Collins แต่บางครั้งจะต่างกับวีธีการของ Kaestner-Burkhardt ในระดับหลักเมตร

 

 

ตัวอย่างที่ 2 กำหนดค่าพิกัดจุด A (1000E, 5300N) ค่าพิกัดจุด  B (2200E, 6300N) และค่าพิกัดจุด C (3100E, 5000N) วัดมุม α = 109°18’16.2″ และมุม β = 115°3’7.2″ หาค่าพิกัดจุด P

การคำนวณผมจะใช้การรวบรัดขอเสนอแต่วิธีการของ Font-Llagunes เท่านั้น ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก Font-Llangunes จาก drop down list กดคีย์ F6 – OK เพื่อออก

ที่เมนูหลักกดคีย์ F2- Coor เพื่อป้อนค่าพิกัด เสร็จแล้วกดคีย์ F6 – OK

ที่เมนูหลักกดคีย์ F3 – Ang เพื่อป้อนมุม เสร็จแล้วกดคีย์ F6 – OK

ที่เมนูหลักกดคีย์ F4 – Calc เพื่อทำคำนวณ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

จะได้ค่าพิกัดจุดเล็งสกัดย้อนมาเท่ากับ (2128.989E, 5575.375N) กดคีย์ F5 – Plot เพื่อดูแผนผัง

ก็ได้รูปการวางตัวของจุดหลักฐาน 3 จุดและจุดที่ต้องการทราบค่าคือจุด P น่าเสียดายที่ผมไม่สามารถวาดมุมแสดง dimension เพราะข้อจำกัดของเครื่องคิดเลข และต้องเตือนกันอีกนิดว่าจุด 4 จุดนี้ต้องไม่อยู่บนวงกลมเดียวกัน (โอกาสเกิดน้อยมากๆครับ) ไม่ว่าสูตรไหนก็ไม่สามารถคำนวณได้ เรียกภาวะอย่างนี้ว่าเกิดเหตุการณ์สภาวะเอกฐาน Singularity

 

สำหรับโปรแกรมในชุดเครื่องคิดเลข Casio fx-9860G II SD คงจะมีอีกหลายตอนครับ สำหรับเครดิตของโปรแกรมในชุดนี้ก็กดคีย์ F5 – Info จากเมนูหลัก ก็ติดตามกันตอนต่อไปครับ 

 

การออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection) ตอนที่ 2 (กรณีศึกษาออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำสำหรับกรุงเทพมหานครและปริมณฑล)

ผมทิ้งช่วงเรื่องการออกแบบและประยุกต์ใช้เส้นโครงแผนที่ความเพี้ยนต่ำเป็นระยะเวลาเนิ่นนานพอสมควรเนื่องจากติดภารกิจไปทำงานต่างประเทศที่หาเวลาว่างนานๆได้ยาก ถ้าผู้อ่านไม่ได้ติดตามเรื่องนี้ตั้งแต่ต้นขอให้กลับไปอ่านหรือศึกษาได้ตามลิ๊งค์ตังต่อไปนี้

แนะนำการใช้เส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection)

และ

การออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection) ตอนที่ 1

เรื่องเส้นโครงแผนที่ความเพี้ยนต่ำเป็นเรื่องใหม่สำหรับประเทศไทย แต่ในต่างประเทศบางประเทศได้ประยุกต์ใช้งานมานานแล้ว ประโยชน์ของเส้นโครงแผนที่ความเพี้ยนต่ำเมื่อประยุกต์ใช้แล้วคือ ความต่างระหว่าง Ground Distance และ Grid Distance จะน้อยมากจนสามารถละเลยไปได้ ไม่เหมือนกับการใช้แผนที่ระบบพิกัดยูทีเอ็ม (UTM) ที่ค่าระยะทางบนพื้นโลกกับระยะทางบนแผนที่ต่างกันมาก (ตัวอย่างระยะทางประมาณ 1 กม. สองระยะทางนี้อาจจะต่างกันประมาณ 40-80 ซม.แล้วแต่พื้นที่) แต่ข้อเสียคือจะต้องมีการกำหนดใช้เส้นโครงแผนที่ความเพี้ยนต่ำแบ่งเป็นพื้นที่หรือเป็นโซน ที่ค่าพิกัดศูนย์กำเนิดจะต่างกันไป อาจจะทำให้ช่างสำรวจหรือผู้ใช้งานสับสนได้ แต่ข้อเสียนี้สามารถลดลงได้ถ้ารัฐหรือหน่วยงานของรัฐได้กำหนดและประกาศใช้เป็นทางการ โดยที่มีเอกสารและไฟล์ projection สำหรับแปลงพิกัดจากระบบพิกัด UTM ไปยังระบบพิกัดที่ใช้ LDP ในแต่ละโซน ผู้ใช้งานสามารถนำค่าพารามิเตอร์นี้หรือนำไฟล์ projection (ตัวอย่างเช่นไฟล์ prj ของ Shape file) ไปแปลงพิกัดได้บนโปรแกรมด้าน GIS หรือนำไปตั้งค่าบนเครื่องมืออุปกรณ์เช่น GNSS RTK ที่สามารถแปลงพิกัดได้แบบ real time

ตัวอย่างการประยุกต์ใช้งาน

ผมขอยกตัวอย่างอีกครั้งเช่นรัฐโอเรกอนของอเมริกาที่มีการออกแบบ LDP และประกาศใช้กันมานานแล้วดังรูปด้านล่าง 

พื้นที่รัฐโอเรกอน ประมาณครึ่งหนึ่งของประเทศไทย (ประมาณ 255,000 ตร.กม.)

ออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำสำหรับกรุงเทพมหานครและปริมณฑล

ก็เป็นกรณีศึกษาก็แล้วกันนะครับ ผมจะออกแบบคร่าวๆให้พอมองเห็นภาพในภาพรวม ผมจะไล่ไปตามขั้นตอนที่ได้กล่าวไว้ในตอนที่ 1 และผมจะตั้งเป้าว่า ความเพี้ยน (Distortion) ไม่เกิน 20 ppm ก็มาดูกันว่าในพื้นที่ศึกษานี้ ค่าความเพี้ยนจะอยู่ในเกณฑ์นี้ไหม 20 ppm ก็คือระยะทางจริงๆบนพื้นโลก (Ground Distance)  1 กม. ระยะทางบนระนาบเส้นโครงแผนที่ LDP (Grid Distance) จะต่างกันไม่เกิน 20 มม. (20 มม. ต่อ 1 ล้านมิลมิเมตร หรือ 1 กม. นั่นเอง)

1.กำหนดพื้นที่ขอบเขตและหาค่าตัวแทนความสูงเฉลี่ยเหนือทรงรี (h0)

สำหรับขอบเขตก็ตามหัวข้อคือประกอบไปด้วยจังหวัดกรุงเทพมหานคร สมุทรปราการ นนทบุรี และปทุมธานี ขนาดพื้นที่ประมาณ 85 กม.ในแนวเหนือใต้ และกว้างประมาณ 75 กม. ในแนวตะวันออกตะวันตก หรือกล่าวโดยย่อพื้นที่ 85 กม. x 75 กม.

ต่อไปจะหาค่าระดับที่เป็นตัวแทนความสูงเฉลี่ยเหนือทรงรี (h0) ข้อมูลที่จะนำมาในการหาค่าเฉลี่ยจะใช้แผนที่ของกรมแผนที่ทหาร ปี 2553 ชื่อ “แผนที่แสดงค่าหมุดระดับในเขตกรุงเทพมหานครและปริมณฑล” เนื่องจากแผนที่ไม่สามารถหาแหล่งดาวน์โหลดทางการได้ จึงได้ดาวน์โหลดจากกระดานสนทนาจากเว็บไซต์ ที่ความคมชัดน้อย บางครั้งตัวเลขค่าระดับอาจจะแตกต่างค่าจริงไปบ้าง แต่ผมคิดว่าคงไม่ได้ทำให้การออกแบบ LDP กรณีศึกษานี้มีความด้อยลง  ผมนำแผนที่ชุดนี้มา ทำ rubber sheet เพื่อขึงพิกัดให้เข้ากับเส้นโครงแผนที่ UTM จากนั้นทำการ digitize จุดแต่ละจุดระดับลง ไม่ได้เอาทุกจุด แต่เลือกจุดประมาณ 10 กม.ต่อหนึ่งจุด ค่าระดับนี้เป็นค่าระดับน้ำทะเลปานกลาง (Orthometric Height) ซึ่งเราจะแปลงค่าระดับนี้ไปเป็นค่าระดับเหนือทรงรี (Ellipsoid Height) ในขั้นตอนต่อไป

จากนั้นทำการจัดเก็บจุดค่าระดับเป็นไฟล์ shape file กำหนดระบบพิกัดเป็นภูมิศาสตร์ (Geographic) เพื่อสะดวกต่อการใช้ค่าพิกัดนี้ในภายหลัง

นำไฟล์รูปที่ขึงแล้วและไฟล์จุดค่าระดับเข้าโปรแกรม QGIS ใช้ฟังก์ชั่น vector ทำการหา Basic Statistics for fields จำนวนจุดทั้งหมด 365 จุด ค่าระดับต่ำสุด 0.000 เมตร ค่าระดับสูงที่สุด  9.956 เมตร ค่าเฉลี่ย Mean 2.988 เมตร ผมจะนำค่าเฉลี่ยนี้ไปใช้งาน ค่านี้ขอใช้ตัวย่อเป็น H0 = 2.988 เมตร

ค่าระดับ H0 = 2.988 เมตร นี้จะนำมาแปลงเป็นความสูงเทียบกับทรงรี (h0) การประยุกต์ใช้ LDP ก็คือนำระนาบมาวางแตะค่าระดับนี้ โดยที่กำหนดโซนความกว้างทางราบ และช่วงค่าระดับความสูงที่ยังสามารถใช้ได้

ไดอะแกรมแสดงเส้นโครงแผนที่ความเพี้ยนต่ำที่ระนาบพิกัดฉากสัมผัสที่ความสูงเฉลี่ย h0

2.เลือกเส้นโครงแผนที่และกำหนด Central Meridian ที่จุดใกล้จุดศูนย์กลางพี้นที่

เส้นโครงแผนที่ที่นิยมนำมาทำ LDP มี 3 ประเภทคือ Transverse Mercator (TM), Lambert Conformal Conic (LCC 1SP) และ Oblique Mercator (OM) โดยที่แนวทางการเลือกถ้าพื้นที่ยาวไปในทิศทางตะวันออกตะวันตกเลือก LCC ถ้าพื้นที่ยาวไปในทิศทางเหนือใต้เลือก TM ถ้าพื้นที่เอียงไปในแนวทะแยงมุมกันทิศเหนือใต้และตะวันออกตะวันตกให้เลือก OM ในเคสนี้พื้นที่ยาวในทิศทางเหนือใต้ก็เลือกเป็น TM ที่เราคุ้นเคยกันดี

ต่อไปกำหนด Central Meridian (CM) ที่จุดกึ่งกลางของพื้นที่ (Centroid) ผมไปดาวน์โหลดไฟล์ shape file ที่รวมเอาเส้นขอบเขตของจังหวัดในประเทศไทย เมื่อโหลดมาแล้วเปิดใน QGIS จากนั้นทำการรวมพื้นที่ 4 จังหวัดนี้แบบ Dissolve เพื่อให้เหลือเส้น polygon เส้นเดียว แล้วจะนำไปหาจุดศูนย์กลางพี้นที่ โดยใช้ฟังก์ชั่นด้าน vector ของ Geometry tools เพื่อหา centroid ได้จุดมาดังรูปด้านล่าง

ได้ค่าพิกัดภูมิศาสตร์ของจุด Centroid ดังนี้ latitude: 13.852166 longitude: 100.629706 หน่วยเป็นดีกรี แปลงเป็นหน่วย DMS ได้ latitude: 31°51’7.8″N longitude: 100°37’46.94″E การวาง central meridian จะนิยมเลือกลิปดา (second) ที่เป็นจำนวนเต็ม ผมเลือก ค่าเต็มๆคือ Latitude = 31°51′ และ  CM =  100°38′

3.คำนวณหาค่าสเกลแฟคเตอร์ k0 ที่แกน Central Meridian

ก่อนหน้านี้ในข้อ 1. ผมได้ค่าระดับน้ำทะเลปานกลางของพื้นที่เฉลี่ย (H0) 2.988 เมตร จะนำค่านี้ไปแปลงเป็นค่าระดับเทียบทรงรี (h0) เตือนกันนิดว่าทรงรีที่เราใช้เป็น WGS84 การหา k0 ไม่ได้ยากตามสูตรนี้

Formula 1: Calculate Axis Scale Factor

การหาค่า h0 ก็ไม่ได้ยุ่งยากอะไรในทูลส์ Surveyor Pocket Tools ก็มีโปรแกรม Geoid Height ให้ใช้งาน h0 = H0 + N  โดยที่ N คือ Geoid Separation แต่จะก่อนคำนวณทีละขั้นตอนแบบนี้แบบแมนวลผมจะขอเสนอวิธีที่สะดวกกว่านั้น ผมจะใช้ทูลส์ชือ Init Design LDP ที่อยู่ใน Surveyor Pocket Tools มาช่วย

 คำนวณหาค่า k0 ด้วย Init Design LDP

เปิดโปรแกรม Init Design LDP จะเห็นหน้าตาโปรแกรมดังรูปด้านล่าง

ป้อนข้อมูลค่าระดับเฉลี่ย H0 = 2.988 เมตร ป้อนค่า Latitude  of project center =  13°51′ และ  Longitude of project center (CM )=  100°38′ ที่ได้จากข้อ 2.

จากนั้นคลิกที่ไอคอนลูกศรชี้ลง (เลข 3) เพื่อทำการคำนวณจะได้ผลลัพธ์ดังนี้

โปรแกรมจะคำนวณหาค่า h0 ให้และนำค่านี้ไปแทนในสมการด้านบน สุดท้ายจะคำนวณหา k0 = 0.999996 (แนะนำให้ใช้ทศนิยม 6 ตำแหน่ง) ผมพยายามขยับ CM ไปทางตะวันออกและตะวันตกครั้ง 15″ แต่ค่า k0 ยังเกาะที่ค่า 0.999996 นี้ ผมเลยเลือก CM = 100°38′

4.ตรวจสอบความเพี้ยน (Distortion) ตลอดทั้งพื้นที่

เป็นขั้นตอนที่สำคัญมาก คือถ้าเราเลือก Central Meridian มาหลายๆอันจะต้องเอาค่า k0 มาคำนวณหา Distortion

ค่า k คือ grid scale factor ค่า k นี้เราสามารถหาได้จาก สูตรด้านล่าง (เครดิตจาก Map Projection – A Working Manual ของ John P. Snyder หน้า 61)

Formula 2: Calculate Grid Scale Factor

สูตรก็เป็นสูตรเดียวที่เราใช้หา grid scale factor สำหรับ UTM เพียงแต่ค่า k0 ที่เราใช้จะเป็นค่า k0 ที่ได้จากการคำนวณด้วยโปรแกรม “Init Design LDP” คือ k0 = 0.99996, แทนค่า λ0 ด้วย 1.756383004 เรเดียน (λ0 คือ Central Meridian = 100°38′)

เลือกจุดทดสอบ

ผมเลือกจุดมาทั้งหมด 10 จุด ตำแหน่งให้อยู่ขอบๆ เป็นที่ทราบกันดี ว่าถ้าเป็นเส้นโครงแผนที่ TM ตัว grid scale factor จะเปลี่ยนจากด้านตะวันออก-ตะวันตกเท่านั้น (ไม่มีผลกับเลื่อนไปทางเหนือ-ใต้)

เรากำหนด CM ค่อนข้างจะกลางของพื้นที่ ดังนั้นจะมาดูกันว่าด้านขอบนั้นมี distortion จะยังอยู่ในเกณฑ์ไหม

การตรวจสอบจุดหาความเพี้ยนจะอาศัยการคำนวณจากสูตรที่ผมลงมาให้ก่อนหน้านี้ ลองมา workshop กันดูสักนิด  เริ่มจากจุดที่ 1  แต่จุดต่อๆไปผมจะใช้ทูลส์อีกตัวในชุด Surveyor Pocket Tools มาช่วย

Point No 1 Lat (ɸ)= 13.5079570 Long (λ) = 100.8674784 H = 6.394

คำนวณหาความสูงเทียบกับทรงรี

ได้ค่าh = H + N แทนค่าในสูตร h = 6.394 -29.9632 = -23.5692 เมตร

คำนวณ RG

RGเป็นค่าที่ขึ้นอยู่กับ latitude และพารามิเตอร์ของทรงรี (Ellipsoid) WGS84 ดังนี้

a = 6378137, e = 0.08181919084262149, e’ = 0.08209443794969568

แทนค่า